15 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
17 REAL(KIND=8) :: err1, s1, s2,norm
22 IF (
SIZE(v,2)==mesh%np)
THEN
24 CALL
ns_l1(mesh , abs(v(nn,:,mode_idx)), s1)
31 CALL
ns_l1(mesh , abs(v(:,nn,mode_idx)), s1)
47 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
48 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
49 REAL(KIND=8) :: err1, s1, s2, norm
53 IF (
SIZE(v,2)==mesh%np)
THEN
54 DO k = 1,
SIZE(list_mode)
58 CALL
ns_0(mesh , v(nn,:,k), s1)
63 IF (list_mode(k) /= 0)
THEN
71 DO k = 1,
SIZE(list_mode)
75 CALL
ns_0(mesh , v(:,nn,k), s1)
80 IF (list_mode(k) /= 0)
THEN
99 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
100 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v, w
102 REAL(KIND=8) :: err1, s1, s2, norm
106 IF (
SIZE(v,2)==mesh%np)
THEN
107 DO k = 1,
SIZE(list_mode)
116 IF (list_mode(k) /= 0)
THEN
124 DO k = 1,
SIZE(list_mode)
133 IF (list_mode(k) /= 0)
THEN
152 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
153 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
155 REAL(KIND=8) :: err1, s1, s2, s0, norm
159 IF (
SIZE(v,2)==mesh%np)
THEN
160 DO k = 1,
SIZE(list_mode)
163 CALL
ns_1(mesh , v(nn,:,k), s1)
164 CALL
ns_0(mesh , v(nn,:,k), s0)
165 s2 = s2+ s1**2 + list_mode(k)**2*s0**2
169 IF (list_mode(k) /= 0)
THEN
176 DO k = 1,
SIZE(list_mode)
179 CALL
ns_1(mesh , v(:,nn,k), s1)
180 CALL
ns_0(mesh , v(:,nn,k), s0)
181 s2 = s2 + s1**2 + list_mode(k)**2*s0**2
185 IF (list_mode(k) /= 0)
THEN
204 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
205 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN),
TARGET :: v
208 INTEGER :: m_max_c, mode, k, m, l, ni, i
209 REAL(KIND=8) :: norm, err, div1, div2, jr, ray
213 m_max_c =
SIZE(list_mode)
215 IF (
SIZE(v,2)==h_mesh%np)
THEN
218 DO l = 1, h_mesh%gauss%l_G
222 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
223 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
225 jr = ray * h_mesh%gauss%rj(l,m)
233 DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
235 div1 = div1 + v(1,i,k)*h_mesh%gauss%ww(ni,l)/ray &
236 + v(1,i,k)*h_mesh%gauss%dw(1,ni,l,m) &
237 + mode/ray*v(4,i,k)*h_mesh%gauss%ww(ni,l) &
238 + v(5,i,k)*h_mesh%gauss%dw(2,ni,l,m)
240 div2 = div2 + v(2,i,k)*h_mesh%gauss%ww(ni,l)/ray &
241 + v(2,i,k)*h_mesh%gauss%dw(1,ni,l,m) &
242 - mode/ray*v(3,i,k)*h_mesh%gauss%ww(ni,l) &
243 + v(6,i,k)*h_mesh%gauss%dw(2,ni,l,m)
247 err = err + (div1**2+div2**2)*jr
254 DO l = 1, h_mesh%gauss%l_G
258 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
259 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
261 jr = ray * h_mesh%gauss%rj(l,m)
269 DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
271 div1 = div1 + v(i,1,k)*h_mesh%gauss%ww(ni,l)/ray &
272 + v(i,1,k)*h_mesh%gauss%dw(1,ni,l,m) &
273 + mode/ray*v(i,4,k)*h_mesh%gauss%ww(ni,l) &
274 + v(i,5,k)*h_mesh%gauss%dw(2,ni,l,m)
276 div2 = div2 + v(i,2,k)*h_mesh%gauss%ww(ni,l)/ray &
277 + v(i,2,k)*h_mesh%gauss%dw(1,ni,l,m) &
278 - mode/ray*v(i,3,k)*h_mesh%gauss%ww(ni,l) &
279 + v(i,6,k)*h_mesh%gauss%dw(2,ni,l,m)
283 err = err + (div1**2+div2**2)*jr
303 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
304 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: v
305 REAL(KIND=8) :: err, norm, jr, ray
306 INTEGER :: m_max_c, mode, k, m, l, ni, i
307 REAL(KIND=8),
DIMENSION(6) :: c
309 m_max_c =
SIZE(list_mode)
312 IF (
SIZE(v,2)==h_mesh%np)
THEN
315 DO l = 1, h_mesh%gauss%l_G
319 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
320 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
322 jr = ray * h_mesh%gauss%rj(l,m)
327 DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
329 c(1) = c(1) + ( mode/ray*v(6,i,k)*h_mesh%gauss%ww(ni,l) &
330 - v(3,i,k)*h_mesh%gauss%dw(2,ni,l,m))
331 c(2) = c(2) + (-mode/ray*v(5,i,k)*h_mesh%gauss%ww(ni,l) &
332 - v(4,i,k)*h_mesh%gauss%dw(2,ni,l,m))
334 c(3) = c(3) + (v(1,i,k)*h_mesh%gauss%dw(2,ni,l,m) &
335 - v(5,i,k)*h_mesh%gauss%dw(1,ni,l,m))
336 c(4) = c(4) + (v(2,i,k)*h_mesh%gauss%dw(2,ni,l,m) &
337 - v(6,i,k)*h_mesh%gauss%dw(1,ni,l,m))
339 c(5) = c(5) + (v(3,i,k)*h_mesh%gauss%dw(1,ni,l,m) &
340 + v(3,i,k)*h_mesh%gauss%ww(ni,l)/ray &
341 - mode/ray*v(2,i,k)*h_mesh%gauss%ww(ni,l))
342 c(6) = c(6) + (v(4,i,k)*h_mesh%gauss%dw(1,ni,l,m) &
343 + v(4,i,k)*h_mesh%gauss%ww(ni,l)/ray &
344 + mode/ray*v(1,i,k)*h_mesh%gauss%ww(ni,l))
346 err = err + sum(c**2)*jr
353 DO l = 1, h_mesh%gauss%l_G
357 DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
358 ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
360 jr = ray * h_mesh%gauss%rj(l,m)
365 DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
367 c(1) = c(1) + ( mode/ray*v(i,6,k)*h_mesh%gauss%ww(ni,l) &
368 - v(i,3,k)*h_mesh%gauss%dw(2,ni,l,m))
369 c(2) = c(2) + (-mode/ray*v(i,5,k)*h_mesh%gauss%ww(ni,l) &
370 - v(i,4,k)*h_mesh%gauss%dw(2,ni,l,m))
372 c(3) = c(3) + (v(i,1,k)*h_mesh%gauss%dw(2,ni,l,m) &
373 - v(i,5,k)*h_mesh%gauss%dw(1,ni,l,m))
374 c(4) = c(4) + (v(i,2,k)*h_mesh%gauss%dw(2,ni,l,m) &
375 - v(i,6,k)*h_mesh%gauss%dw(1,ni,l,m))
377 c(5) = c(5) + (v(i,3,k)*h_mesh%gauss%dw(1,ni,l,m) &
378 + v(i,3,k)*h_mesh%gauss%ww(ni,l)/ray &
379 - mode/ray*v(i,2,k)*h_mesh%gauss%ww(ni,l))
380 c(6) = c(6) + (v(i,4,k)*h_mesh%gauss%dw(1,ni,l,m) &
381 + v(i,4,k)*h_mesh%gauss%ww(ni,l)/ray &
382 + mode/ray*v(i,1,k)*h_mesh%gauss%ww(ni,l))
384 err = err + sum(c**2)*jr
401 TYPE(mesh_type),
INTENT(IN) :: h_mesh, phi_mesh
403 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: h
404 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: phi
405 INTEGER,
INTENT(IN) :: mode
406 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
407 REAL(KIND=8),
INTENT(IN) :: mu_phi
408 REAL(KIND=8),
INTENT(OUT) :: x
409 REAL(KIND=8),
DIMENSION(2) :: rgauss
410 REAL(KIND=8),
DIMENSION(phi_mesh%gauss%n_ws,phi_mesh%gauss%l_Gs) :: w_cs
413 INTEGER :: ms, ls, ms2, n_ws1, n_ws2, m, i, ni
414 REAL(KIND=8),
DIMENSION(6) :: b, mub, grd
415 REAL(KIND=8) ::
z, zmu, ray, err, muhl
418 IF (h_mesh%gauss%n_ws == n_ws)
THEN
422 w_cs(1,ls)= wws(1,ls)+0.5*wws(3,ls)
423 w_cs(2,ls)= wws(2,ls)+0.5*wws(3,ls)
428 n_ws1 = h_mesh%gauss%n_ws
429 n_ws2 = phi_mesh%gauss%n_ws
436 IF (
SIZE(h,2)==h_mesh%np)
THEN
437 DO ms = 1, interface%mes
438 ms2 = interface%mesh2(ms)
439 m = phi_mesh%neighs(ms2)
445 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
446 ray = ray + phi_mesh%rr(1,i)* wws(ni,ls)
449 rgauss(1) = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* wws(:,ls))
450 rgauss(2) = sum(phi_mesh%rr(2,phi_mesh%jjs(:,ms2))* wws(:,ls))
453 b(i) = sum(h(i,interface%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
454 mub(i) = sum(mu_h_field(interface%jjs1(1:n_ws1,ms))*h(i,interface%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
457 grd(1) = sum(phi(1,phi_mesh%jj(:,m))*dw_s(1,:,ls,ms2))
458 grd(2) = sum(phi(2,phi_mesh%jj(:,m))*dw_s(1,:,ls,ms2))
459 grd(3) = mode/ray * sum(phi(2,interface%jjs2(:,ms))*wws(:,ls))
460 grd(4) = -mode/ray * sum(phi(1,interface%jjs2(:,ms))*wws(:,ls))
461 grd(5) = sum(phi(1,phi_mesh%jj(:,m))*dw_s(2,:,ls,ms2))
462 grd(6) = sum(phi(2,phi_mesh%jj(:,m))*dw_s(2,:,ls,ms2))
464 z =
z + sum(b(:)**2)*rjs(ls,ms2)*ray
465 zmu = zmu + sum(mub(:)**2)*rjs(ls,ms2)*ray
468 x = x + ray* rjs(ls,ms2)*( &
469 ((b(5)-grd(5))*rnorms(1,ls,ms2) - &
470 (b(1)-grd(1))*rnorms(2,ls,ms2))**2 + &
471 ((b(6)-grd(6))*rnorms(1,ls,ms2) - &
472 (b(2)-grd(2))*rnorms(2,ls,ms2))**2 + &
473 (b(3)-grd(3))**2 + (b(4)-grd(4))**2)
476 err = err + ray* rjs(ls,ms2)*( &
477 ((mub(1)-mu_phi*grd(1))*rnorms(1,ls,ms2) +&
478 (mub(5)-mu_phi*grd(5))*rnorms(2,ls,ms2))**2 + &
479 ((mub(2)-mu_phi*grd(2))*rnorms(1,ls,ms2) +&
480 (mub(6)-mu_phi*grd(6))*rnorms(2,ls,ms2))**2)
485 DO ms = 1, interface%mes
486 ms2 = interface%mesh2(ms)
487 m = phi_mesh%neighs(ms2)
491 muhl = sum(mu_h_field(interface%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
495 DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
496 ray = ray + phi_mesh%rr(1,i)* wws(ni,ls)
499 rgauss(1) = sum(phi_mesh%rr(1,phi_mesh%jjs(:,ms2))* wws(:,ls))
500 rgauss(2) = sum(phi_mesh%rr(2,phi_mesh%jjs(:,ms2))* wws(:,ls))
503 b(i) = sum(h(interface%jjs1(1:n_ws1,ms),i)*w_cs(1:n_ws1,ls))
507 grd(1) = sum(phi(phi_mesh%jj(:,m),1)*dw_s(1,:,ls,ms2))
508 grd(2) = sum(phi(phi_mesh%jj(:,m),2)*dw_s(1,:,ls,ms2))
509 grd(3) = mode/ray * sum(phi(interface%jjs2(:,ms),2)*wws(:,ls))
510 grd(4) =-mode/ray * sum(phi(interface%jjs2(:,ms),1)*wws(:,ls))
511 grd(5) = sum(phi(phi_mesh%jj(:,m),1)*dw_s(2,:,ls,ms2))
512 grd(6) = sum(phi(phi_mesh%jj(:,m),2)*dw_s(2,:,ls,ms2))
514 z =
z + sum(b(:)**2)*rjs(ls,ms2)*ray
515 zmu = zmu + sum(mub(:)**2)*rjs(ls,ms2)*ray
518 x = x + ray* rjs(ls,ms2)*( &
519 ((b(5)-grd(5))*rnorms(1,ls,ms2) - &
520 (b(1)-grd(1))*rnorms(2,ls,ms2))**2 + &
521 ((b(6)-grd(6))*rnorms(1,ls,ms2) - &
522 (b(2)-grd(2))*rnorms(2,ls,ms2))**2 + &
523 (b(3)-grd(3))**2 + (b(4)-grd(4))**2)
526 err = err + ray* rjs(ls,ms2)*( &
527 ((mub(1)-mu_phi*grd(1))*rnorms(1,ls,ms2) +&
528 (mub(5)-mu_phi*grd(5))*rnorms(2,ls,ms2))**2 + &
529 ((mub(2)-mu_phi*grd(2))*rnorms(1,ls,ms2) +&
530 (mub(6)-mu_phi*grd(6))*rnorms(2,ls,ms2))**2)
535 WRITE(*,
'(A,e12.5)')
'Collage tangent Sigma_phi ', sqrt(x)/(sqrt(
z)+1.d-16)
536 WRITE(*,
'(A,e12.5)')
'Collage normal Sigma_phi ', sqrt(err)/(sqrt(zmu)+1.d-16)
537 x = sqrt(x)/(sqrt(
z)+1.d-16) + sqrt(err)/(sqrt(zmu)+1.d-16)
550 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: h
551 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: mu_h_field
552 REAL(KIND=8),
INTENT(OUT) :: x
553 REAL(KIND=8),
DIMENSION(2) :: rgauss
554 REAL(KIND=8),
DIMENSION(H_mesh%gauss%n_ws,H_mesh%gauss%l_Gs) :: w_1s, w_2s
556 INTEGER :: ms, ls, ms2, n_ws1, n_ws2, m, i, ni
557 REAL(KIND=8),
DIMENSION(6) :: b1, mub1, b2, mub2
558 REAL(KIND=8) ::
z, zmu, ray, err
564 n_ws1 = h_mesh%gauss%n_ws
565 n_ws2 = h_mesh%gauss%n_ws
572 DO ms = 1, interface%mes
573 ms2 = interface%mesh2(ms)
574 m = h_mesh%neighs(ms2)
580 DO ni = 1, n_ws2; i = h_mesh%jjs(ni,ms2)
581 ray = ray + h_mesh%rr(1,i)* wws(ni,ls)
584 rgauss(1) = sum(h_mesh%rr(1,h_mesh%jjs(:,ms2))* wws(:,ls))
585 rgauss(2) = sum(h_mesh%rr(2,h_mesh%jjs(:,ms2))* wws(:,ls))
588 b1(i) = sum(h(interface%jjs1(1:n_ws1,ms),i)*w_1s(1:n_ws1,ls))
589 mub1(i) = sum(mu_h_field(interface%jjs1(1:n_ws1,ms))*h(interface%jjs1(1:n_ws1,ms),i)*w_1s(1:n_ws1,ls))
590 b2(i) = sum(h(interface%jjs2(1:n_ws2,ms),i)*w_2s(1:n_ws2,ls))
591 mub2(i) = sum(mu_h_field(interface%jjs2(1:n_ws2,ms))*h(interface%jjs2(1:n_ws2,ms),i)*w_2s(1:n_ws2,ls))
594 z =
z + (sum(b2(:)**2)+sum(b1(:)**2))*rjs(ls,ms2)*ray
595 zmu = zmu + (sum(mub1(:)**2)+sum(mub2(:)**2))*rjs(ls,ms2)*ray
598 x = x + ray* rjs(ls,ms2)*( &
599 ((b1(5)-b2(5))*rnorms(1,ls,ms2) - &
600 (b1(1)-b2(1))*rnorms(2,ls,ms2))**2 + &
601 ((b1(6)-b2(6))*rnorms(1,ls,ms2) - &
602 (b1(2)-b2(2))*rnorms(2,ls,ms2))**2 + &
603 (b1(3)-b2(3))**2 + (b1(4)-b2(4))**2)
605 err = err + ray* rjs(ls,ms2)*( &
606 ((mub1(1)-mub2(1))*rnorms(1,ls,ms2) +&
607 (mub1(5)-mub2(5))*rnorms(2,ls,ms2))**2 + &
608 ((mub1(2)-mub2(2))*rnorms(1,ls,ms2) +&
609 (mub1(6)-mub2(6))*rnorms(2,ls,ms2))**2)
613 WRITE(*,
'(A,e12.5)')
'Collage tangent Sigma_mu ', sqrt(x)/(sqrt(
z)+1.d-16)
614 WRITE(*,
'(A,e12.5)')
'Collage normal Sigma_mu ', sqrt(err)/(sqrt(zmu)+1.d-16)
615 x = sqrt(x)/(sqrt(
z)+1.d-16) + sqrt(err)/(sqrt(zmu)+1.d-16)
real(kind=8) function norme_h1_champ(mesh, list_mode, v)
real(kind=8) function norme_div(H_mesh, list_mode, v)
subroutine ns_1(mesh, ff, t)
real(kind=8) function dot_product_champ(mesh, list_mode, v, w)
real(kind=8) function norme_l2_champ(mesh, list_mode, v)
subroutine ns_0(mesh, ff, t)
real(kind=8) function norme_curl(H_mesh, list_mode, v)
subroutine dot_product(mesh, ff, gg, t)
subroutine norme_interface_h_mu(H_mesh, INTERFACE, mu_H_field, H, x)
subroutine norme_interface(H_mesh, phi_mesh, INTERFACE, mu_H_field, mu_phi, H, phi, mode, x)
subroutine ns_l1(mesh, ff, t)
real(kind=8) function norme_l1_one_mode(mesh, mode_idx, v)
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic z