SFEMaNS  version 4.1 (work in progress) Reference documentation for SFEMaNS
fem_tn_axi.f90
Go to the documentation of this file.
1 !
2 !Authors: Jean-Luc Guermond, Raphael Laguerre, Luigi Quartapelle, Copyrights 1996, 2000, 2004
3 !
4 MODULE fem_tn_axi
5 CONTAINS
6
7
8 SUBROUTINE dot_product (mesh, ff, gg, t)
9 !===============================
10
11 ! sqrt(< f^2 >) ===> t
12
13  USE gauss_points
14
15  IMPLICIT NONE
16
17  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff, gg
18  REAL(KIND=8), INTENT(OUT) :: t
19
20  INTEGER :: m, l ,i ,ni
21
22  TYPE(mesh_type), TARGET :: mesh
23  INTEGER, DIMENSION(:,:), POINTER :: jj
24  INTEGER, POINTER :: me
25  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
26  REAL(KIND=8) :: ray
27
28  CALL gauss(mesh)
29  jj => mesh%jj
30  me => mesh%me
31  rr=> mesh%rr
32  t = 0
33
34  DO m = 1, me
35  DO l = 1, l_g
36
37 !===Compute radius of Gauss point
38  ray = 0
39  DO ni = 1, n_w; i = jj(ni,m)
40  ray = ray + rr(1,i)*ww(ni,l)
41  END DO
42
43  t = t + sum(ff(jj(:,m)) * ww(:,l))*sum(gg(jj(:,m)) * ww(:,l))*rj(l,m)*ray
44
45  ENDDO
46  ENDDO
47
48 END SUBROUTINE dot_product
49
50
51 SUBROUTINE ns_0 (mesh, ff, t)
52 !===============================
53
54 ! sqrt(< f^2 >) ===> t
55
56  USE gauss_points
57
58  IMPLICIT NONE
59
60  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
61  REAL(KIND=8), INTENT(OUT) :: t
62
63  INTEGER :: m, l ,i ,ni
64
65  TYPE(mesh_type), TARGET :: mesh
66  INTEGER, DIMENSION(:,:), POINTER :: jj
67  INTEGER, POINTER :: me
68  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
69  REAL(KIND=8) :: ray
70
71  CALL gauss(mesh)
72  jj => mesh%jj
73  me => mesh%me
74  rr=> mesh%rr
75  t = 0
76
77  DO m = 1, me
78  DO l = 1, l_g
79
80 !===Compute radius of Gauss point
81  ray = 0
82  DO ni = 1, n_w; i = jj(ni,m)
83  ray = ray + rr(1,i)*ww(ni,l)
84  END DO
85
86  t = t + sum(ff(jj(:,m)) * ww(:,l))**2 * rj(l,m)*ray
87
88  ENDDO
89  ENDDO
90
91  t = sqrt(t)
92
93 END SUBROUTINE ns_0
94
95 SUBROUTINE ns_l1 (mesh, ff, t)
96 !===============================
97
98 ! < f > ===> t
99
100  USE gauss_points
101
102  IMPLICIT NONE
103
104  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
105  REAL(KIND=8), INTENT(OUT) :: t
106
107  INTEGER :: m, l ,i ,ni
108
109  TYPE(mesh_type), TARGET :: mesh
110  INTEGER, DIMENSION(:,:), POINTER :: jj
111  INTEGER, POINTER :: me
112  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
113  REAL(KIND=8) :: ray
114
115  CALL gauss(mesh)
116  jj => mesh%jj
117  me => mesh%me
118  rr=> mesh%rr
119  t = 0
120
121  DO m = 1, me
122  DO l = 1, l_g
123
124 !===Compute radius of Gauss point
125  ray = 0
126  DO ni = 1, n_w; i = jj(ni,m)
127  ray = ray + rr(1,i)*ww(ni,l)
128  END DO
129
130  t = t + sum(ff(jj(:,m)) * ww(:,l))* rj(l,m)*ray
131
132  ENDDO
133  ENDDO
134
135 END SUBROUTINE ns_l1
136
137 SUBROUTINE average (mesh, ff, t)
138  !===============================
139
140  ! < f >/vol ===> t
141
142  USE gauss_points
143
144  IMPLICIT NONE
145
146  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
147  REAL(KIND=8), INTENT(OUT) :: t
148
149  INTEGER :: m, l ,i ,ni
150
151  TYPE(mesh_type), TARGET :: mesh
152  INTEGER, DIMENSION(:,:), POINTER :: jj
153  INTEGER, POINTER :: me
154  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
155  REAL(KIND=8) :: ray, vol
156
157
158  CALL gauss(mesh)
159  jj => mesh%jj
160  me => mesh%me
161  rr=> mesh%rr
162  t = 0
163  vol = 0.d0
164  IF (SIZE(ff)==mesh%np) THEN
165  DO m = 1, me
166  DO l = 1, l_g
167
168  !===Compute radius of Gauss point
169  ray = 0
170  DO ni = 1, n_w; i = jj(ni,m)
171  ray = ray + rr(1,i)*ww(ni,l)
172  END DO
173
174  t = t + sum(ff(jj(:,m)) * ww(:,l))* rj(l,m)*ray
175  vol = vol + rj(l,m)*ray
176  ENDDO
177  ENDDO
178  ELSE IF (SIZE(ff)==mesh%me) THEN
179  DO m = 1, me
180  DO l = 1, l_g
181
182  !===Compute radius of Gauss point
183  ray = 0
184  DO ni = 1, n_w; i = jj(ni,m)
185  ray = ray + rr(1,i)*ww(ni,l)
186  END DO
187
188  t = t + ff(m)* rj(l,m)*ray
189  vol = vol + rj(l,m)*ray
190  ENDDO
191  ENDDO
192  ELSE
193  WRITE(*,*) ' BUG in average '
194  stop
195  END IF
196  t = t / vol
197
198 END SUBROUTINE average
199
200 SUBROUTINE ns_anal_0 (mesh, ff, ff_anal, t)
201 !===============================
202
203 ! sqrt(< f^2 >) ===> t
204
205  USE gauss_points
206
207  IMPLICIT NONE
208
209  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
210  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff_anal
211  REAL(KIND=8), INTENT(OUT) :: t
212
213  INTEGER :: m, l ,i , ni, index
214
215  TYPE(mesh_type), TARGET :: mesh
216  INTEGER, DIMENSION(:,:), POINTER :: jj
217  INTEGER, POINTER :: me
218  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
219  REAL(KIND=8) :: ray, fl
220
221  CALL gauss(mesh)
222  jj => mesh%jj
223  me => mesh%me
224  rr=> mesh%rr
225  t = 0
226
227  DO m = 1, me
228  index = (m-1)*l_g
229  DO l = 1, l_g; index = index + 1
230
231 !===Compute radius of Gauss point
232  ray = 0
233  fl = 0
234  DO ni = 1, n_w; i = jj(ni,m)
235  ray = ray + rr(1,i)*ww(ni,l)
236  fl = fl + ff(i) * ww(ni,l)
237  END DO
238
239  t = t + (fl - ff_anal(index))**2 * ray* rj(l,m)
240
241  ENDDO
242  ENDDO
243
244  t = sqrt(t)
245
246 END SUBROUTINE ns_anal_0
247
248 SUBROUTINE ns_1 (mesh, ff, t)
249 !===============================
250
251 ! sqrt(<< (Df).(Df) >>) ===> t
252
253  USE gauss_points
254
255  IMPLICIT NONE
256
257  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
258  REAL(KIND=8), INTENT(OUT) :: t
259
260  INTEGER :: m, l, k,i,ni
261  REAL(KIND=8) :: s
262
263  TYPE(mesh_type), TARGET :: mesh
264  INTEGER, DIMENSION(:,:), POINTER :: jj
265  INTEGER, POINTER :: me
266  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
267  REAL(KIND=8) :: ray
268
269  CALL gauss(mesh)
270  jj => mesh%jj
271  me => mesh%me
272  rr=> mesh%rr
273
274  t = 0
275
276  DO m = 1, me
277  DO l = 1, l_g
278  s = 0
279  DO k = 1, k_d
280
281  !===Compute radius of Gauss point
282  ray = 0
283  DO ni = 1, n_w; i = jj(ni,m)
284  ray = ray + rr(1,i)*ww(ni,l)
285  END DO
286
287  s = s + sum(ff(jj(:,m)) * dw(k,:,l,m))**2*ray
288
289  ENDDO
290
291  t = t + s * rj(l,m)
292
293  ENDDO
294  ENDDO
295
296  t = sqrt(t)
297
298 END SUBROUTINE ns_1
299
300 SUBROUTINE nv_1(mesh, mod_max, ff, t)
301 !semi-norme vectoriel H1 pour une representation de Fourier
302 !(v_rc, v_rs, v_thc, v_ths, v_zc, v_zs)
303 !sqrt(<< (Dv).(Dv) >>) ===> t
304
305  USE gauss_points
306
307  IMPLICIT NONE
308
309  TYPE(mesh_type), TARGET :: mesh
310  INTEGER, INTENT(IN) :: mod_max
311  REAL(KIND=8), DIMENSION(6,mesh%np,0:mod_max), INTENT(IN) :: ff
312  REAL(KIND=8), INTENT(OUT) :: t
313
314  INTEGER :: m, l,i,ni ,j
315
316  INTEGER, DIMENSION(:,:), POINTER :: jj
317  INTEGER, POINTER :: me
318  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
319  REAL(KIND=8) :: ray
320  REAL(KIND=8), DIMENSION(2) :: div
321
322  CALL gauss(mesh)
323  jj => mesh%jj
324  me => mesh%me
325  rr=> mesh%rr
326
327  t = 0.d0
328  div = 0.d0
329
330  DO j = 0, mod_max
331  DO m = 1, me
332  DO l = 1, l_g
333 !===Compute radius of Gauss point
334  ray = 0
335  DO ni = 1, n_w; i = jj(ni,m)
336  ray = ray + rr(1,i)*ww(ni,l)
337  END DO
338
339  div(1) = (sum(ff(1,jj(:,m),j) * dw(1,:,l,m)) + &
340  sum(ff(1,jj(:,m),j) * ww(:,l))/ray + &
341  j/ray * sum(ff(4,jj(:,m),j) * ww(:,l)) + &
342  sum(ff(5,jj(:,m),j) * dw(2,:,l,m)))* ray*rj(l,m)
343
344  IF (j > 0) THEN
345  div(2) = (sum(ff(2,jj(:,m),j) * dw(1,:,l,m)) + &
346  sum(ff(2,jj(:,m),j) * ww(:,l))/ray - &
347  j/ray * sum(ff(3,jj(:,m),j) * ww(:,l)) + &
348  sum(ff(6,jj(:,m),j) * dw(2,:,l,m)))* ray*rj(l,m)
349  ENDIF
350
351  t = t +(div(1)**2+div(2)**2)
352
353  ENDDO
354  ENDDO
355
356  ENDDO
357
358  t = sqrt(t)
359
360
361 END SUBROUTINE nv_1
362
363 SUBROUTINE nv_0_cn(mesh, ff, p, t)
364  !semi-norme vectoriel H1 pour une representation de Fourier
365  !(v_rc, v_rs, v_zc, v_zs)
366
367  USE gauss_points
368
369  IMPLICIT NONE
370
371  TYPE(mesh_type), TARGET :: mesh
372  REAL(KIND=8), DIMENSION(mesh%np,6), INTENT(IN) :: ff
373  REAL(KIND=8), INTENT(OUT) :: t, p
374
375  INTEGER :: m, l, i, ni
376  REAL(KIND=8) :: s, rp, rt
377
378  INTEGER, DIMENSION(:,:), POINTER :: jj
379  INTEGER, POINTER :: me
380  REAL(KIND=8) ,DIMENSION(:,:), POINTER :: rr
381  REAL(KIND=8) :: ray
382
383  CALL gauss(mesh)
384  jj => mesh%jj
385  me => mesh%me
386  rr=> mesh%rr
387
388  rp = 0.d0
389  rt = 0.d0
390  s = 0.d0
391  DO m = 1, me
392  DO l = 1, l_g
393  !===Compute radius of Gauss point
394  ray = 0
395  DO ni = 1, n_w; i = jj(ni,m)
396  ray = ray + rr(1,i)*ww(ni,l)
397  END DO
398
399  rp = rp + sqrt(sum(ff(jj(:,m),1)* ww(:,l))**2 + sum(ff(jj(:,m),5)* ww(:,l))**2)*ray*rj(l,m)
400  rt = rt + abs(sum(ff(jj(:,m),3)* ww(:,l)))*ray*rj(l,m)
401  s = s + ray*rj(l,m)
402
403  ENDDO
404  ENDDO
405
406  p = rp/s
407  t = rt/s
408
409 END SUBROUTINE nv_0_cn
410
411 END MODULE fem_tn_axi
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 t
Definition: doc_intro.h:199
subroutine ns_1(mesh, ff, t)
Definition: fem_tn_axi.f90:248
subroutine nv_1(mesh, mod_max, ff, t)
Definition: fem_tn_axi.f90:300
subroutine ns_0(mesh, ff, t)
Definition: fem_tn_axi.f90:51
subroutine dot_product(mesh, ff, gg, t)
Definition: fem_tn_axi.f90:8
subroutine nv_0_cn(mesh, ff, p, t)
Definition: fem_tn_axi.f90:363
subroutine gauss(mesh)
subroutine ns_anal_0(mesh, ff, ff_anal, t)
Definition: fem_tn_axi.f90:200
subroutine ns_l1(mesh, ff, t)
Definition: fem_tn_axi.f90:95
subroutine average(mesh, visc)