SFEMaNS  version 5.3
Reference documentation for SFEMaNS
fem_tn_navier_mhd.f90
Go to the documentation of this file.
1 !
2 !Authors: Jean-Luc Guermond Copyrights 2005
3 !
4 MODULE fem_tn_ns_mhd
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
def_type_mesh::interface_type
Definition: def_type_mesh.f90:116
fem_tn_axi
Definition: fem_tn_axi.f90:4
gauss_points
Definition: associate_gauss.f90:4
fem_tn_ns_mhd::dot_product_champ
real(kind=8) function dot_product_champ(mesh, list_mode, v, w)
Definition: fem_tn_navier_mhd.f90:93
gauss_points::l_gs
integer, public l_gs
Definition: associate_gauss.f90:9
gauss_points::n_ws
integer, public n_ws
Definition: associate_gauss.f90:9
fem_tn_ns_mhd::norme_l2_champ
real(kind=8) function norme_l2_champ(mesh, list_mode, v)
Definition: fem_tn_navier_mhd.f90:43
fem_tn_axi::ns_l1
subroutine ns_l1(mesh, ff, t)
Definition: fem_tn_axi.f90:96
fem_tn_ns_mhd::norme_l1_one_mode
real(kind=8) function norme_l1_one_mode(mesh, mode_idx, v)
Definition: fem_tn_navier_mhd.f90:11
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
gauss_points::rnorms
real(kind=8), dimension(:,:,:), pointer rnorms
Definition: associate_gauss.f90:13
fem_tn_ns_mhd::norme_interface
subroutine norme_interface(H_mesh, phi_mesh, INTERFACE, mu_H_field, mu_phi, H, phi, mode, x)
Definition: fem_tn_navier_mhd.f90:395
fem_tn_ns_mhd::norme_h1_champ
real(kind=8) function norme_h1_champ(mesh, list_mode, v)
Definition: fem_tn_navier_mhd.f90:146
fem_tn_ns_mhd::norme_div
real(kind=8) function norme_div(H_mesh, list_mode, v)
Definition: fem_tn_navier_mhd.f90:198
fem_tn_ns_mhd::norme_curl
real(kind=8) function norme_curl(H_mesh, list_mode, v)
Definition: fem_tn_navier_mhd.f90:297
fem_tn_axi::ns_1
subroutine ns_1(mesh, ff, t)
Definition: fem_tn_axi.f90:249
fem_tn_axi::dot_product
subroutine dot_product(mesh, ff, gg, t)
Definition: fem_tn_axi.f90:9
fem_tn_axi::ns_0
subroutine ns_0(mesh, ff, t)
Definition: fem_tn_axi.f90:52
fem_tn_ns_mhd
Definition: fem_tn_navier_mhd.f90:4
gauss_points::rjs
real(kind=8), dimension(:,:), pointer rjs
Definition: associate_gauss.f90:15
gauss_points::wws
real(kind=8), dimension(:,:), pointer wws
Definition: associate_gauss.f90:11
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:24
fem_tn_ns_mhd::norme_interface_h_mu
subroutine norme_interface_h_mu(H_mesh, INTERFACE, mu_H_field, H, x)
Definition: fem_tn_navier_mhd.f90:542
dir_nodes
Definition: dir_nodes.f90:4
def_type_mesh
Definition: def_type_mesh.f90:4
gauss_points::dw_s
real(kind=8), dimension(:,:,:,:), pointer dw_s
Definition: associate_gauss.f90:16