SFEMaNS  version 4.1 (work in progress) Reference documentation for SFEMaNS
tn_axi.f90
Go to the documentation of this file.
1 !
2 !Authors Jean-Luc Guermond, Raphael Laguerre, Copyrights 2005
3 !Revised: Jean-Luc Guermond Jan 2011
4 !
5 MODULE tn_axi
6
7 CONTAINS
8
9  FUNCTION dot_product_sf(communicator, mesh, list_mode, v, w) RESULT(norm)
10  USE def_type_mesh
12  USE fem_tn_ns_mhd
13  IMPLICIT NONE
14  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
15  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
16  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v, w
17  REAL(KIND=8) :: norm_loc, norm_tot, norm
18  INTEGER :: code
19  REAL(KIND=8) :: pi
20 #include "petsc/finclude/petsc.h"
21  mpi_comm, DIMENSION(2) :: communicator
22
23  norm = 0.d0
24  norm_tot = 0.d0
25  IF (mesh%me==0) THEN
26  norm_loc = 0.d0
27  ELSE
28  norm_loc = dot_product_champ(mesh, list_mode, v, w)
29  END IF
30
31  CALL mpi_allreduce(norm_loc,norm_tot,1,mpi_double_precision, mpi_sum, communicator(2), code)
32  CALL mpi_allreduce(norm_tot,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
33  pi = acos(-1.d0)
34  norm = norm*2*pi
35
36  END FUNCTION dot_product_sf
37
38  FUNCTION norm_sf(communicator, norm_type, mesh, list_mode, v) RESULT(norm)
39  USE def_type_mesh
41  USE fem_tn_ns_mhd
42  IMPLICIT NONE
43  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
44  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
45  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
46  CHARACTER(*), INTENT(IN) :: norm_type
47  REAL(KIND=8) :: norm_loc, norm_tot, norm
48  INTEGER :: deb, fin, code
49  REAL(KIND=8) :: pi
50 #include "petsc/finclude/petsc.h"
51  mpi_comm, DIMENSION(2) :: communicator
52
53  norm = 0.d0
54  norm_tot = 0.d0
55  IF (mesh%me==0) THEN
56  norm_loc = 0.d0
57  ELSE
58
59  deb = start_of_string(norm_type)
60  fin = last_of_string(norm_type)
61  IF (norm_type(deb:fin)=='L2') THEN
62  norm_loc = norme_l2_champ(mesh, list_mode, v)
63  ELSE IF (norm_type(deb:fin)=='H1') THEN
64  norm_loc = sqrt(norme_l2_champ(mesh, list_mode, v)**2 &
65  + norme_h1_champ(mesh, list_mode, v)**2)
66  ELSE IF (norm_type(deb:fin)=='sH1') THEN
67  norm_loc = norme_h1_champ(mesh, list_mode, v)
68  ELSE IF (norm_type(deb:fin)=='div') THEN
69  norm_loc = norme_div(mesh, list_mode, v)
70  ELSE IF (norm_type(deb:fin)=='curl') THEN
71  norm_loc = norme_curl(mesh, list_mode, v)
72  ELSE
73  WRITE(*,*) ' BUG in norm, norm_type not programmed yet' , norm_type(deb:fin)
74  stop
75  END IF
76  END IF
77  CALL mpi_allreduce(norm_loc**2,norm_tot,1,mpi_double_precision, mpi_sum, communicator(2), code)
78  CALL mpi_allreduce(norm_tot,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
79  pi = acos(-1.d0)
80  norm = sqrt(norm*2*pi)
81
82  END FUNCTION norm_sf
83
84
85  FUNCTION norm_s(communicator, norm_type, mesh, list_mode, v) RESULT(norm)
86  USE def_type_mesh
88  USE fem_tn_ns_mhd
89  IMPLICIT NONE
90  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
91  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
92  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
93  CHARACTER(*), INTENT(IN) :: norm_type
94  REAL(KIND=8) :: norm_loc, norm_tot, norm
95  INTEGER :: deb, fin, code
96  REAL(KIND=8) :: pi
97 #include "petsc/finclude/petsc.h"
98  mpi_comm, DIMENSION(2) :: communicator
99
100  norm = 0.d0
101  norm_tot = 0.d0
102  IF (mesh%me==0) THEN
103  norm_loc = 0.d0
104  ELSE
105
106  deb = start_of_string(norm_type)
107  fin = last_of_string(norm_type)
108  IF (norm_type(deb:fin)=='L2') THEN
109  norm_loc = norme_l2_champ(mesh, list_mode, v)
110  ELSE IF (norm_type(deb:fin)=='H1') THEN
111  norm_loc = sqrt(norme_l2_champ(mesh, list_mode, v)**2 &
112  + norme_h1_champ(mesh, list_mode, v)**2)
113  ELSE IF (norm_type(deb:fin)=='sH1') THEN
114  norm_loc = norme_h1_champ(mesh, list_mode, v)
115  ELSE IF (norm_type(deb:fin)=='div') THEN
116  norm_loc = norme_div(mesh, list_mode, v)
117  ELSE IF (norm_type(deb:fin)=='curl') THEN
118  norm_loc = norme_curl(mesh, list_mode, v)
119  ELSE
120  WRITE(*,*) ' BUG in norm, norm_type not programmed yet' , norm_type(deb:fin)
121  stop
122  END IF
123  END IF
124  CALL mpi_allreduce(norm_loc**2,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
125  pi = acos(-1.d0)
126  norm = sqrt(norm*2*pi)
127
128  END FUNCTION norm_s
129
130  !DCQ, compute S_L1_norm using only the zero mode
131  FUNCTION norm_s_l1_zero_mode(communicator,mesh, list_mode, v) RESULT(norm)
132  USE def_type_mesh
133  USE chaine_caractere
134  USE fem_tn_ns_mhd
135  IMPLICIT NONE
136  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
137  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
138  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
139  REAL(KIND=8) :: norm_loc, norm_tot, norm
140  INTEGER :: code
141  INTEGER, DIMENSION(1) :: zero_mode
142  REAL(KIND=8) :: pi
143 #include "petsc/finclude/petsc.h"
144  mpi_comm, DIMENSION(2) :: communicator
145
146  norm = 0.d0
147  norm_tot = 0.d0
148  IF (mesh%me==0) THEN
149  norm_loc = 0.d0
150  ELSE
151  IF (minval(list_mode)==0) THEN !Just mode zero
152  zero_mode = minloc(list_mode)
153  norm_loc = norme_l1_one_mode(mesh, zero_mode(1), v)
154  ELSE
155  norm_loc=0;
156  END IF
157  END IF
158  CALL mpi_allreduce(norm_loc,norm,1,mpi_double_precision, mpi_sum, communicator(1), code)
159  pi = acos(-1.d0)
160  norm =(norm*2*pi)
161  END FUNCTION norm_s_l1_zero_mode
162
163  SUBROUTINE integration_mode(communicator,norm_loc,norm_tot)
164  IMPLICIT NONE
165  REAL(KIND=8), INTENT(IN) :: norm_loc
166  REAL(KIND=8), INTENT(OUT) :: norm_tot
167  INTEGER :: code
168  REAL(KIND=8) :: pi
169 #include "petsc/finclude/petsc.h"
170  mpi_comm :: communicator
171
172  pi = acos(-1.d0)
173  CALL mpi_allreduce(norm_loc,norm_tot,1,mpi_double_precision, mpi_sum, communicator, code)
174  !CN-AR Tue Jan 13 2009
175  norm_tot = norm_tot*(2*pi)
176  END SUBROUTINE integration_mode
177
178
179  FUNCTION norme_l2_champ_par(communicator, mesh, list_mode, v) RESULT(norm)
180  USE def_type_mesh
181  USE fem_tn_ns_mhd
182  IMPLICIT NONE
183  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
184  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
185  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
186  REAL(KIND=8) :: norm_loc, norm
187 #include "petsc/finclude/petsc.h"
188  mpi_comm :: communicator
189
190  IF (mesh%me==0) THEN
191  norm_loc = 0.d0
192  ELSE
193  norm_loc = norme_l2_champ(mesh, list_mode, v)
194  ENDIF
195  CALL integration_mode(communicator,norm_loc**2, norm)
196  norm = sqrt(norm)
197  END FUNCTION norme_l2_champ_par
198
199  FUNCTION norme_h1_champ_par(communicator, mesh, list_mode, v) RESULT(norm)
200  USE def_type_mesh
201  USE fem_tn_ns_mhd
202  IMPLICIT NONE
203  TYPE(mesh_type), INTENT(IN) :: mesh !type de maillage
204  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
205  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
206  REAL(KIND=8) :: norm_loc, norm
207 #include "petsc/finclude/petsc.h"
208  mpi_comm :: communicator
209
210  norm_loc = norme_h1_champ(mesh, list_mode, v)
211  CALL integration_mode(communicator,norm_loc**2, norm)
212  norm = sqrt(norm)
213  END FUNCTION norme_h1_champ_par
214
215  FUNCTION norme_div_par(communicator, H_mesh, list_mode, v) RESULT(norm)
216  USE def_type_mesh
217  USE fem_tn_ns_mhd
218  IMPLICIT NONE
219  TYPE(mesh_type), INTENT(IN) :: h_mesh !type de maillage
220  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
221  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
222  REAL(KIND=8) :: norm_loc, norm
223 #include "petsc/finclude/petsc.h"
224  mpi_comm :: communicator
225
226  norm_loc = norme_div(h_mesh, list_mode, v)
227  CALL integration_mode(communicator,norm_loc**2, norm)
228  norm = sqrt(norm)
229  END FUNCTION norme_div_par
230
231  FUNCTION norme_curl_par(communicator, H_mesh, list_mode, v) RESULT(norm)
232  USE def_type_mesh
233  USE fem_tn_ns_mhd
234  IMPLICIT NONE
235  TYPE(mesh_type), INTENT(IN) :: h_mesh !type de maillage
236  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
237  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
238  REAL(KIND=8) :: norm_loc, norm
239 #include "petsc/finclude/petsc.h"
240  mpi_comm :: communicator
241
242  norm_loc = norme_curl(h_mesh, list_mode, v)
243  CALL integration_mode(communicator,norm_loc**2, norm)
244  norm = sqrt(norm)
245  END FUNCTION norme_curl_par
246
247  SUBROUTINE angular_momentum(mesh, list_mode, field, moments_out)
248  USE def_type_mesh
249  IMPLICIT NONE
250  TYPE(mesh_type), INTENT(IN) :: mesh
251  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
252  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: field
253  REAL(KIND=8), DIMENSION(3), INTENT(OUT) :: moments_out
254
255  REAL(KIND=8), DIMENSION(3) :: moments_loc
256  REAL(KIND=8) :: m_x, m_y, ray, zed, urc, urs, utc, uts, uzc, uzs
257  INTEGER :: m, l, code, i
258 #include "petsc/finclude/petsc.h"
259
260  moments_loc=0.d0
261  moments_out=0.d0
262  DO i = 1, SIZE(list_mode)
263  IF (list_mode(i)==0) THEN
264  DO m = 1, mesh%me
265  DO l = 1, mesh%gauss%l_G
266  ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
267  utc = sum(field(mesh%jj(:,m),3,i)*mesh%gauss%ww(:,l))
268
269  moments_loc(3) = moments_loc(3) + ray**2*utc*mesh%gauss%rj(l,m)
270  END DO
271  END DO
272  ELSE IF (list_mode(i)==1) THEN
273  DO m = 1, mesh%me
274  DO l = 1, mesh%gauss%l_G
275  ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
276  zed = sum(mesh%rr(2,mesh%jj(:,m))*mesh%gauss%ww(:,l))
277
278  urc = sum(field(mesh%jj(:,m),1,i)*mesh%gauss%ww(:,l))
279  urs = sum(field(mesh%jj(:,m),2,i)*mesh%gauss%ww(:,l))
280  utc = sum(field(mesh%jj(:,m),3,i)*mesh%gauss%ww(:,l))
281  uts = sum(field(mesh%jj(:,m),4,i)*mesh%gauss%ww(:,l))
282  uzc = sum(field(mesh%jj(:,m),5,i)*mesh%gauss%ww(:,l))
283  uzs = sum(field(mesh%jj(:,m),6,i)*mesh%gauss%ww(:,l))
284
285  m_x = ray*(-zed*(urs+utc)+ray*uzs)
286  m_y = ray*(zed*(urc-uts)-ray*uzc)
287
288  moments_loc(1) = moments_loc(1) + m_x*mesh%gauss%rj(l,m)
289  moments_loc(2) = moments_loc(2) + m_y*mesh%gauss%rj(l,m)
290  END DO
291  END DO
292  END IF
293  END DO
294
295  CALL mpi_allreduce(moments_loc, moments_out, 3, mpi_double_precision, mpi_sum, petsc_comm_world, code)
296  moments_out(1) = acos(-1.d0)*moments_out(1)
297  moments_out(2) = acos(-1.d0)*moments_out(2)
298  moments_out(3) = 2*acos(-1.d0)*moments_out(3)
299
300  END SUBROUTINE angular_momentum
301
302  SUBROUTINE moments(communicator, H_mesh, list_mode, v, dipole_out, quadripole_out)
303  USE def_type_mesh
304  IMPLICIT NONE
305  TYPE(mesh_type), INTENT(IN) :: h_mesh
306  INTEGER, DIMENSION(:), INTENT(IN) :: list_mode
307  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: v
308  REAL(KIND=8), DIMENSION(3), INTENT(OUT):: dipole_out
310  REAL(KIND=8), DIMENSION(3) :: dipole
312  REAL(KIND=8) :: jr, ray, zed, pi
313  INTEGER :: m_max_c, mode, k, m, l, ni, i
314  INTEGER :: code
315  REAL(KIND=8), DIMENSION(6) :: c
316 #include "petsc/finclude/petsc.h"
317  mpi_comm :: communicator
318
319  pi = acos(-1.d0)
320  m_max_c = SIZE(list_mode)
321  dipole = 0
323
324  IF (SIZE(v,1)/=h_mesh%np .OR. SIZE(v,2)/=6 .OR. SIZE(v,3)/=m_max_c ) THEN
325  WRITE(*,*) ' BUG in MOMENTS', SIZE(v,1), h_mesh%np
326  stop
327  END IF
328  DO m = 1, h_mesh%me
329  DO l = 1, h_mesh%gauss%l_G
330
331  !--------On calcule le rayon et z du point gauss
332  ray = 0
333  zed = 0
334  DO ni = 1, h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
335  ray = ray + h_mesh%rr(1,i)*h_mesh%gauss%ww(ni,l)
336  zed = zed + h_mesh%rr(2,i)*h_mesh%gauss%ww(ni,l)
337  END DO
338  jr = ray * h_mesh%gauss%rj(l,m)
339
340  DO k=1, m_max_c
341  mode = list_mode(k)
342  IF (mode /=0 .AND. mode /=1 .AND. mode /=2) cycle
343
344  !--------Compute Curl------
345  c = 0
346  DO ni = 1,h_mesh%gauss%n_w; i = h_mesh%jj(ni,m)
347  !--------Composante r------
348  c(1) = c(1) + ( mode/ray*v(i,6,k)*h_mesh%gauss%ww(ni,l) &
349  - v(i,3,k)*h_mesh%gauss%dw(2,ni,l,m))
350  c(2) = c(2) + (-mode/ray*v(i,5,k)*h_mesh%gauss%ww(ni,l) &
351  - v(i,4,k)*h_mesh%gauss%dw(2,ni,l,m))
352  !--------Composante theta------
353  c(3) = c(3) + (v(i,1,k)*h_mesh%gauss%dw(2,ni,l,m) &
354  - v(i,5,k)*h_mesh%gauss%dw(1,ni,l,m))
355  c(4) = c(4) + (v(i,2,k)*h_mesh%gauss%dw(2,ni,l,m) &
356  - v(i,6,k)*h_mesh%gauss%dw(1,ni,l,m))
357  !--------Composante z------
358  c(5) = c(5) + (v(i,3,k)*h_mesh%gauss%dw(1,ni,l,m) &
359  + v(i,3,k)*h_mesh%gauss%ww(ni,l)/ray &
360  - mode/ray*v(i,2,k)*h_mesh%gauss%ww(ni,l))
361  c(6) = c(6) + (v(i,4,k)*h_mesh%gauss%dw(1,ni,l,m) &
362  + v(i,4,k)*h_mesh%gauss%ww(ni,l)/ray &
363  + mode/ray*v(i,1,k)*h_mesh%gauss%ww(ni,l))
364  ENDDO
365
367  IF (mode == 0) THEN
368  dipole(3) = dipole(3) + 2*pi*ray*c(3)*jr
374  ELSE IF (mode == 1) THEN
375  dipole(1) = dipole(1) + pi*(-zed*c(3) -zed*c(2) +ray*c(6))*jr
376  dipole(2) = dipole(2) + pi*(-zed*c(4) +zed*c(1) -ray*c(5))*jr
381  ELSE IF (mode == 2) THEN
386  END IF
387
388  END DO
389  END DO
390  END DO
391
392  !--------Collect from everybody------
393  CALL mpi_allreduce(dipole, dipole_out, 3,mpi_double_precision, &
394  mpi_sum, communicator, code)
396  mpi_sum, communicator, code)
397  RETURN
398  END SUBROUTINE moments
399
400  END MODULE tn_axi
real(kind=8) function norm_s_l1_zero_mode(communicator, mesh, list_mode, v)
Definition: tn_axi.f90:131
real(kind=8) function norme_h1_champ(mesh, list_mode, v)
subroutine angular_momentum(mesh, list_mode, field, moments_out)
Definition: tn_axi.f90:247
real(kind=8) function norme_div(H_mesh, list_mode, v)
Definition: tn_axi.f90:5
subroutine integration_mode(communicator, norm_loc, norm_tot)
Definition: tn_axi.f90:163
real(kind=8) function dot_product_champ(mesh, list_mode, v, w)
real(kind=8) function norme_curl_par(communicator, H_mesh, list_mode, v)
Definition: tn_axi.f90:231
real(kind=8) function norme_l2_champ(mesh, list_mode, v)
integer function, public start_of_string(string)
real(kind=8) function norme_h1_champ_par(communicator, mesh, list_mode, v)
Definition: tn_axi.f90:199
real(kind=8) function norme_curl(H_mesh, list_mode, v)
real(kind=8) function norm_s(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:85
integer function, public last_of_string(string)
real(kind=8) function norme_div_par(communicator, H_mesh, list_mode, v)
Definition: tn_axi.f90:215
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
Definition: tn_axi.f90:38
subroutine moments(communicator, H_mesh, list_mode, v, dipole_out, quadripole_out)
Definition: tn_axi.f90:302
real(kind=8) function norme_l1_one_mode(mesh, mode_idx, v)
real(kind=8) function norme_l2_champ_par(communicator, mesh, list_mode, v)
Definition: tn_axi.f90:179
real(kind=8) function dot_product_sf(communicator, mesh, list_mode, v, w)
Definition: tn_axi.f90:9