SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
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)