SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
fem_tn_navier_mhd.f90
Go to the documentation of this file.
1 !
2 !Authors: Jean-Luc Guermond Copyrights 2005
3 !
5  IMPLICIT NONE
6 
7 CONTAINS
8 
9  !DCQ, compute L1_norm using only list_mode(mode_idx) mode
10  FUNCTION norme_l1_one_mode(mesh, mode_idx, v) RESULT(norm)
11  USE def_type_mesh
12  USE fem_tn_axi
13  IMPLICIT NONE
14  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
15  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
16  INTEGER :: mode_idx
17  REAL(KIND=8) :: err1, s1, s2,norm
18  INTEGER :: nn
19 
20  err1 = 0.d0
21  s2 = 0.d0
22  IF (SIZE(v,2)==mesh%np) THEN
23  DO nn = 1,SIZE(v,1)
24  CALL ns_l1(mesh , abs(v(nn,:,mode_idx)), s1)
25  s2 = s2 + s1
26  END DO
27  err1 = err1 + s2
28  ! CN-AR Tue Jan 13 2009
29  ELSE
30  DO nn = 1,SIZE(v,2)
31  CALL ns_l1(mesh , abs(v(:,nn,mode_idx)), s1)
32  s2 = s2 + s1
33  END DO
34  ! CN-AR Tue Jan 13 2009
35  ! JLG/CN correction du bug CN-AR April 7, 2010
36  err1 = err1 + s2
37  ! CN-AR Tue Jan 13 2009
38  END IF
39  norm = err1
40  END FUNCTION norme_l1_one_mode
41 
42  FUNCTION norme_l2_champ(mesh, list_mode, v) RESULT(norm)
43  USE def_type_mesh
44  USE fem_tn_axi
45  IMPLICIT NONE
46  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
47  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
48  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
49  REAL(KIND=8) :: err1, s1, s2, norm
50  INTEGER :: k, nn
51 
52  err1 = 0.d0
53  IF (SIZE(v,2)==mesh%np) THEN
54  DO k = 1, SIZE(list_mode)
55  s2 = 0.d0
56  DO nn = 1,SIZE(v,1)
57  s1 = 0.d0
58  CALL ns_0(mesh , v(nn,:,k), s1)
59  s2 = s2 + s1**2
60  END DO
61  ! CN-AR Tue Jan 13 2009
62  ! JLG/CN correction du bug CN-AR April 7, 2010
63  IF (list_mode(k) /= 0) THEN
64  s2 = s2/2.d0
65  END IF
66  err1 = err1 + s2
67  ! CN-AR Tue Jan 13 2009
68 
69  ENDDO
70  ELSE
71  DO k = 1, SIZE(list_mode)
72  s2 = 0.d0
73  DO nn = 1,SIZE(v,2)
74  s1 = 0.d0
75  CALL ns_0(mesh , v(:,nn,k), s1)
76  s2 = s2 + s1**2
77  END DO
78  ! CN-AR Tue Jan 13 2009
79  ! JLG/CN correction du bug CN-AR April 7, 2010
80  IF (list_mode(k) /= 0) THEN
81  s2 = s2/2.d0
82  END IF
83  err1 = err1 + s2
84  ! CN-AR Tue Jan 13 2009
85 
86  ENDDO
87  END IF
88  norm = sqrt(err1)
89 
90  END FUNCTION norme_l2_champ
91 
92  FUNCTION dot_product_champ(mesh, list_mode, v, w) RESULT(norm)
93 
94  USE def_type_mesh
95  USE fem_tn_axi
96 
97  IMPLICIT NONE
98  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
99  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
100  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v, w
101 
102  REAL(KIND=8) :: err1, s1, s2, norm
103  INTEGER :: k, nn
104 
105  err1 = 0.d0
106  IF (SIZE(v,2)==mesh%np) THEN
107  DO k = 1, SIZE(list_mode)
108  s2 = 0.d0
109  DO nn = 1,SIZE(v,1)
110  s1 = 0.d0
111  CALL dot_product(mesh , v(nn,:,k), w(nn,:,k), s1)
112  s2 = s2 + s1
113  END DO
114  ! CN-AR Tue Jan 13 2009
115  ! JLG/CN correction du bug CN-AR April 7, 2010
116  IF (list_mode(k) /= 0) THEN
117  s2 = s2/2.d0
118  END IF
119  err1 = err1 + s2
120  ! CN-AR Tue Jan 13 2009
121 
122  ENDDO
123  ELSE
124  DO k = 1, SIZE(list_mode)
125  s2 = 0.d0
126  DO nn = 1,SIZE(v,2)
127  s1 = 0.d0
128  CALL dot_product(mesh , v(:,nn,k), w(:,nn,k), s1)
129  s2 = s2 + s1
130  END DO
131  ! CN-AR Tue Jan 13 2009
132  ! JLG/CN correction du bug CN-AR April 7, 2010
133  IF (list_mode(k) /= 0) THEN
134  s2 = s2/2.d0
135  END IF
136  err1 = err1 + s2
137  ! CN-AR Tue Jan 13 2009
138 
139  ENDDO
140  END IF
141  norm = err1
142 
143  END FUNCTION dot_product_champ
144 
145  FUNCTION norme_h1_champ(mesh, list_mode, v) RESULT(norm)
146 
147  USE def_type_mesh
148  USE fem_tn_axi
149 
150  IMPLICIT NONE
151  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
152  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
153  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
154 
155  REAL(KIND=8) :: err1, s1, s2, s0, norm
156  INTEGER :: k,nn
157 
158  err1 = 0.d0
159  IF (SIZE(v,2)==mesh%np) THEN
160  DO k = 1, SIZE(list_mode)
161  s2 = 0.d0
162  DO nn = 1,SIZE(v,1)
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
166  END DO
167  ! CN-AR Tue Jan 13 2009
168  ! JLG/CN correction du bug CN-AR April 7, 2010
169  IF (list_mode(k) /= 0) THEN
170  s2 = s2/2.d0
171  END IF
172  err1 = err1 + s2
173  ! CN-AR Tue Jan 13 2009
174  ENDDO
175  ELSE
176  DO k = 1, SIZE(list_mode)
177  s2=0.d0
178  DO nn = 1,SIZE(v,2)
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
182  END DO
183  ! CN-AR Tue Jan 13 2009
184  ! JLG/CN correction du bug CN-AR April 7, 2010
185  IF (list_mode(k) /= 0) THEN
186  s2 = s2/2.d0
187  END IF
188  err1 = err1 + s2
189  ! CN-AR Tue Jan 13 2009
190  ENDDO
191  END IF
192 
193  norm = sqrt(err1)
194 
195  END FUNCTION norme_h1_champ
196 
197  FUNCTION norme_div(H_mesh, list_mode, v) RESULT(norm)
198 
199  USE def_type_mesh
200  USE fem_tn_axi
201 
202  IMPLICIT NONE
203  TYPE(mesh_type), INTENT(IN) :: h_mesh !type de maillage
204  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
205  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN), TARGET :: v
206 
207 
208  INTEGER :: m_max_c, mode, k, m, l, ni, i
209  REAL(KIND=8) :: norm, err, div1, div2, jr, ray
210 
211  err = 0.d0
212 
213  m_max_c = SIZE(list_mode)
214 
215  IF (SIZE(v,2)==h_mesh%np) THEN
216 
217  DO m = 1, h_mesh%me
218  DO l = 1, h_mesh%gauss%l_G
219 
220  !===Compute radius of Gauss point
221  ray = 0.d0
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)
224  END DO
225  jr = ray * h_mesh%gauss%rj(l,m)
226 
227 
228  DO k=1, m_max_c
229  mode = list_mode(k)
230  div1 = 0.d0
231  div2 = 0.d0
232 
233  DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
234 
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)
239 
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)
244 
245  ENDDO
246 
247  err = err + (div1**2+div2**2)*jr
248  ENDDO
249 
250  END DO
251  END DO
252  ELSE
253  DO m = 1, h_mesh%me
254  DO l = 1, h_mesh%gauss%l_G
255 
256  !===Compute radius of Gauss point
257  ray = 0.d0
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)
260  END DO
261  jr = ray * h_mesh%gauss%rj(l,m)
262 
263 
264  DO k=1, m_max_c
265  mode = list_mode(k)
266  div1 = 0.d0
267  div2 = 0.d0
268 
269  DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
270 
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)
275 
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)
280 
281  ENDDO
282 
283  err = err + (div1**2+div2**2)*jr
284  ENDDO
285 
286  END DO
287  END DO
288 
289  END IF
290 
291 
292  norm = sqrt(err)
293 
294  END FUNCTION norme_div
295 
296  FUNCTION norme_curl(H_mesh, list_mode, v) RESULT(norm)
297 
298  USE def_type_mesh
299  USE fem_tn_axi
300 
301  IMPLICIT NONE
302  TYPE(mesh_type), INTENT(IN) :: h_mesh !type de maillage
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
308 
309  m_max_c = SIZE(list_mode)
310  err = 0.d0
311 
312  IF (SIZE(v,2)==h_mesh%np) THEN
313 
314  DO m = 1, h_mesh%me
315  DO l = 1, h_mesh%gauss%l_G
316 
317  !===Compute radius of Gauss point
318  ray = 0
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)
321  END DO
322  jr = ray * h_mesh%gauss%rj(l,m)
323 
324  DO k=1, m_max_c
325  mode = list_mode(k)
326  c = 0
327  DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
328  !--------Composante r------
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))
333  !--------Composante theta------
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))
338  !--------Composante z------
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))
345  ENDDO
346  err = err + sum(c**2)*jr
347  END DO
348  END DO
349  END DO
350 
351  ELSE
352  DO m = 1, h_mesh%me
353  DO l = 1, h_mesh%gauss%l_G
354 
355  !===Compute radius of Gauss point
356  ray = 0
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)
359  END DO
360  jr = ray * h_mesh%gauss%rj(l,m)
361 
362  DO k=1, m_max_c
363  mode = list_mode(k)
364  c = 0
365  DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
366  !--------Composante r------
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))
371  !--------Composante theta------
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))
376  !--------Composante z------
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))
383  ENDDO
384  err = err + sum(c**2)*jr
385  END DO
386  END DO
387  END DO
388  END IF
389 
390  norm = sqrt(err)
391 
392  END FUNCTION norme_curl
393 
394  SUBROUTINE norme_interface(H_mesh,phi_mesh,INTERFACE,mu_H_field,mu_phi,H,phi,mode,x)
395  USE def_type_mesh
396  USE dir_nodes
397  USE gauss_points
398 
399  IMPLICIT NONE
400 
401  TYPE(mesh_type), INTENT(IN) :: h_mesh, phi_mesh
402  TYPE(interface_type), INTENT(IN) :: interface
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
411 
412 
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
416  CALL gauss(phi_mesh)
417 
418  IF (h_mesh%gauss%n_ws == n_ws) THEN
419  w_cs = wws
420  ELSE
421  DO ls = 1, l_gs
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)
424  w_cs(3,ls)=0
425  END DO
426  END IF
427 
428  n_ws1 = h_mesh%gauss%n_ws
429  n_ws2 = phi_mesh%gauss%n_ws
430 
431  err = 0
432  x = 0
433  z = 0
434  zmu = 0
435 
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)
440 
441  DO ls = 1, l_gs
442 
443  !===Compute radius of Gauss point
444  ray = 0.d0
445  DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
446  ray = ray + phi_mesh%rr(1,i)* wws(ni,ls)
447  END DO
448 
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))
451 
452  DO i = 1, 6
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))
455  ENDDO
456 
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))
463 
464  z = z + sum(b(:)**2)*rjs(ls,ms2)*ray
465  zmu = zmu + sum(mub(:)**2)*rjs(ls,ms2)*ray
466 
467  !Error on tangential component (magnetic induction)
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)
474 
475  !Error on normal component (magnetic field)
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)
481 
482  END DO
483  END DO
484  ELSE
485  DO ms = 1, interface%mes
486  ms2 = interface%mesh2(ms)
487  m = phi_mesh%neighs(ms2)
488 
489  DO ls = 1, l_gs
490  !June 6 2008, muhl
491  muhl = sum(mu_h_field(interface%jjs1(1:n_ws1,ms))*w_cs(1:n_ws1,ls))
492  !June 6 2008, muhl
493  !===Compute radius of Gauss point
494  ray = 0.d0
495  DO ni = 1, n_ws2; i = phi_mesh%jjs(ni,ms2)
496  ray = ray + phi_mesh%rr(1,i)* wws(ni,ls)
497  END DO
498 
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))
501 
502  DO i = 1, 6
503  b(i) = sum(h(interface%jjs1(1:n_ws1,ms),i)*w_cs(1:n_ws1,ls))
504  mub(i) = b(i)*muhl
505  ENDDO
506 
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))
513 
514  z = z + sum(b(:)**2)*rjs(ls,ms2)*ray
515  zmu = zmu + sum(mub(:)**2)*rjs(ls,ms2)*ray
516 
517  !Error on tangential component (magnetic induction)
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)
524 
525  !Error on normal component (magnetic field)
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)
531 
532  END DO
533  END DO
534  END IF
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)
538 
539  END SUBROUTINE norme_interface
540 
541  SUBROUTINE norme_interface_h_mu(H_mesh,INTERFACE,mu_H_field,H,x)
542  USE def_type_mesh
543  USE dir_nodes
544  USE gauss_points
545 
546  IMPLICIT NONE
547 
548  TYPE(mesh_type), INTENT(IN) :: h_mesh
549  TYPE(interface_type), INTENT(IN) :: interface
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
555 
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
559  CALL gauss(h_mesh)
560 
561  w_1s = wws
562  w_2s = wws
563 
564  n_ws1 = h_mesh%gauss%n_ws
565  n_ws2 = h_mesh%gauss%n_ws
566 
567  err = 0
568  x = 0
569  z = 0
570  zmu = 0
571 
572  DO ms = 1, interface%mes
573  ms2 = interface%mesh2(ms)
574  m = h_mesh%neighs(ms2)
575 
576  DO ls = 1, l_gs
577 
578  !===Compute radius of Gauss point
579  ray = 0.d0
580  DO ni = 1, n_ws2; i = h_mesh%jjs(ni,ms2)
581  ray = ray + h_mesh%rr(1,i)* wws(ni,ls)
582  END DO
583 
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))
586 
587  DO i = 1, 6
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))
592  ENDDO
593 
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
596 
597  !Error on tangential component (magnetic field)
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)
604  !Error on normal component (magnetic induction)
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)
610 
611  END DO
612  END DO
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)
616 
617  END SUBROUTINE norme_interface_h_mu
618 
619 END MODULE fem_tn_ns_mhd
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)
Definition: fem_tn_axi.f90:248
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)
Definition: fem_tn_axi.f90:51
real(kind=8) function norme_curl(H_mesh, list_mode, v)
subroutine dot_product(mesh, ff, gg, t)
Definition: fem_tn_axi.f90:8
subroutine gauss(mesh)
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)
Definition: fem_tn_axi.f90:95
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
Definition: doc_intro.h:193