SFEMaNS  version 5.3
Reference documentation for SFEMaNS
fem_s_direct_axi.f90
Go to the documentation of this file.
1 !
2 !
3 !
4 !Authors: Jean-Luc Guermond, Raphael Laguerre, Copyrights 2000, 2004
5 !
6 MODULE fem_s_axi
7 
8 CONTAINS
9 
10 
11  !----------------------------------------------------------------------
12  !----------------------------------------------------------------------
13  !Procedure pour le cas scalaire
14  !----------------------------------------------------------------------
15  !----------------------------------------------------------------------
16 
17  SUBROUTINE dirichlet (jjs, us, ff)
18  !=================================
19 
20  IMPLICIT NONE
21  INTEGER, DIMENSION(:), INTENT(IN) :: jjs
22  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: us
23  REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: ff
24 
25  INTEGER :: n, i
26 
27  DO n = 1, SIZE(jjs); i = jjs(n)
28  ff(i) = us(n)
29  END DO
30 
31  END SUBROUTINE dirichlet
32 
33 
34  SUBROUTINE moy(mesh,p,RESULT)
35  !===========================
36  !moyenne
37 
38  USE gauss_points
39 
40  IMPLICIT NONE
41  REAL(KIND=8), DIMENSION(:) , INTENT(IN) :: p
42  REAL(KIND=8) , INTENT(OUT) :: RESULT
43  REAL(KIND=8) :: vol
44  INTEGER :: m, l , i , ni,k
45 
46 
47  TYPE(mesh_type), TARGET :: mesh
48  INTEGER, DIMENSION(:,:), POINTER :: jj
49  INTEGER, POINTER :: me
50  REAL(KIND=8) :: ray
51 
52  CALL gauss(mesh)
53  jj => mesh%jj
54  me => mesh%me
55 
56 
57  result = 0.d0
58  vol = 0.d0
59 
60  DO m = 1, me
61  DO l = 1, l_g
62 
63  !===Compute radius of Gauss point
64  ray = 0
65  DO ni = 1, n_w; i = jj(ni,m)
66  ray = ray + mesh%rr(1,i)*ww(ni,l)
67  END DO
68 
69  result = result + sum(p(jj(:,m)) * ww(:,l))* ray* rj(l,m)
70  vol = vol + ray* rj(l,m)
71 
72  ENDDO
73  ENDDO
74 
75  result = result / vol
76 
77  END SUBROUTINE moy
78 
79 
80 
81  SUBROUTINE qs_00_ssr(mesh, ff, u0) !sans le rayon dans l'integration
82  !=================================
83 
84  !second membre simple
85  ! < w, f > ===> u0
86 
87  USE gauss_points
88 
89  IMPLICIT NONE
90  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
91  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
92 
93  INTEGER :: m, l , i , ni
94  REAL(KIND=8) :: fl
95 
96  TYPE(mesh_type), TARGET :: mesh
97  INTEGER, DIMENSION(:,:), POINTER :: jj
98  INTEGER, POINTER :: me
99  REAL(KIND=8) :: ray
100 
101  CALL gauss(mesh)
102  jj => mesh%jj
103  me => mesh%me
104 
105  u0 = 0
106 
107  DO m = 1, me
108  DO l = 1, l_g
109 
110  !===Compute radius of Gauss point
111  ray = 0
112  DO ni = 1, n_w; i = jj(ni,m)
113  ray = ray + mesh%rr(1,i)*ww(ni,l)
114  END DO
115 
116  fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
117 
118  DO ni = 1, n_w
119  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl
120  ENDDO
121 
122  ENDDO
123  ENDDO
124 
125  END SUBROUTINE qs_00_ssr
126 
127  SUBROUTINE qs_00(mesh, ff, u0) !avec le rayon dans l'integration
128  !=================================
129  !second membre simple
130  ! < w, f > ===> u0
131 
132  USE gauss_points
133 
134  IMPLICIT NONE
135  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
136  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
137 
138  INTEGER :: m, l , i , ni
139  REAL(KIND=8) :: fl
140 
141  TYPE(mesh_type), TARGET :: mesh
142  INTEGER, DIMENSION(:,:), POINTER :: jj
143  INTEGER, POINTER :: me
144  REAL(KIND=8) :: ray
145 
146  CALL gauss(mesh)
147  jj => mesh%jj
148  me => mesh%me
149 
150  u0 = 0
151 
152  DO m = 1, me
153  DO l = 1, l_g
154 
155  !===Compute radius of Gauss point
156  ray = 0
157  DO ni = 1, n_w; i = jj(ni,m)
158  ray = ray + mesh%rr(1,i)*ww(ni,l)
159  END DO
160 
161  fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
162 
163  DO ni = 1, n_w
164  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl *ray
165  ENDDO
166 
167  ENDDO
168  ENDDO
169 
170  END SUBROUTINE qs_00
171 
172 
173  SUBROUTINE qs_00_inst_3d_ssr(mesh,ff,T1,T2,dt,u0) !sans le rayon
174  !=================================
175  !second membre pour l'equation de la diffusion avec BDF2
176  ! < w, f > + <w,(4T(n)-T(n-1))/2dt> ===> u0 pour les cos(1) et sin(2)
177 
178  USE gauss_points
179 
180  IMPLICIT NONE
181  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
182  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
183  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T1 !T(noeud) ,Tn
184  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T2 !T(noeud) ,Tn-1
185  REAL(KIND=8), INTENT(IN) :: dt
186  INTEGER :: m, l , i , ni
187  REAL(KIND=8) :: fl
188 
189 
190  TYPE(mesh_type), TARGET :: mesh
191  INTEGER, DIMENSION(:,:), POINTER :: jj
192  INTEGER, POINTER :: me
193  REAL(KIND=8) :: ray
194 
195  CALL gauss(mesh)
196  jj => mesh%jj
197  me => mesh%me
198 
199  u0 = 0
200 
201  DO m = 1, me
202  DO l = 1, l_g
203 
204  !===Compute radius of Gauss point
205  ray = 0
206  DO ni = 1, n_w; i = jj(ni,m)
207  ray = ray + mesh%rr(1,i)*ww(ni,l)
208  END DO
209 
210  fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
211  + sum(4 * t1(jj(:,m)) * ww(:,l) - t2(jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)*ray
212 
213  DO ni = 1, n_w
214  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl
215  ENDDO
216 
217  ENDDO
218  ENDDO
219 
220  END SUBROUTINE qs_00_inst_3d_ssr
221 
222  SUBROUTINE qs_00_inst_3d(mesh,ff,T1,T2,dt,u0) !avec le rayon
223  !=================================
224 
225  !second membre pour l'equation de la diffusion avec BDF2
226  ! < w, f > + <w,(4T(n)-T(n-1))/2dt> ===> u0 pour les cos(1) et sin(2)
227 
228  USE gauss_points
229 
230  IMPLICIT NONE
231  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
232  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
233  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T1 !T(noeud) ,Tn
234  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T2 !T(noeud) ,Tn-1
235  REAL(KIND=8), INTENT(IN) :: dt
236  INTEGER :: m, l , i , ni
237  REAL(KIND=8) :: fl
238 
239 
240  TYPE(mesh_type), TARGET :: mesh
241  INTEGER, DIMENSION(:,:), POINTER :: jj
242  INTEGER, POINTER :: me
243  REAL(KIND=8) :: ray
244 
245  CALL gauss(mesh)
246  jj => mesh%jj
247  me => mesh%me
248 
249  u0 = 0
250 
251  DO m = 1, me
252  DO l = 1, l_g
253 
254  !===Compute radius of Gauss point
255  ray = 0
256  DO ni = 1, n_w; i = jj(ni,m)
257  ray = ray + mesh%rr(1,i)*ww(ni,l)
258  END DO
259 
260  fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
261  + sum(4 * t1(jj(:,m)) * ww(:,l) - t2(jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
262 
263  DO ni = 1, n_w
264  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl *ray
265  ENDDO
266 
267  ENDDO
268  ENDDO
269 
270  END SUBROUTINE qs_00_inst_3d
271 
272 
273  SUBROUTINE qs_00_inst_init_3d(mesh,ff,T1,dt,u0) !avec le rayon
274  !=================================
275  !second membre pour l'equation de la diffusion avec Euler retarde
276  ! < w, f > +<w,T(n)/dt> ===> u0 pour les cos(1) et sin(2)
277 
278  USE gauss_points
279 
280  IMPLICIT NONE
281  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
282  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
283  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T1 !T(noeud) ,Tn
284  REAL(KIND=8), INTENT(IN) :: dt
285  INTEGER :: m, l , i , ni
286  REAL(KIND=8) :: fl
287 
288 
289  TYPE(mesh_type), TARGET :: mesh
290  INTEGER, DIMENSION(:,:), POINTER :: jj
291  INTEGER, POINTER :: me
292  REAL(KIND=8) :: ray
293 
294  CALL gauss(mesh)
295  jj => mesh%jj
296  me => mesh%me
297 
298  u0 = 0
299 
300  DO m = 1, me
301  DO l = 1, l_g
302 
303  !===Compute radius of Gauss point
304  ray = 0
305  DO ni = 1, n_w; i = jj(ni,m)
306  ray = ray + mesh%rr(1,i)*ww(ni,l)
307  END DO
308 
309  fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
310  +sum(t1(jj(:,m)) * ww(:,l)) * rj(l,m) / dt
311 
312 
313  DO ni = 1, n_w
314  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * ray * fl
315  ENDDO
316 
317  ENDDO
318  ENDDO
319 
320  END SUBROUTINE qs_00_inst_init_3d
321 
322  SUBROUTINE qs_00_inst_init_3d_ssr(mesh,ff,T1,dt,u0) !sans le rayon
323  !=================================
324  !second membre pour l'equation de la diffusion avec Euler retarde
325  ! < w, f > +<w,T(n)/dt> ===> u0 pour les cos(1) et sin(2)
326 
327  USE gauss_points
328 
329  IMPLICIT NONE
330  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
331  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
332  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T1 !T(noeud) ,Tn
333  REAL(KIND=8), INTENT(IN) :: dt
334  INTEGER :: m, l , i , ni
335  REAL(KIND=8) :: fl
336 
337 
338  TYPE(mesh_type), TARGET :: mesh
339  INTEGER, DIMENSION(:,:), POINTER :: jj
340  INTEGER, POINTER :: me
341  REAL(KIND=8) :: ray
342 
343  CALL gauss(mesh)
344  jj => mesh%jj
345  me => mesh%me
346 
347  u0 = 0
348 
349  DO m = 1, me
350  DO l = 1, l_g
351 
352  !===Compute radius of Gauss point
353  ray = 0.d0
354  DO ni = 1, n_w; i = jj(ni,m)
355  ray = ray + mesh%rr(1,i)*ww(ni,l)
356  END DO
357 
358  fl = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m) &
359  +sum(t1(jj(:,m)) * ww(:,l)) * rj(l,m) / dt *ray
360 
361 
362  DO ni = 1, n_w
363  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * fl
364  ENDDO
365 
366  ENDDO
367  ENDDO
368 
369  END SUBROUTINE qs_00_inst_init_3d_ssr
370 
371 
372  SUBROUTINE qs_00_adv_diff_3d(eps,mode,mesh,ff,T1,T2,gg,dt,u0)
373  !=================================
374  !second membre pour l'advection diffusion d'un scalaire avec un BDF2
375  ! < w, f > + <w,(4T(n)-T(n-1))/2dt> + <w,eps* mV_theta /r (2T(n)-T(n-1)>
376  ! ===> u0 pour les cos(1) et sin(2)
377  ! le +ou- est donne par eps
378  !eps = 1 si on traite les cos et -1 sinon
379 
380 
381  USE gauss_points
382 
383  IMPLICIT NONE
384  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
385  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
386  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T1 !T(noeud) ,Tn
387  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T2 !T(noeud) ,Tn-1
388  REAL(KIND=8), INTENT(IN) :: dt
389  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg !champ de vitesse V_theta
390  REAL(KIND=8), INTENT(IN) :: eps
391  INTEGER , INTENT(IN) :: mode
392  INTEGER :: m, l , i , ni
393  REAL(KIND=8) :: fs , ft , fexp
394 
395 
396  TYPE(mesh_type), TARGET :: mesh
397  INTEGER, DIMENSION(:,:), POINTER :: jj
398  INTEGER, POINTER :: me
399  REAL(KIND=8) :: ray
400 
401  CALL gauss(mesh)
402  jj => mesh%jj
403  me => mesh%me
404 
405  u0 = 0
406 
407  DO m = 1, me
408  DO l = 1, l_g
409 
410  !===Compute radius of Gauss point
411  ray = 0
412  DO ni = 1, n_w; i = jj(ni,m)
413  ray = ray + mesh%rr(1,i)*ww(ni,l)
414  END DO
415 
416  fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
417  ft = sum(4 * t1(jj(:,m)) * ww(:,l) - t2(jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
418  fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(jj(:,m)) * ww(:,l))*rj(l,m) / ray
419 
420 
421  DO ni = 1, n_w
422  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * ray *( fs + ft + fexp)
423  ENDDO
424 
425  ENDDO
426  ENDDO
427 
428  END SUBROUTINE qs_00_adv_diff_3d
429 
430  SUBROUTINE qs_00_adv_diff_3d_ssr(eps,mode,mesh,ff,T1,T2,gg,dt,u0) !sans le rayon dans l'integral
431  !=================================
432  !second membre pour l'advection diffusion d'un scalaire avec un BDF2
433  ! (< w, f > + <w,(4T(n)-T(n-1))/2dt> + <w,esp* mV_theta /r (2T(n)-T(n-1)> )
434  ! ===> u0 pour les cos(1) et sin(2)
435  ! le +ou- est donne par eps
436  !eps = 1 si on traite les cos et -1 sinon
437 
438 
439  USE gauss_points
440 
441  IMPLICIT NONE
442  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
443  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
444  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: T1 !T(type,noeud) ,Tn
445  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: T2 !T(type,noeud) ,Tn-1
446  REAL(KIND=8), INTENT(IN) :: dt
447  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg !champ de vitesse V_theta
448  REAL(KIND=8), INTENT(IN) :: eps
449  INTEGER , INTENT(IN) :: mode
450  INTEGER :: m, l , i , ni
451  REAL(KIND=8) :: fs , ft , fexp
452 
453 
454  TYPE(mesh_type), TARGET :: mesh
455  INTEGER, DIMENSION(:,:), POINTER :: jj
456  INTEGER, POINTER :: me
457  REAL(KIND=8) :: ray
458 
459  CALL gauss(mesh)
460  jj => mesh%jj
461  me => mesh%me
462 
463  u0 = 0
464 
465  DO m = 1, me
466  DO l = 1, l_g
467 
468  !===Compute radius of Gauss point
469  ray = 0
470  DO ni = 1, n_w; i = jj(ni,m)
471  ray = ray + mesh%rr(1,i)*ww(ni,l)
472  END DO
473 
474  fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
475 
476  IF (eps.EQ.(-1.d0)) THEN !cas des cosinus
477  ft = sum(4 * t1(1,jj(:,m)) * ww(:,l) - t2(1,jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
478  fexp = eps * mode * sum(gg(3,jj(:,m)) * (2*t1(2,jj(:,m))-t2(2,jj(:,m))) &
479  * ww(:,l))*rj(l,m)/ray
480  ELSEIF (eps.EQ.(1.d0)) THEN !cas des sinus
481  ft = sum(4 * t1(2,jj(:,m)) * ww(:,l) - t2(2,jj(:,m))*ww(:,l))*rj(l,m) /(2* dt)
482  fexp = eps * mode * sum(gg(3,jj(:,m)) * (2*t1(1,jj(:,m))-t2(1,jj(:,m))) &
483  * ww(:,l))*rj(l,m)/ray
484  ELSE
485  WRITE(*,*) 'probleme ds calcul second membre avec epsilon'
486  ENDIF
487 
488  DO ni = 1, n_w
489  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * (fs + ray * (ft + fexp))
490  ENDDO
491 
492  ENDDO
493  ENDDO
494 
495  END SUBROUTINE qs_00_adv_diff_3d_ssr
496 
497 
498 
499  SUBROUTINE qs_00_adv_diff_init_3d(eps,mode,mesh,ff,T1,dt,gg,u0)
500  !=================================
501  !second membre pour l'advection diffusion d'un scalaire avec un Euler retarde
502  ! < w, f > +<w,T(n)/dt> + <w,eps* mV_theta /r T(n)>
503  ! ===> u0 pour les cos(1) et sin(2)
504 
505  USE gauss_points
506 
507  IMPLICIT NONE
508  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
509  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
510  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: T1 !T(noeud) ,Tn
511  REAL(KIND=8), INTENT(IN) :: dt
512  REAL(KIND=8), DIMENSION(:,:), INTENT(IN):: gg !champ de vitesse V_theta
513  REAL(KIND=8), INTENT(IN) :: eps
514  INTEGER , INTENT(IN) :: mode
515  INTEGER :: m, l , i , ni
516  REAL(KIND=8) :: fs , ft , fexp
517 
518 
519  TYPE(mesh_type), TARGET :: mesh
520  INTEGER, DIMENSION(:,:), POINTER :: jj
521  INTEGER, POINTER :: me
522  REAL(KIND=8) :: ray
523 
524  CALL gauss(mesh)
525  jj => mesh%jj
526  me => mesh%me
527 
528  u0 = 0
529 
530  DO m = 1, me
531  DO l = 1, l_g
532 
533  !===Compute radius of Gauss point
534  ray = 0
535  DO ni = 1, n_w; i = jj(ni,m)
536  ray = ray + mesh%rr(1,i)*ww(ni,l)
537  END DO
538 
539  fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
540  ft = sum(t1(jj(:,m)) * ww(:,l)) * rj(l,m) / dt
541  fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(jj(:,m)) * ww(:,l))*rj(l,m)/ray
542 
543  DO ni = 1, n_w
544  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * ray * (fs + ft + fexp)
545  ENDDO
546 
547  ENDDO
548  ENDDO
549 
550  END SUBROUTINE qs_00_adv_diff_init_3d
551 
552  SUBROUTINE qs_00_adv_diff_init_3d_ssr(eps,mode,mesh,ff,T1,dt,gg,u0) !sans le rayon ds l'integrale
553  !=================================
554  !second membre pour l'advection diffusion d'un scalaire avec un Euler retarde
555  ! < w, f > +<w,T(n)/dt> + <w,eps* mV_theta /r T(n)>
556  ! ===> u0 pour les cos(1) et sin(2)
557 
558  USE gauss_points
559 
560  IMPLICIT NONE
561  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff !analytique (noeud)
562  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: u0
563  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: T1 !T(noeud) ,Tn
564  REAL(KIND=8), INTENT(IN) :: dt
565  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg !champ de vitesse V_theta
566  REAL(KIND=8), INTENT(IN) :: eps
567  INTEGER , INTENT(IN) :: mode
568  INTEGER :: m, l , i , ni
569  REAL(KIND=8) :: fs , ft , fexp
570 
571 
572  TYPE(mesh_type), TARGET :: mesh
573  INTEGER, DIMENSION(:,:), POINTER :: jj
574  INTEGER, POINTER :: me
575  REAL(KIND=8) :: ray
576 
577  CALL gauss(mesh)
578  jj => mesh%jj
579  me => mesh%me
580 
581  u0 = 0
582 
583  DO m = 1, me
584  DO l = 1, l_g
585 
586  !===Compute radius of Gauss point
587  ray = 0
588  DO ni = 1, n_w; i = jj(ni,m)
589  ray = ray + mesh%rr(1,i)*ww(ni,l)
590  END DO
591 
592  IF (eps.EQ.(-1.d0)) THEN !cas des cosinus
593  fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
594  ft = sum(t1(1,jj(:,m)) * ww(:,l)) * rj(l,m) / dt
595  fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(2,jj(:,m)) * ww(:,l))*rj(l,m)/ray
596  ELSEIF (eps.EQ.(1.d0)) THEN !cas des sinus
597  fs = sum(ff(jj(:,m)) * ww(:,l)) * rj(l,m)
598  ft = sum(t1(2,jj(:,m)) * ww(:,l)) * rj(l,m) / dt
599  fexp = eps * mode * sum(gg(3,jj(:,m)) * t1(1,jj(:,m))* ww(:,l))*rj(l,m)/ray
600  ELSE
601  WRITE(*,*) 'problem calc second membre init avec epsilon'
602  ENDIF
603 
604  DO ni = 1, n_w
605  u0(jj(ni,m)) = u0(jj(ni,m)) + ww(ni,l) * (fs + ray *(ft + fexp))
606  ENDDO
607 
608  ENDDO
609  ENDDO
610 
611  END SUBROUTINE qs_00_adv_diff_init_3d_ssr
612 
613 
614  SUBROUTINE qs_00_lap(mesh,alpha,mode,ff,V2,V1,dt,u0) !avec le rayon
615  !=================================================
616  !calcul du second membre en entier pour le laplacien scalaire
617  !et un BDF2
618 
619 
620 
621  USE gauss_points
622 
623  IMPLICIT NONE
624  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
625  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
626  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !V(type,noeud) ,Vn
627  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !V(type,noeud) ,V(n-1)
628  REAL(KIND=8), INTENT(IN) :: dt , alpha
629  INTEGER , INTENT(IN) :: mode
630  INTEGER :: m, l , i , ni , j
631  REAL(KIND=8), DIMENSION(6) :: fs , ft
632 
633  !fs pour le terme de forcage
634  !ft pour le terme temporelle
635 
636  TYPE(mesh_type), TARGET :: mesh
637  INTEGER, DIMENSION(:,:), POINTER :: jj
638  INTEGER, POINTER :: me
639  REAL(KIND=8) :: ray
640 
641  CALL gauss(mesh)
642  jj => mesh%jj
643  me => mesh%me
644 
645  u0 = 0
646 
647  DO m = 1, me
648  DO l = 1, l_g
649 
650  !===Compute radius of Gauss point
651  ray = 0
652  DO ni = 1, n_w; i = jj(ni,m)
653  ray = ray + mesh%rr(1,i)*ww(ni,l)
654  END DO
655 
656  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
657 
658  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
659  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
660 
661  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
662 
663  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
664  ft(2) = 1/(2*dt)* sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
665 
666  !---------- calcul du second membre total
667 
668  DO j=5,6
669  DO ni = 1, n_w
670  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + ft(j))
671  ENDDO
672  ENDDO
673 
674 
675  ENDDO
676  ENDDO
677 
678  END SUBROUTINE qs_00_lap
679 
680  SUBROUTINE qs_00_lap_init(mesh,alpha,mode,ff,V,dt,u0) !avec le rayon
681  !=================================================
682  !calcul du second membre en entier pour le laplacien scalaire
683  !et un Euler retarde
684 
685 
686 
687  USE gauss_points
688 
689  IMPLICIT NONE
690  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
691  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
692  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V !V(type,noeud) ,Vn
693  REAL(KIND=8), INTENT(IN) :: dt , alpha
694  INTEGER , INTENT(IN) :: mode
695  INTEGER :: m, l , i , ni , j
696  REAL(KIND=8), DIMENSION(6) :: fs , ft
697 
698  !fs pour le terme de forcage
699  !ft pour le terme temporelle
700 
701  TYPE(mesh_type), TARGET :: mesh
702  INTEGER, DIMENSION(:,:), POINTER :: jj
703  INTEGER, POINTER :: me
704  REAL(KIND=8) :: ray
705 
706  CALL gauss(mesh)
707  jj => mesh%jj
708  me => mesh%me
709 
710  u0 = 0
711 
712  DO m = 1, me
713  DO l = 1, l_g
714 
715  !===Compute radius of Gauss point
716  ray = 0
717  DO ni = 1, n_w; i = jj(ni,m)
718  ray = ray + mesh%rr(1,i)*ww(ni,l)
719  END DO
720 
721  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
722 
723  fs(1) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
724  ft(1) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
725 
726  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
727 
728  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
729  ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
730 
731  !---------- calcul du second membre total
732 
733  DO j=5,6
734  DO ni = 1, n_w
735  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + ft(j))
736  ENDDO
737  ENDDO
738 
739 
740  ENDDO
741  ENDDO
742 
743  END SUBROUTINE qs_00_lap_init
744 
745  SUBROUTINE qs_00_rot(mesh, mode, V, Rot)
746 
747  USE gauss_points
748  IMPLICIT NONE
749  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: Rot !Rot(type,noeud)
750  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V !V(type,noeud) ,Vn
751  INTEGER , INTENT(IN) :: mode
752  INTEGER :: m, l , i , ni , j
753  REAL(KIND=8), DIMENSION(6) :: f
754  TYPE(mesh_type), TARGET :: mesh
755  INTEGER, DIMENSION(:,:), POINTER :: jj
756  INTEGER, POINTER :: me
757  REAL(KIND=8) :: ray, rayj
758  REAL(KIND=8), DIMENSION(6,mesh%gauss%n_w) :: V_loc
759  REAL(KIND=8), DIMENSION(2,mesh%gauss%n_w) :: dw_loc
760  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
761  CALL gauss(mesh)
762  jj => mesh%jj
763  me => mesh%me
764 
765  rot = 0.d0
766 
767  DO m = 1, me
768  j_loc = jj(:,m)
769  v_loc = v(:,j_loc)
770  DO l = 1, l_g
771  dw_loc = dw(:,:,l,m)
772 
773  !===Compute radius of Gauss point
774  ray = 0
775  DO ni = 1, n_w; i = jj(ni,m)
776  ray = ray + mesh%rr(1,i)*ww(ni,l)
777  END DO
778  rayj = ray*rj(l,m)
779  !----------------------
780 
781  f(1) = ( mode/ray*sum(v_loc(6,:)*ww(:,l))-sum(v_loc(3,:)*dw_loc(2,:)))*rayj
782  f(2) = (-mode/ray*sum(v_loc(5,:)*ww(:,l))-sum(v_loc(4,:)*dw_loc(2,:)))*rayj
783  f(3) = (sum(v_loc(1,:)*dw_loc(2,:))-sum(v_loc(5,:)*dw_loc(1,:)))*rayj
784  f(4) = (sum(v_loc(2,:)*dw_loc(2,:))-sum(v_loc(6,:)*dw_loc(1,:)))*rayj
785  f(5) = (1/ray*sum(v_loc(3,:)*ww(:,l))+sum(v_loc(3,:)*dw_loc(1,:))&
786  -mode/ray*sum(v_loc(2,:)*ww(:,l)))*rayj
787  f(6) = (1/ray*sum(v_loc(4,:)*ww(:,l))+sum(v_loc(4,:)*dw_loc(1,:))&
788  +mode/ray*sum(v_loc(1,:)*ww(:,l)))*rayj
789 
790  DO ni = 1, n_w
791  DO j= 1,6
792  rot(j,j_loc(ni)) = rot(j,j_loc(ni)) + ww(ni,l)*f(j)
793  ENDDO
794  ENDDO
795 
796  ENDDO
797  ENDDO
798 
799  END SUBROUTINE qs_00_rot
800 
801 
802 
803  SUBROUTINE qs_00_div(mesh,mode,V,u0)
804  !===============================================
805  !calucl de la divergence du champ V(6 composantes : cos,sin,cos,sin ...)
806  !en spectral pour theta
807 
808  USE gauss_points
809 
810  IMPLICIT NONE
811 
812 
813  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
814  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V !V(type,noeud) ,Vn
815  INTEGER , INTENT(IN) :: mode
816  INTEGER :: m, l , i , ni , j
817  REAL(KIND=8), DIMENSION(2) :: fr , fth , fz
818 
819  !fr pour le terme de forcage
820  !fth pour le terme temporelle
821  !fz
822 
823  TYPE(mesh_type), TARGET :: mesh
824  INTEGER, DIMENSION(:,:), POINTER :: jj
825  INTEGER, POINTER :: me
826  REAL(KIND=8) :: ray
827 
828  CALL gauss(mesh)
829  jj => mesh%jj
830  me => mesh%me
831 
832  u0 = 0
833  fr = 0.d0
834  fth = 0.d0
835  fz = 0.d0
836 
837  DO m = 1, me
838  DO l = 1, l_g
839 
840  !===Compute radius of Gauss point
841  ray = 0
842  DO ni = 1, n_w; i = jj(ni,m)
843  ray = ray + mesh%rr(1,i)*ww(ni,l)
844  END DO
845  !----------------------
846 
847  !fr(1) = -SUM(V(1,jj(:,m))*ww(:,l))* rj(l,m)
848  !fth(1) = mode*SUM(V(4,jj(:,m))*ww(:,l))* rj(l,m)
849  !fz(1) = -SUM(V(5,jj(:,m))*ww(:,l))* rj(l,m)
850 
851  !fr(2) = -SUM(V(2,jj(:,m))*ww(:,l))* rj(l,m)
852  !fth(2) = - mode*SUM(V(3,jj(:,m))*ww(:,l))* rj(l,m)
853  !fz(2) = -SUM(V(6,jj(:,m))*ww(:,l))* rj(l,m)
854 
855  fr(1) = -sum(v(1,jj(:,m))*dw(1,:,l,m))* rj(l,m)
856  fth(1) = mode*sum(v(4,jj(:,m))*ww(:,l))* rj(l,m)/ray
857  fz(1) = -sum(v(5,jj(:,m))*dw(2,:,l,m))* rj(l,m)
858 
859  fr(2) = -sum(v(2,jj(:,m))*dw(1,:,l,m))* rj(l,m)
860  fth(2) = - mode*sum(v(3,jj(:,m))*ww(:,l))* rj(l,m)/ray
861  fz(2) = -sum(v(6,jj(:,m))*dw(2,:,l,m))* rj(l,m)
862 
863 
864  !---------- calcul du second membre total
865 
866  DO j=1,2
867  DO ni = 1, n_w
868  !u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ray*(ww(ni,l)*fth(j) + dw(1,ni,l,m)*fr(j) &
869  ! +dw(2,ni,l,m)*fz(j))
870  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ray*ww(ni,l)*(fth(j) + fr(j) + fz(j))
871  ENDDO
872  ENDDO
873 
874 
875  ENDDO
876  ENDDO
877 
878  END SUBROUTINE qs_00_div
879 
880 
881 
882  !------------------------------------------------------------------------------
883  SUBROUTINE qs_01_div_hybrid_old (uu_mesh, pp_mesh,mode, gg, u0_c)
884  !================================================
885 
886  ! < w_c, D.g > ===> u0_c
887 
888  USE gauss_points
889 
890  IMPLICIT NONE
891 
892  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
893  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0_c
894  INTEGER , INTENT(IN) :: mode
895  TYPE(mesh_type), TARGET :: uu_mesh, pp_mesh
896  INTEGER, DIMENSION(:,:), POINTER :: jj, jj_c
897  INTEGER, POINTER :: me, nwc
898 
899  REAL(KIND=8), DIMENSION(pp_mesh%gauss%n_w, uu_mesh%gauss%l_G) :: w_c
900  INTEGER :: m, l, n, k,i,j
901  REAL(KIND=8) :: dgl
902  REAL(KIND=8), DIMENSION(6,uu_mesh%gauss%n_w) :: ggkn
903  REAL(KIND=8), DIMENSION(uu_mesh%gauss%k_d,uu_mesh%gauss%n_w) :: dwkn
904  REAL(KIND=8), DIMENSION(2,uu_mesh%gauss%n_w) :: u0_cn
905  REAL(KIND=8), DIMENSION(2) :: fr , fth , fz
906  REAL(KIND=8), DIMENSION(3,2) :: f
907  REAL(KIND=8) :: ray
908  INTEGER, DIMENSION(uu_mesh%gauss%n_w) :: j_loc
909  REAL(KIND=8), DIMENSION(2) :: smb
910  REAL(KIND=8):: user_time
911  !EXTERNAL :: user_time
912  REAL(KIND=8) :: tps, dummy, tt,tps1
913 
914 
915  CALL gauss(uu_mesh)
916  jj => uu_mesh%jj
917  jj_c => pp_mesh%jj
918  me => uu_mesh%me
919  nwc => pp_mesh%gauss%n_w
920 
921  DO l = 1, l_g
922  w_c(1,l) = ww(1,l) + 0.5*(ww(n_w-1,l) + ww(n_w,l))
923  w_c(2,l) = ww(2,l) + 0.5*(ww(n_w,l) + ww(4,l))
924  w_c(3,l) = ww(3,l) + 0.5*(ww(4,l) + ww(n_w-1,l))
925  END DO
926 
927  !write(*,*) 'me=',me
928  !write(*,*) 'lg=',l_G
929  !write(*,*) 'n_w=',n_w
930 
931  u0_c = 0.d0
932  tps = 0.d0
933 
934  DO m = 1, me
935  j_loc(:) = jj(:,m)
936  !u0_cn = 0.d0
937  DO l = 1, l_g
938 
939  !===Compute radius of P2 Gauss point
940  ray = 0
941  DO n = 1, n_w; i = jj(n,m)
942  ray = ray + uu_mesh%rr(1,i)*ww(n,l)
943  END DO
944  !----------------------
945 
946  !tt = user_time(dummy)
947  !calcul de la divergence sur les cos
948  f(1,1) = (ray*sum(gg(1,j_loc(:))*dw(1,:,l,m)) &
949  + sum(gg(1,j_loc(:))*ww(:,l)))
950  f(2,1) = mode*sum(gg(4,j_loc(:))*ww(:,l))
951  f(3,1) = sum(gg(5,j_loc(:))*dw(2,:,l,m)) * ray
952 
953  !calcul de la divergence sur les sin
954  f(1,2) = (ray*sum(gg(2,j_loc(:))*dw(1,:,l,m)) &
955  + sum(gg(2,j_loc(:))*ww(:,l)))
956  f(2,2) = - mode*sum(gg(3,j_loc(:))*ww(:,l))
957  f(3,2) = sum(gg(6,j_loc(:))*dw(2,:,l,m)) * ray
958 
959  f = f *rj(l,m)
960 
961 
962 
963  !tps = tps +user_time(dummy) -tt
964 
965  DO j=1,2
966  DO n=1, nwc
967  u0_c(j,jj_c(n,m)) = u0_c(j,jj_c(n,m)) + w_c(n,l)*(f(1,j)+f(2,j)+f(3,j))
968  ENDDO
969  ENDDO
970 
971  !calcul de la divergence sur les cos
972  !fr(1) = (ray*SUM(gg(1,jj(:,m))*dw(1,:,l,m)) &
973  ! + SUM(gg(1,jj(:,m))*ww(:,l)))
974  !fth(1) = mode*SUM(gg(4,jj(:,m))*ww(:,l))
975  !fz(1) = SUM(gg(5,jj(:,m))*dw(2,:,l,m)) * ray
976  !
977  ! !calcul de la divergence sur les sin
978  ! fr(2) = (ray*SUM(gg(2,jj(:,m))*dw(1,:,l,m)) &
979  ! + SUM(gg(2,jj(:,m))*ww(:,l)))
980  ! fth(2) = - mode*SUM(gg(3,jj(:,m))*ww(:,l))
981  ! fz(2) = SUM(gg(6,jj(:,m))*dw(2,:,l,m)) * ray
982 
983  ! fr(:) = fr(:) * rj(l,m)
984  ! fth(:) = fth(:) * rj(l,m)
985  ! fz(:) = fz(:) * rj(l,m)
986  !integration sur les P1
987  !j=1 pour les cos et j=2 pour les sin
988 
989  !DO n = 1, nwc
990  ! DO j = 1, 2
991  ! u0_cn(j,n) = u0_cn(j,n) + w_c(n,l)*(fr(j) + fz(j) + fth(j))
992  ! ENDDO
993  !END DO
994 
995  ENDDO
996 
997  !DO n = 1, nwc
998  ! DO j=1, 2
999  ! u0_c(j,jj_c(n,m)) = u0_c(j,jj_c(n,m)) + u0_cn(j,n)
1000  ! ENDDO
1001  !END DO
1002 
1003 
1004  ENDDO
1005  !write(*,*) 'tps=',tps
1006 
1007  END SUBROUTINE qs_01_div_hybrid_old
1008 
1009 
1010  !------------------------------------------------------------------------------
1011  SUBROUTINE qs_01_div_hybrid (uu_mesh, pp_mesh,mode, gg, u0_c)
1012  !================================================
1013 
1014  ! < w_c, D.g > ===> u0_c
1015 
1016  USE gauss_points
1017 
1018  IMPLICIT NONE
1019 
1020  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
1021  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0_c
1022  INTEGER , INTENT(IN) :: mode
1023  TYPE(mesh_type), TARGET :: uu_mesh, pp_mesh
1024  INTEGER, DIMENSION(:,:), POINTER :: jj, jj_c
1025  INTEGER, POINTER :: me, nwc
1026 
1027  REAL(KIND=8), DIMENSION(pp_mesh%gauss%n_w, uu_mesh%gauss%l_G) :: w_c
1028  INTEGER :: m, l, n, k,i,j
1029  REAL(KIND=8) :: dgl
1030  REAL(KIND=8), DIMENSION(6,uu_mesh%gauss%n_w) :: ggkn
1031  REAL(KIND=8), DIMENSION(uu_mesh%gauss%k_d,uu_mesh%gauss%n_w) :: dwkn
1032  REAL(KIND=8), DIMENSION(2,uu_mesh%gauss%n_w) :: u0_cn
1033  REAL(KIND=8), DIMENSION(2) :: fr , fth , fz
1034  REAL(KIND=8), DIMENSION(3,2) :: f
1035  REAL(KIND=8) :: ray
1036  INTEGER, DIMENSION(uu_mesh%gauss%n_w) :: j_loc
1037  REAL(KIND=8), DIMENSION(2) :: smb
1038  REAL(KIND=8):: user_time
1039  !EXTERNAL :: user_time
1040  REAL(KIND=8) :: tps, dummy, tt,tps1
1041 
1042 
1043  CALL gauss(uu_mesh)
1044  jj => uu_mesh%jj
1045  jj_c => pp_mesh%jj
1046  me => uu_mesh%me
1047  nwc => pp_mesh%gauss%n_w
1048 
1049  DO l = 1, l_g
1050  w_c(1,l) = ww(1,l) + 0.5*(ww(n_w-1,l) + ww(n_w,l))
1051  w_c(2,l) = ww(2,l) + 0.5*(ww(n_w,l) + ww(4,l))
1052  w_c(3,l) = ww(3,l) + 0.5*(ww(4,l) + ww(n_w-1,l))
1053  END DO
1054 
1055  u0_c = 0.d0
1056  tps = 0.d0
1057 
1058  DO m = 1, me
1059  j_loc(:) = jj(:,m)
1060  !u0_cn = 0.d0
1061  DO l = 1, l_g
1062 
1063  !===Compute radius of P2 Gauss point
1064  ray = 0
1065  DO n = 1, n_w; i = jj(n,m)
1066  ray = ray + uu_mesh%rr(1,i)*ww(n,l)
1067  END DO
1068  !----------------------
1069  tt = user_time(dummy)
1070 
1071  !calcul de la divergence sur les cos
1072  f(1,1) = (ray*sum(gg(1,j_loc(:))*dw(1,:,l,m)) &
1073  + sum(gg(1,j_loc(:))*ww(:,l)))
1074  f(2,1) = mode*sum(gg(4,j_loc(:))*ww(:,l))
1075  f(3,1) = sum(gg(5,j_loc(:))*dw(2,:,l,m)) * ray
1076 
1077  !calcul de la divergence sur les sin
1078  f(1,2) = (ray*sum(gg(2,j_loc(:))*dw(1,:,l,m)) &
1079  + sum(gg(2,j_loc(:))*ww(:,l)))
1080  f(2,2) = - mode*sum(gg(3,j_loc(:))*ww(:,l))
1081  f(3,2) = sum(gg(6,j_loc(:))*dw(2,:,l,m)) * ray
1082 
1083  f = f *rj(l,m)
1084 
1085  tps = tps +user_time(dummy) -tt
1086 
1087  DO j=1,2
1088  DO n=1, nwc
1089  u0_c(j,jj_c(n,m)) = u0_c(j,jj_c(n,m)) + w_c(n,l)*(f(1,j)+f(2,j)+f(3,j))
1090  ENDDO
1091  ENDDO
1092 
1093  ENDDO
1094 
1095  ENDDO
1096 
1097 
1098  END SUBROUTINE qs_01_div_hybrid
1099 
1100  SUBROUTINE qs_01_div_hybrid_2006 (uu_mesh, pp_mesh,mode, gg, u0_c)
1101  !================================================
1102 
1103  ! < w_c, D.g > ===> u0_c
1104 
1105  USE gauss_points
1106 
1107  IMPLICIT NONE
1108 
1109  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
1110  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0_c
1111  INTEGER , INTENT(IN) :: mode
1112  TYPE(mesh_type), TARGET :: uu_mesh, pp_mesh
1113  INTEGER, DIMENSION(:,:), POINTER :: jj, jj_c
1114  INTEGER, POINTER :: me, nwc
1115  REAL(KIND=8) :: x
1116  REAL(KIND=8), DIMENSION(pp_mesh%gauss%n_w, uu_mesh%gauss%l_G) :: w_c
1117  INTEGER :: m, l, n, k,i,j
1118  REAL(KIND=8), DIMENSION(3,2) :: f
1119  REAL(KIND=8) :: ray
1120  INTEGER, DIMENSION(uu_mesh%gauss%n_w) :: j_loc
1121 
1122 
1123  CALL gauss(uu_mesh)
1124  jj => uu_mesh%jj
1125  jj_c => pp_mesh%jj
1126  me => uu_mesh%me
1127  nwc => pp_mesh%gauss%n_w
1128 
1129 
1130 
1131  DO l = 1, l_g
1132  !cas P1 P2
1133  w_c(1,l) = ww(1,l) + 0.5*(ww(n_w-1,l) + ww(n_w,l))
1134  w_c(2,l) = ww(2,l) + 0.5*(ww(n_w,l) + ww(4,l))
1135  w_c(3,l) = ww(3,l) + 0.5*(ww(4,l) + ww(n_w-1,l))
1136  !cas P1 P1
1137  !w_c(1,l) = ww(1,l)
1138  !w_c(2,l) = ww(2,l)
1139  !w_c(3,l) = ww(3,l)
1140  END DO
1141 
1142  u0_c = 0.d0
1143 
1144  DO m = 1, me
1145  j_loc(:) = jj(:,m)
1146  DO l = 1, l_g
1147 
1148  !===Compute radius of P2 Gauss point
1149  ray = 0
1150  DO n = 1, n_w; i = jj(n,m)
1151  ray = ray + uu_mesh%rr(1,i)*ww(n,l)
1152  END DO
1153  !----------------------
1154 
1155  !calcul de la divergence sur les cos
1156  f(1,1) = (ray*sum(gg(j_loc,1)*dw(1,:,l,m)) &
1157  + sum(gg(j_loc,1)*ww(:,l)))
1158  f(2,1) = mode*sum(gg(j_loc,4)*ww(:,l))
1159  f(3,1) = sum(gg(j_loc,5)*dw(2,:,l,m)) * ray
1160 
1161  !calcul de la divergence sur les sin
1162  f(1,2) = (ray*sum(gg(j_loc,2)*dw(1,:,l,m)) &
1163  + sum(gg(j_loc,2)*ww(:,l)))
1164  f(2,2) =-mode*sum(gg(j_loc,3)*ww(:,l))
1165  f(3,2) = sum(gg(j_loc,6)*dw(2,:,l,m)) * ray
1166 
1167  f = f *rj(l,m)
1168 
1169  DO j=1,2
1170  x = f(1,j)+f(2,j)+f(3,j)
1171  DO n=1, nwc
1172  u0_c(jj_c(n,m),j) = u0_c(jj_c(n,m),j) + w_c(n,l)*x
1173  ENDDO
1174  ENDDO
1175 
1176  ENDDO
1177 
1178  ENDDO
1179 
1180  END SUBROUTINE qs_01_div_hybrid_2006
1181 
1182  !-------------------------------------
1183  SUBROUTINE qs_01_grad(mesh,mode,pp,u0)
1184  !=====================================
1185  !calcul le second membre pour le calcul du gradient pour un mode
1186  USE gauss_points
1187 
1188  IMPLICIT NONE
1189 
1190  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: pp
1191  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0
1192  INTEGER , INTENT(IN) :: mode
1193  INTEGER :: m, l , i , ni , j
1194  REAL(KIND=8),DIMENSION(6) :: fp
1195 
1196 
1197  TYPE(mesh_type), TARGET :: mesh
1198  INTEGER, DIMENSION(:,:), POINTER :: jj
1199  INTEGER, POINTER :: me
1200  REAL(KIND=8) :: ray
1201 
1202  CALL gauss(mesh)
1203  jj => mesh%jj
1204  me => mesh%me
1205 
1206  fp = 0.d0
1207  u0 = 0.d0
1208 
1209  DO m = 1, me
1210  DO l = 1, l_g
1211  !===Compute radius of Gauss point
1212  ray = 0
1213  DO ni = 1, n_w; i = jj(ni,m)
1214  ray = ray + mesh%rr(1,i)*ww(ni,l)
1215  END DO
1216 
1217 
1218  fp(1) = sum(pp(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
1219  fp(2) = sum(pp(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
1220  !JLG CN, 27 May 2009
1221  !Il y avait une inversion pp(1,...) et pp(2,...)
1222  fp(3) = mode/ray * sum(pp(2,jj(:,m))*ww(:,l)) * rj(l,m)
1223  fp(4) = -mode/ray * sum(pp(1,jj(:,m))*ww(:,l)) * rj(l,m)
1224  !JLG CN, 27 May 2009
1225  fp(5) = sum(pp(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
1226  fp(6) = sum(pp(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
1227 
1228  DO ni = 1, n_w
1229  DO j=1, 6
1230  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * fp(j)
1231  ENDDO
1232  ENDDO
1233 
1234  ENDDO
1235  ENDDO
1236 
1237  END SUBROUTINE qs_01_grad
1238 
1239  SUBROUTINE qs_01_grad_gl(mesh,mod_max,pp,u0)
1240  !=====================================
1241  !calcul le second membre pour le calcul du gradient pour m=0 -> mod_max
1242  USE gauss_points
1243 
1244  IMPLICIT NONE
1245 
1246  TYPE(mesh_type), TARGET :: mesh
1247  INTEGER , INTENT(IN) :: mod_max
1248  REAL(KIND=8), DIMENSION(2,mesh%np,0:mod_max), INTENT(IN) :: pp
1249  REAL(KIND=8), DIMENSION(6,mesh%np,0:mod_max), INTENT(OUT) :: u0
1250  INTEGER :: m, l , i , ni , j, k
1251  REAL(KIND=8),DIMENSION(6) :: fp
1252 
1253 
1254 
1255  INTEGER, DIMENSION(:,:), POINTER :: jj
1256  INTEGER, POINTER :: me
1257  REAL(KIND=8) :: ray
1258 
1259  CALL gauss(mesh)
1260  jj => mesh%jj
1261  me => mesh%me
1262 
1263  fp = 0.d0
1264  u0 = 0.d0
1265 
1266  DO k=0, mod_max
1267 
1268  DO m = 1, me
1269  DO l = 1, l_g
1270  !===Compute radius of Gauss point
1271  ray = 0
1272  DO ni = 1, n_w; i = jj(ni,m)
1273  ray = ray + mesh%rr(1,i)*ww(ni,l)
1274  END DO
1275 
1276  !DO k=0, mod_max
1277 
1278  !je multiplie par le rayon ici, car ray peut etre tres petit
1279  fp(1) = sum(pp(1,jj(:,m),k)*dw(1,:,l,m)) * rj(l,m)*ray
1280  fp(2) = sum(pp(2,jj(:,m),k)*dw(1,:,l,m)) * rj(l,m)*ray
1281 
1282  fp(3) = k* sum(pp(2,jj(:,m),k)*ww(:,l)) * rj(l,m) !/ray * ray
1283  fp(4) = -k* sum(pp(1,jj(:,m),k)*ww(:,l)) * rj(l,m) !/ray * ray
1284 
1285  fp(5) = sum(pp(1,jj(:,m),k)*dw(2,:,l,m)) * rj(l,m)*ray
1286  fp(6) = sum(pp(2,jj(:,m),k)*dw(2,:,l,m)) * rj(l,m)*ray
1287 
1288  DO ni = 1, n_w
1289  DO j=1, 6
1290  u0(j,jj(ni,m),k) = u0(j,jj(ni,m),k) + ww(ni,l) * fp(j)
1291  ENDDO
1292  ENDDO
1293 
1294  ENDDO
1295 
1296  ENDDO
1297 
1298  ENDDO
1299 
1300  END SUBROUTINE qs_01_grad_gl
1301  !---------------------------------------------------------------------
1302  !---------------------------------------------------------------------
1303  !procedure pour le cas vectoriel
1304  !---------------------------------------------------------------------
1305  !---------------------------------------------------------------------
1306 
1307 
1308  SUBROUTINE qs_00_inst_vect_3d(mesh,alpha,mode,ff,V1,V2,dt,u0) !avec le rayon
1309  !=================================
1310  !calcul le second membre en entier pour la diffusion d'un vecteur
1311  !avec un BDF2
1312  ! < w, f > +ou-alpha* <w , 2m/r^2*(2V(n)-V(n-1))> + (<w,(4V(n)-V(n-1))/2dt>
1313  ! ===> u0
1314  !cela depend des composantes et il y a en plus couplage entre les composantes
1315  !1 et 4, 2 et 3, 5 et 6
1316 
1317 
1318  USE gauss_points
1319 
1320  IMPLICIT NONE
1321  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
1322  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
1323  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !V(type,noeud) ,Vn-1
1324  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !V(type,noeud) ,Vn
1325  REAL(KIND=8), INTENT(IN) :: dt , alpha
1326  INTEGER , INTENT(IN) :: mode
1327  INTEGER :: m, l , i , ni , j
1328  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft
1329 
1330  !fs pour le terme de forcage
1331  !fexp pour le terme de vitesse explicite
1332  !ft pour le terme temporelle
1333 
1334 
1335  TYPE(mesh_type), TARGET :: mesh
1336  INTEGER, DIMENSION(:,:), POINTER :: jj
1337  INTEGER, POINTER :: me
1338  REAL(KIND=8) :: ray
1339 
1340  CALL gauss(mesh)
1341  jj => mesh%jj
1342  me => mesh%me
1343 
1344  u0 = 0
1345  fs = 0.d0
1346  fexp = 0.d0
1347  ft = 0.d0
1348 
1349  DO m = 1, me
1350  DO l = 1, l_g
1351 
1352  !===Compute radius of Gauss point
1353  ray = 0
1354  DO ni = 1, n_w; i = jj(ni,m)
1355  ray = ray + mesh%rr(1,i)*ww(ni,l)
1356  END DO
1357 
1358  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1359 
1360  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1361  fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
1362  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
1363 
1364  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1365 
1366  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1367  fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l)) * rj(l,m)
1368  ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
1369 
1370  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1371 
1372  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1373  fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
1374  ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
1375 
1376  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1377 
1378  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1379  fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
1380  ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
1381 
1382  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1383 
1384  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1385  fexp(5) = 0.d0
1386  ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
1387 
1388  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1389 
1390  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1391  fexp(6) = 0.d0
1392  ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
1393 
1394 
1395  !---------- calcul du second membre total
1396 
1397  DO j=1,6
1398  DO ni = 1, n_w
1399  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + fexp(j) + ft(j))
1400  ENDDO
1401  ENDDO
1402 
1403 
1404  ENDDO
1405  ENDDO
1406 
1407  END SUBROUTINE qs_00_inst_vect_3d
1408 
1409 
1410 
1411  SUBROUTINE qs_00_inst_vect_3d_ssr(mesh,alpha,mode,ff,V1,V2,dt,u0) !sans le rayon
1412  !pour le calculde l'integral du au forcage
1413  !=================================
1414  !calcul le second membre en entier pour la diffusion d'un vecteur
1415  !avec un BDF2
1416  ! < w, f > +ou-alpha* <w , 2m/r^2*(2V(n)-V(n-1))> + (<w,(4V(n)-V(n-1))/2dt>
1417  ! ===> u0
1418  !cela depend des composantes et il y a en plus couplage
1419 
1420 
1421  USE gauss_points
1422 
1423  IMPLICIT NONE
1424  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
1425  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
1426  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !T(type,noeud) ,Tn-1
1427  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !T(type,noeud) ,Tn
1428  REAL(KIND=8), INTENT(IN) :: dt , alpha
1429  INTEGER , INTENT(IN) :: mode
1430  INTEGER :: m, l , i , ni , j
1431  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft
1432 
1433  !fs pour le terme de forcage
1434  !fexp pour le terme de vitesse explicite
1435  !ft pour le terme temporelle
1436 
1437 
1438  TYPE(mesh_type), TARGET :: mesh
1439  INTEGER, DIMENSION(:,:), POINTER :: jj
1440  INTEGER, POINTER :: me
1441  REAL(KIND=8) :: ray
1442 
1443  CALL gauss(mesh)
1444  jj => mesh%jj
1445  me => mesh%me
1446 
1447  u0 = 0
1448 
1449  DO m = 1, me
1450  DO l = 1, l_g
1451 
1452  !===Compute radius of Gauss point
1453  ray = 0
1454  DO ni = 1, n_w; i = jj(ni,m)
1455  ray = ray + mesh%rr(1,i)*ww(ni,l)
1456  END DO
1457 
1458  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1459 
1460  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1461  fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
1462  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
1463 
1464  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1465 
1466  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1467  fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l)) * rj(l,m)
1468  ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
1469 
1470  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1471 
1472  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1473  fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
1474  ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
1475 
1476  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1477 
1478  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1479  fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
1480  ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
1481 
1482  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1483 
1484  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1485  fexp(5) = 0.d0
1486  ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
1487 
1488  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1489 
1490  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1491  fexp(6) = 0.d0
1492  ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
1493 
1494 
1495  !---------- calcul du second membre total
1496 
1497 
1498  DO j=1,6
1499  DO ni = 1, n_w
1500  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j)))
1501  ENDDO
1502  ENDDO
1503 
1504 
1505  ENDDO
1506  ENDDO
1507 
1508  END SUBROUTINE qs_00_inst_vect_3d_ssr
1509 
1510 
1511  SUBROUTINE qs_00_inst_vect_3d_init(mesh,alpha,mode,ff,V2,dt,u0) !avec le rayon
1512  !=================================
1513  !calcul le second membre en entier pour la diffusion d'un vecteur
1514  !avec un Euler retarde
1515  !calcul le second membre en entier avce les termes explicites
1516  ! < w, f > +ou-alpha* <w , 2m/r^2*T(n)> + (<w,V(n)/dt>
1517  ! ===> u0
1518  !cela depend des composantes et il y a en plus couplage
1519 
1520 
1521  USE gauss_points
1522 
1523  IMPLICIT NONE
1524  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
1525  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
1526  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !T(type,noeud) ,Tn
1527  REAL(KIND=8), INTENT(IN) :: dt , alpha
1528  INTEGER , INTENT(IN) :: mode
1529  INTEGER :: m, l , i , ni , j
1530  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft
1531 
1532  !fs pour le terme de forcage
1533  !fexp pour le terme de vitesse explicite
1534  !ft pour le terme temporelle
1535 
1536 
1537  TYPE(mesh_type), TARGET :: mesh
1538  INTEGER, DIMENSION(:,:), POINTER :: jj
1539  INTEGER, POINTER :: me
1540  REAL(KIND=8) :: ray
1541 
1542  CALL gauss(mesh)
1543  jj => mesh%jj
1544  me => mesh%me
1545 
1546  u0 = 0
1547 
1548  DO m = 1, me
1549  DO l = 1, l_g
1550 
1551  !===Compute radius of Gauss point
1552  ray = 0
1553  DO ni = 1, n_w; i = jj(ni,m)
1554  ray = ray + mesh%rr(1,i)*ww(ni,l)
1555  END DO
1556 
1557  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1558 
1559  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1560  fexp(1) = (-alpha*2*mode/ray**2)*sum(v2(4,jj(:,m))* ww(:,l)) * rj(l,m)
1561  ft(1) = 1/(dt) * sum(v2(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1562 
1563  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1564 
1565  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1566  fexp(2) = (alpha*2*mode/ray**2)*sum(v2(3,jj(:,m))* ww(:,l)) * rj(l,m)
1567  ft(2) = 1/(dt) * sum(v2(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1568 
1569  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1570 
1571  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1572  fexp(3) = (alpha*2*mode/ray**2)*sum(v2(2,jj(:,m))* ww(:,l)) * rj(l,m)
1573  ft(3) = 1/(dt) * sum(v2(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1574 
1575  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1576 
1577  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1578  fexp(4) = (-alpha*2*mode/ray**2)*sum(v2(1,jj(:,m))* ww(:,l)) * rj(l,m)
1579  ft(4) = 1/(dt) * sum(v2(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1580 
1581  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1582 
1583  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1584  fexp(5) = 0.d0
1585  ft(5) = 1/(dt) * sum(v2(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1586 
1587  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1588 
1589  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1590  fexp(6) = 0.d0
1591  ft(6) = 1/(dt) * sum(v2(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1592 
1593  !---------- calcul du second membre total
1594 
1595  DO j=1,6
1596  DO ni = 1, n_w
1597  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + fexp(j) + ft(j))
1598  ENDDO
1599  ENDDO
1600 
1601 
1602  ENDDO
1603  ENDDO
1604 
1605  END SUBROUTINE qs_00_inst_vect_3d_init
1606 
1607 
1608  SUBROUTINE qs_00_inst_vect_3d_init_ssr(mesh,alpha,mode,ff,V,dt,u0) !sans le rayon
1609  !dans l'intgration du forcage
1610  !=================================
1611  !calcul le second membre en entier pour la diffusion d'un vecteur
1612  !avec un Euler retarde
1613  !calcul le second membre en entier avce les termes explicites
1614  ! < w, f > +ou-alpha* <w , 2m/r^2*T(n)> + (<w,V(n)/dt>
1615  ! ===> u0
1616  !cela depend des composantes et il y a en plus couplage
1617 
1618 
1619  USE gauss_points
1620 
1621  IMPLICIT NONE
1622  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
1623  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
1624  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V !T(type,noeud) ,Tn
1625  REAL(KIND=8), INTENT(IN) :: dt , alpha
1626  INTEGER , INTENT(IN) :: mode
1627  INTEGER :: m, l , i , ni , j
1628  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft
1629 
1630  !fs pour le terme de forcage
1631  !fexp pour le terme de vitesse explicite
1632  !ft pour le terme temporelle
1633 
1634 
1635  TYPE(mesh_type), TARGET :: mesh
1636  INTEGER, DIMENSION(:,:), POINTER :: jj
1637  INTEGER, POINTER :: me
1638  REAL(KIND=8) :: ray
1639 
1640  CALL gauss(mesh)
1641  jj => mesh%jj
1642  me => mesh%me
1643 
1644  u0 = 0
1645  fs = 0.d0
1646  ft = 0.d0
1647  fexp = 0.d0
1648 
1649  DO m = 1, me
1650  DO l = 1, l_g
1651  !===Compute radius of Gauss point
1652  ray = 0
1653  DO ni = 1, n_w; i = jj(ni,m)
1654  ray = ray + mesh%rr(1,i)*ww(ni,l)
1655  END DO
1656 
1657  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1658 
1659  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1660  fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
1661  ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1662 
1663  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1664 
1665  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1666  fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
1667  ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1668 
1669  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1670 
1671  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1672  fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
1673  ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1674 
1675  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1676 
1677  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1678  fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
1679  ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1680 
1681  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1682 
1683  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1684  fexp(5) = 0.d0
1685  ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1686 
1687  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1688 
1689  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1690  fexp(6) = 0.d0
1691  ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1692 
1693  !---------- calcul du second membre total
1694 
1695 
1696 
1697  DO j=1,6
1698  DO ni = 1, n_w
1699  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j)))
1700  ENDDO
1701  ENDDO
1702 
1703 
1704  ENDDO
1705  ENDDO
1706 
1707  END SUBROUTINE qs_00_inst_vect_3d_init_ssr
1708 
1709 
1710  SUBROUTINE qs_00_adv_diff_vect_3d_init_ssr(mesh,alpha,mode,gg,ff,V,dt,u0) !sans le rayon
1711  !dans l'intgration du forcage
1712  !=================================
1713 
1714 
1715 
1716  USE gauss_points
1717 
1718  IMPLICIT NONE
1719  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
1720  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
1721  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V !T(type,noeud) ,Tn
1722  REAL(KIND=8), INTENT(IN) :: dt , alpha
1723  INTEGER , INTENT(IN) :: mode
1724  INTEGER :: m, l , i , ni , j
1725  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft , fv
1726  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg !champ de vitesse V_theta
1727  !fs pour le terme de forcage
1728  !fexp pour le terme de vitesse explicite
1729  !ft pour le terme temporelle
1730 
1731 
1732  TYPE(mesh_type), TARGET :: mesh
1733  INTEGER, DIMENSION(:,:), POINTER :: jj
1734  INTEGER, POINTER :: me
1735  REAL(KIND=8) :: ray
1736 
1737  CALL gauss(mesh)
1738  jj => mesh%jj
1739  me => mesh%me
1740 
1741  u0 = 0
1742  fs = 0.d0
1743  ft = 0.d0
1744  fexp = 0.d0
1745  fv = 0.d0
1746  DO m = 1, me
1747  DO l = 1, l_g
1748  !===Compute radius of Gauss point
1749  ray = 0
1750  DO ni = 1, n_w; i = jj(ni,m)
1751  ray = ray + mesh%rr(1,i)*ww(ni,l)
1752  END DO
1753 
1754  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1755 
1756  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1757  fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
1758  ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1759  fv(1) = -mode/ray*sum(gg(3,jj(:,m))*v(2,jj(:,m))* ww(:,l))* rj(l,m)
1760 
1761  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1762 
1763  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1764  fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
1765  ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1766  fv(2) = mode/ray*sum(gg(3,jj(:,m))*v(1,jj(:,m))* ww(:,l))* rj(l,m)
1767 
1768  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1769 
1770  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1771  fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
1772  ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1773  fv(3) = -mode/ray*sum(gg(3,jj(:,m))*v(4,jj(:,m))* ww(:,l))* rj(l,m)
1774 
1775  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1776 
1777  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1778  fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
1779  ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1780  fv(4) = mode/ray*sum(gg(3,jj(:,m))*v(3,jj(:,m))* ww(:,l))* rj(l,m)
1781 
1782  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1783 
1784  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1785  fexp(5) = 0.d0
1786  ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1787  fv(5) = -mode/ray*sum(gg(3,jj(:,m))*v(6,jj(:,m))* ww(:,l))* rj(l,m)
1788 
1789  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1790 
1791  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1792  fexp(6) = 0.d0
1793  ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1794  fv(6) = mode/ray*sum(gg(3,jj(:,m))*v(5,jj(:,m))* ww(:,l))* rj(l,m)
1795 
1796  !---------- calcul du second membre total
1797 
1798 
1799 
1800  DO j=1,6
1801  DO ni = 1, n_w
1802  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j) + fv(j)))
1803  ENDDO
1804  ENDDO
1805 
1806 
1807  ENDDO
1808  ENDDO
1809 
1810  END SUBROUTINE qs_00_adv_diff_vect_3d_init_ssr
1811 
1812 
1813  SUBROUTINE qs_00_adv_diff_vect_3d_init(mesh,alpha,mode,gg,ff,V,dt,u0) !avec le rayon
1814  !dans l'intgration du forcage
1815  !=================================
1816 
1817 
1818 
1819  USE gauss_points
1820 
1821  IMPLICIT NONE
1822  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
1823  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
1824  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V !T(type,noeud) ,Tn
1825  REAL(KIND=8), INTENT(IN) :: dt , alpha
1826  INTEGER , INTENT(IN) :: mode
1827  INTEGER :: m, l , i , ni , j
1828  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft , fv
1829  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg !champ de vitesse V_theta
1830  !fs pour le terme de forcage
1831  !fexp pour le terme de vitesse explicite
1832  !ft pour le terme temporelle
1833 
1834 
1835  TYPE(mesh_type), TARGET :: mesh
1836  INTEGER, DIMENSION(:,:), POINTER :: jj
1837  INTEGER, POINTER :: me
1838  REAL(KIND=8) :: ray
1839 
1840  CALL gauss(mesh)
1841  jj => mesh%jj
1842  me => mesh%me
1843 
1844  u0 = 0
1845  fs = 0.d0
1846  ft = 0.d0
1847  fexp = 0.d0
1848  fv = 0.d0
1849  DO m = 1, me
1850  DO l = 1, l_g
1851  !===Compute radius of Gauss point
1852  ray = 0
1853  DO ni = 1, n_w; i = jj(ni,m)
1854  ray = ray + mesh%rr(1,i)*ww(ni,l)
1855  END DO
1856 
1857  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1858 
1859  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1860  fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
1861  ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1862  fv(1) = -mode/ray*sum(gg(3,jj(:,m))*v(2,jj(:,m))* ww(:,l))* rj(l,m)
1863 
1864  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1865 
1866  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1867  fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
1868  ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1869  fv(2) = mode/ray*sum(gg(3,jj(:,m))*v(1,jj(:,m))* ww(:,l))* rj(l,m)
1870 
1871  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1872 
1873  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1874  fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
1875  ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1876  fv(3) = -mode/ray*sum(gg(3,jj(:,m))*v(4,jj(:,m))* ww(:,l))* rj(l,m)
1877 
1878  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1879 
1880  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1881  fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
1882  ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1883  fv(4) = mode/ray*sum(gg(3,jj(:,m))*v(3,jj(:,m))* ww(:,l))* rj(l,m)
1884 
1885  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1886 
1887  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1888  fexp(5) = 0.d0
1889  ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1890  fv(5) = -mode/ray*sum(gg(3,jj(:,m))*v(6,jj(:,m))* ww(:,l))* rj(l,m)
1891 
1892  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1893 
1894  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1895  fexp(6) = 0.d0
1896  ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
1897  fv(6) = mode/ray*sum(gg(3,jj(:,m))*v(5,jj(:,m))* ww(:,l))* rj(l,m)
1898 
1899  !---------- calcul du second membre total
1900 
1901 
1902 
1903  DO j=1,6
1904  DO ni = 1, n_w
1905  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray*(fs(j) + fexp(j) + ft(j) + fv(j))
1906  ENDDO
1907  ENDDO
1908 
1909 
1910  ENDDO
1911  ENDDO
1912 
1913  END SUBROUTINE qs_00_adv_diff_vect_3d_init
1914 
1915  SUBROUTINE qs_00_adv_diff_vect_3d_ssr(mesh,alpha,mode,gg,ff,V1,V2,dt,u0) !sans le rayon
1916  !pour le calculde l'integral du au forcage
1917  !=================================
1918 
1919 
1920  USE gauss_points
1921 
1922  IMPLICIT NONE
1923  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
1924  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
1925  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !T(type,noeud) ,Tn-1
1926  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !T(type,noeud) ,Tn
1927  REAL(KIND=8), INTENT(IN) :: dt , alpha
1928  INTEGER , INTENT(IN) :: mode
1929  INTEGER :: m, l , i , ni , j
1930  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft , fv
1931 
1932  !fs pour le terme de forcage
1933  !fexp pour le terme de vitesse explicite
1934  !ft pour le terme temporelle
1935 
1936 
1937  TYPE(mesh_type), TARGET :: mesh
1938  INTEGER, DIMENSION(:,:), POINTER :: jj
1939  INTEGER, POINTER :: me
1940  REAL(KIND=8) :: ray
1941  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg !champ de vitesse V_theta
1942  CALL gauss(mesh)
1943  jj => mesh%jj
1944  me => mesh%me
1945 
1946  u0 = 0
1947 
1948  fs = 0.d0
1949  fexp = 0.d0
1950  ft = 0.d0
1951  fv = 0.d0
1952 
1953  DO m = 1, me
1954  DO l = 1, l_g
1955 
1956  !===Compute radius of Gauss point
1957  ray = 0
1958  DO ni = 1, n_w; i = jj(ni,m)
1959  ray = ray + mesh%rr(1,i)*ww(ni,l)
1960  END DO
1961 
1962  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1963 
1964  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
1965  fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
1966  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
1967  fv(1) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(2,jj(:,m)) - v1(2,jj(:,m)))* ww(:,l))*rj(l,m)
1968 
1969  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1970 
1971  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
1972  fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
1973  ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) *rj(l,m)
1974  fv(2) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(1,jj(:,m)) - v1(1,jj(:,m)))* ww(:,l))*rj(l,m)
1975 
1976  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1977 
1978  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
1979  fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
1980  ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
1981  fv(3) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(4,jj(:,m)) - v1(4,jj(:,m)))* ww(:,l))*rj(l,m)
1982 
1983  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1984 
1985  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
1986  fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
1987  ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
1988  fv(4) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(3,jj(:,m)) - v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
1989 
1990  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1991 
1992  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
1993  fexp(5) = 0.d0
1994  ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
1995  fv(5) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(6,jj(:,m)) - v1(6,jj(:,m)))* ww(:,l))*rj(l,m)
1996 
1997  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1998 
1999  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2000  fexp(6) =0.d0
2001  ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2002  fv(6) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(5,jj(:,m)) - v1(5,jj(:,m)))* ww(:,l))*rj(l,m)
2003 
2004  !---------- calcul du second membre total
2005 
2006 
2007  DO j=1,6
2008  DO ni = 1, n_w
2009  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(fexp(j) + ft(j) + fv(j)))
2010  ENDDO
2011  ENDDO
2012 
2013 
2014  ENDDO
2015  ENDDO
2016 
2017  END SUBROUTINE qs_00_adv_diff_vect_3d_ssr
2018 
2019  SUBROUTINE qs_00_adv_diff_vect_3d(mesh,alpha,mode,gg,ff,V1,V2,dt,u0) !sans le rayon
2020  !pour le calculde l'integral du au forcage
2021  !=================================
2022 
2023 
2024  USE gauss_points
2025 
2026  IMPLICIT NONE
2027  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
2028  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
2029  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !T(type,noeud) ,Tn-1
2030  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !T(type,noeud) ,Tn
2031  REAL(KIND=8), INTENT(IN) :: dt , alpha
2032  INTEGER , INTENT(IN) :: mode
2033  INTEGER :: m, l , i , ni , j
2034  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft , fv
2035 
2036  !fs pour le terme de forcage
2037  !fexp pour le terme de vitesse explicite
2038  !ft pour le terme temporelle
2039 
2040 
2041  TYPE(mesh_type), TARGET :: mesh
2042  INTEGER, DIMENSION(:,:), POINTER :: jj
2043  INTEGER, POINTER :: me
2044  REAL(KIND=8) :: ray
2045  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg !champ de vitesse V_theta
2046  CALL gauss(mesh)
2047  jj => mesh%jj
2048  me => mesh%me
2049 
2050  u0 = 0
2051 
2052  fs = 0.d0
2053  fexp = 0.d0
2054  ft = 0.d0
2055  fv = 0.d0
2056 
2057  DO m = 1, me
2058  DO l = 1, l_g
2059 
2060  !===Compute radius of Gauss point
2061  ray = 0
2062  DO ni = 1, n_w; i = jj(ni,m)
2063  ray = ray + mesh%rr(1,i)*ww(ni,l)
2064  END DO
2065 
2066  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
2067 
2068  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2069  fexp(1) = -(alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
2070  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2071  fv(1) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(2,jj(:,m)) - v1(2,jj(:,m)))* ww(:,l))*rj(l,m)
2072 
2073  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
2074 
2075  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2076  fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
2077  ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) *rj(l,m)
2078  fv(2) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(1,jj(:,m)) - v1(1,jj(:,m)))* ww(:,l))*rj(l,m)
2079 
2080  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
2081 
2082  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2083  fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
2084  ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2085  fv(3) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(4,jj(:,m)) - v1(4,jj(:,m)))* ww(:,l))*rj(l,m)
2086 
2087  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
2088 
2089  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2090  fexp(4) = -(alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
2091  ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2092  fv(4) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(3,jj(:,m)) - v1(3,jj(:,m)))* ww(:,l))*rj(l,m)
2093 
2094  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
2095 
2096  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2097  fexp(5) = 0.d0
2098  ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2099  fv(5) = -mode/ray*sum(gg(3,jj(:,m))*(2*v2(6,jj(:,m)) - v1(6,jj(:,m)))* ww(:,l))*rj(l,m)
2100 
2101  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
2102 
2103  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2104  fexp(6) =0.d0
2105  ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2106  fv(6) = mode/ray*sum(gg(3,jj(:,m))*(2*v2(5,jj(:,m)) - v1(5,jj(:,m)))* ww(:,l))*rj(l,m)
2107 
2108  !---------- calcul du second membre total
2109 
2110 
2111  DO j=1,6
2112  DO ni = 1, n_w
2113  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray* (fs(j) + fexp(j) + ft(j) + fv(j))
2114  ENDDO
2115  ENDDO
2116 
2117 
2118  ENDDO
2119  ENDDO
2120 
2121  END SUBROUTINE qs_00_adv_diff_vect_3d
2122 
2123 
2124  SUBROUTINE qs_00_stokes_3d_init(mesh,alpha,mode,ff,P,V,dt,u0) !sans le rayon
2125  !dans l'intgration du forcage
2126  !=================================
2127  !calcul le second membre en entier pour le probleme de Stokes avec un Euler retarde
2128  ! < w, f > +ou-alpha* <w , 2m/r^2*T(n)> + (<w,(V(n))/dt> + <w,grad(P)>
2129  ! ===> u0
2130  !cela depend des composantes et il y a en plus couplage
2131 
2132 
2133  USE gauss_points
2134 
2135  IMPLICIT NONE
2136  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
2137  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
2138  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V !T(type,noeud) ,Tn
2139  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
2140  REAL(KIND=8), INTENT(IN) :: dt , alpha
2141  INTEGER , INTENT(IN) :: mode
2142  INTEGER :: m, l , i , ni , j
2143  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft , fp
2144 
2145  !fs pour le terme de forcage
2146  !fexp pour le terme de vitesse explicite
2147  !ft pour le terme temporelle
2148  !fp gradient de pression
2149 
2150  TYPE(mesh_type), TARGET :: mesh
2151  INTEGER, DIMENSION(:,:), POINTER :: jj
2152  INTEGER, POINTER :: me
2153  REAL(KIND=8) :: ray
2154 
2155  CALL gauss(mesh)
2156  jj => mesh%jj
2157  me => mesh%me
2158 
2159  u0 = 0
2160  fs = 0.d0
2161  ft = 0.d0
2162  fexp = 0.d0
2163  fp = 0.d0
2164 
2165  DO m = 1, me
2166  DO l = 1, l_g
2167  !===Compute radius of Gauss point
2168  ray = 0
2169  DO ni = 1, n_w; i = jj(ni,m)
2170  ray = ray + mesh%rr(1,i)*ww(ni,l)
2171  END DO
2172 
2173  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
2174 
2175  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2176  fexp(1) = -(alpha*2*mode/ray**2)*sum(v(4,jj(:,m))* ww(:,l)) * rj(l,m)
2177  ft(1) = 1/(dt) * sum(v(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2178  !fp(1) = SUM(P(1,jj(:,m))*ww(:,l)) * rj(l,m)
2179  fp(1) = sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2180  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
2181 
2182  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2183  fexp(2) = (alpha*2*mode/ray**2)*sum(v(3,jj(:,m))* ww(:,l)) * rj(l,m)
2184  ft(2) = 1/(dt) * sum(v(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2185  !fp(2) = SUM(P(2,jj(:,m))*ww(:,l)) * rj(l,m)
2186  fp(2) = sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2187  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
2188 
2189  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2190  fexp(3) = (alpha*2*mode/ray**2)*sum(v(2,jj(:,m))* ww(:,l)) * rj(l,m)
2191  ft(3) = 1/(dt) * sum(v(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2192  fp(3) = sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2193 
2194  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
2195 
2196  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2197  fexp(4) = -(alpha*2*mode/ray**2)*sum(v(1,jj(:,m))* ww(:,l)) * rj(l,m)
2198  ft(4) = 1/(dt) * sum(v(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2199  fp(4) = -sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2200 
2201  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
2202 
2203  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2204  fexp(5) = 0.d0
2205  ft(5) = 1/(dt) * sum(v(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2206  !fp(5) = SUM(P(1,jj(:,m))*ww(:,l)) * rj(l,m)
2207  fp(5) = sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2208  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
2209 
2210  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2211  fexp(6) = 0.d0
2212  ft(6) = 1/(dt) * sum(v(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2213  !fp(6) = SUM(P(2,jj(:,m))*ww(:,l)) * rj(l,m)
2214  fp(6) = sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2215  !---------- calcul du second membre total
2216 
2217 
2218 
2219  DO j=1,6
2220  DO ni = 1, n_w
2221  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray* (fs(j) + fexp(j) + ft(j) + fp(j))
2222 
2223  !IF (j == 3.or.j == 4) THEN
2224  ! u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray *fp(j)
2225  !ELSEIF (j == 1.or. j == 2) THEN
2226  ! u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * fp(j)
2227  !ELSEIF (j == 5.or. j == 6) THEN
2228  ! u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * fp(j)
2229  !ELSE
2230  !ELSEIF (j == 1.or. j == 2) THEN
2231  ! u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + dw(1,ni,l,m) * ray * fp(j)
2232  !ELSEIF (j == 5.or. j == 6) THEN
2233  ! u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + dw(2,ni,l,m) * ray * fp(j)
2234  !ELSE
2235  ! WRITE(*,*) 'probleme ds calcul second membre avec j,subrout stokes'
2236  !ENDIF
2237  ENDDO
2238  ENDDO
2239 
2240 
2241  ENDDO
2242  ENDDO
2243 
2244  END SUBROUTINE qs_00_stokes_3d_init
2245 
2246 
2247  SUBROUTINE qs_00_stokes_3d(mesh,alpha,mode,ff,V1,V2,P,dt,u0) !avec le rayon
2248  !=================================
2249 
2250 
2251 
2252  USE gauss_points
2253 
2254  IMPLICIT NONE
2255  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
2256  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
2257  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !V(type,noeud) ,Vn-1
2258  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !V(type,noeud) ,Vn
2259  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
2260  REAL(KIND=8), INTENT(IN) :: dt , alpha
2261  INTEGER , INTENT(IN) :: mode
2262  INTEGER :: m, l , i , ni , j
2263  REAL(KIND=8), DIMENSION(6) :: fs , fexp , ft , fp
2264 
2265  !fs pour le terme de forcage
2266  !fexp pour le terme de vitesse explicite
2267  !ft pour le terme temporelle
2268  !fp pour le gradient de pression
2269 
2270  TYPE(mesh_type), TARGET :: mesh
2271  INTEGER, DIMENSION(:,:), POINTER :: jj
2272  INTEGER, POINTER :: me
2273  REAL(KIND=8) :: ray
2274 
2275  CALL gauss(mesh)
2276  jj => mesh%jj
2277  me => mesh%me
2278 
2279  u0 = 0.d0
2280  fs = 0.d0
2281  fexp = 0.d0
2282  ft = 0.d0
2283  fp = 0.d0
2284  DO m = 1, me
2285  DO l = 1, l_g
2286 
2287  !===Compute radius of Gauss point
2288  ray = 0
2289  DO ni = 1, n_w; i = jj(ni,m)
2290  ray = ray + mesh%rr(1,i)*ww(ni,l)
2291  END DO
2292 
2293  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
2294 
2295  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2296  fexp(1) = (-alpha*2*mode/ray**2)*sum((2*v2(4,jj(:,m))-v1(4,jj(:,m)))* ww(:,l)) * rj(l,m)
2297  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2298  fp(1) = -sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2299 
2300  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
2301 
2302  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2303  fexp(2) = (alpha*2*mode/ray**2)*sum((2*v2(3,jj(:,m))-v1(3,jj(:,m)))* ww(:,l)) * rj(l,m)
2304  ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
2305  fp(2) = -sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2306 
2307  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
2308 
2309  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2310  fexp(3) = (alpha*2*mode/ray**2)*sum((2*v2(2,jj(:,m))-v1(2,jj(:,m)))* ww(:,l)) * rj(l,m)
2311  ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2312  fp(3) = -sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2313 
2314  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
2315 
2316  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2317  fexp(4) = (-alpha*2*mode/ray**2)*sum((2*v2(1,jj(:,m))-v1(1,jj(:,m)))* ww(:,l)) * rj(l,m)
2318  ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2319  fp(4) = sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2320 
2321  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
2322 
2323  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2324  fexp(5) = 0.d0
2325  ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2326  fp(5) = -sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2327 
2328  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
2329 
2330  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2331  fexp(6) = 0.d0
2332  ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2333  fp(6) = -sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2334 
2335  !---------- calcul du second membre total
2336 
2337  DO j=1,6
2338  DO ni = 1, n_w
2339  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + fexp(j) + ft(j) + fp(j))
2340  ENDDO
2341  ENDDO
2342 
2343 
2344  ENDDO
2345  ENDDO
2346 
2347  END SUBROUTINE qs_00_stokes_3d
2348 
2349  SUBROUTINE qs_00_stokes_3d_new_ssr(mesh,alpha,mode,ff,V1,V2,P,dt,u0)
2350  !=================================
2351  !sans le terme de couplage
2352 
2353  USE gauss_points
2354 
2355  IMPLICIT NONE
2356  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
2357  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
2358  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !V(type,noeud) ,Vn-1
2359  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !V(type,noeud) ,Vn
2360  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
2361  REAL(KIND=8), INTENT(IN) :: dt , alpha
2362  INTEGER , INTENT(IN) :: mode
2363  INTEGER :: m, l , i , ni , j
2364  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp
2365 
2366  !fs pour le terme de forcage
2367  !fexp pour le terme de vitesse explicite
2368  !ft pour le terme temporelle
2369  !fp pour le gradient de pression
2370 
2371  TYPE(mesh_type), TARGET :: mesh
2372  INTEGER, DIMENSION(:,:), POINTER :: jj
2373  INTEGER, POINTER :: me
2374  REAL(KIND=8) :: ray
2375 
2376  CALL gauss(mesh)
2377  jj => mesh%jj
2378  me => mesh%me
2379 
2380  u0 = 0.d0
2381  fs = 0.d0
2382  ft = 0.d0
2383  fp = 0.d0
2384 
2385  DO m = 1, me
2386  DO l = 1, l_g
2387 
2388  !===Compute radius of Gauss point
2389  ray = 0
2390  DO ni = 1, n_w; i = jj(ni,m)
2391  ray = ray + mesh%rr(1,i)*ww(ni,l)
2392  END DO
2393 
2394  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
2395 
2396  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2397  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2398  fp(1) = -sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2399 
2400  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
2401 
2402  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2403  ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
2404  fp(2) = -sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2405 
2406  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
2407 
2408  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2409  ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2410  fp(3) = -sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2411 
2412  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
2413 
2414  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2415  ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2416  fp(4) = sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2417 
2418  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
2419 
2420  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2421  ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2422  fp(5) = -sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2423 
2424  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
2425 
2426  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2427  ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2428  fp(6) = -sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2429 
2430  !---------- calcul du second membre total
2431 
2432  DO j=1,6
2433  DO ni = 1, n_w
2434  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * (fs(j) + ray*(ft(j) + fp(j)))
2435  ENDDO
2436  ENDDO
2437 
2438 
2439  ENDDO
2440  ENDDO
2441 
2442  END SUBROUTINE qs_00_stokes_3d_new_ssr
2443 
2444 
2445  SUBROUTINE qs_00_stokes_3d_new(mesh,alpha,mode,ff,V1,V2,P,dt,u0)
2446  !=================================
2447  !sans le terme de couplage
2448 
2449 
2450  USE gauss_points
2451 
2452  IMPLICIT NONE
2453  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
2454  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
2455  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1 !V(type,noeud) ,Vn-1
2456  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V2 !V(type,noeud) ,Vn
2457  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
2458  REAL(KIND=8), INTENT(IN) :: dt , alpha
2459  INTEGER , INTENT(IN) :: mode
2460  INTEGER :: m, l , i , ni , j
2461  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp
2462 
2463 
2464  !fs pour le terme de forcage
2465  !fexp pour le terme de vitesse explicite
2466  !ft pour le terme temporelle
2467  !fp pour le gradient de pression
2468 
2469  TYPE(mesh_type), TARGET :: mesh
2470  INTEGER, DIMENSION(:,:), POINTER :: jj
2471  INTEGER, POINTER :: me
2472  REAL(KIND=8) :: ray
2473 
2474  CALL gauss(mesh)
2475  jj => mesh%jj
2476  me => mesh%me
2477 
2478  u0 = 0.d0
2479  fs = 0.d0
2480  ft = 0.d0
2481  fp = 0.d0
2482 
2483  DO m = 1, me
2484  DO l = 1, l_g
2485 
2486  !===Compute radius of Gauss point
2487  ray = 0
2488  DO ni = 1, n_w; i = jj(ni,m)
2489  ray = ray + mesh%rr(1,i)*ww(ni,l)
2490  END DO
2491 
2492  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
2493 
2494  fs(1) = sum(ff(1,jj(:,m)) * ww(:,l)) * rj(l,m)
2495  ft(1) = 1/(2*dt) * sum((4*v2(1,jj(:,m))-v1(1,jj(:,m))) * ww(:,l)) * rj(l,m)
2496  fp(1) = -sum(p(1,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2497 
2498  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
2499 
2500  fs(2) = sum(ff(2,jj(:,m)) * ww(:,l)) * rj(l,m)
2501  ft(2) = 1/(2*dt) * sum((4*v2(2,jj(:,m))-v1(2,jj(:,m))) * ww(:,l)) * rj(l,m)
2502  fp(2) = -sum(p(2,jj(:,m))*dw(1,:,l,m)) * rj(l,m)
2503 
2504  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
2505 
2506  ! fs(3) = SUM(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)/ray
2507  fs(3) = sum(ff(3,jj(:,m)) * ww(:,l)) * rj(l,m)
2508  ft(3) = 1/(2*dt) * sum((4*v2(3,jj(:,m))-v1(3,jj(:,m))) * ww(:,l)) * rj(l,m)
2509  fp(3) = -sum(p(2,jj(:,m))*ww(:,l)) * rj(l,m)/ray*mode
2510 
2511  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
2512 
2513  ! fs(4) = SUM(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)/ray
2514  fs(4) = sum(ff(4,jj(:,m)) * ww(:,l)) * rj(l,m)
2515  ft(4) = 1/(2*dt) * sum((4*v2(4,jj(:,m))-v1(4,jj(:,m))) * ww(:,l)) * rj(l,m)
2516  fp(4) = sum(p(1,jj(:,m))*ww(:,l)) * rj(l,m)/ray *mode
2517 
2518  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
2519 
2520  fs(5) = sum(ff(5,jj(:,m)) * ww(:,l)) * rj(l,m)
2521  ft(5) = 1/(2*dt) * sum((4*v2(5,jj(:,m))-v1(5,jj(:,m))) * ww(:,l)) * rj(l,m)
2522  fp(5) = -sum(p(1,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2523 
2524  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
2525 
2526  fs(6) = sum(ff(6,jj(:,m)) * ww(:,l)) * rj(l,m)
2527  ft(6) = 1/(2*dt) * sum((4*v2(6,jj(:,m))-v1(6,jj(:,m))) * ww(:,l)) * rj(l,m)
2528  fp(6) = -sum(p(2,jj(:,m))*dw(2,:,l,m)) * rj(l,m)
2529 
2530  !---------- calcul du second membre total
2531 
2532  DO j=1,6
2533  DO ni = 1, n_w
2534  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l) * ray * (fs(j) + ft(j) + fp(j))
2535  ENDDO
2536  ENDDO
2537 
2538 
2539  ENDDO
2540  ENDDO
2541 
2542  END SUBROUTINE qs_00_stokes_3d_new
2543 
2544 
2545  SUBROUTINE qs_00_ns_inline(mesh,mode,m_max,ff,V1,V2,P,dt,u0,meth,tps_sft)
2546  !=================================
2547  !Calcul du second membre pour le probleme de Navier_Stokes
2548  !avec terme non lineaire et source (ff)
2549  !Le terme non lineaire est calcule sur les points de Gauss
2550  !directement dans cette procedure. Le calcul peut s'effectuer avec
2551  !Matmul ou avec des multiplications classiques
2552 
2553 
2554  USE gauss_points
2555 
2556  IMPLICIT NONE
2557  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
2558  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
2559  TYPE(mesh_type), TARGET :: mesh
2560  INTEGER , INTENT(IN) :: m_max
2561  REAL(KIND=8), DIMENSION(6,mesh%np,0:m_max), INTENT(IN) :: V1 !V(type,noeud) ,Vn-1
2562  REAL(KIND=8), DIMENSION(6,mesh%np,0:m_max), INTENT(IN) :: V2 !V(type,noeud) ,Vn
2563  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
2564  REAL(KIND=8), INTENT(IN) :: dt
2565  INTEGER , INTENT(IN) :: mode
2566  CHARACTER(len=200), INTENT(IN) :: meth
2567  REAL(KIND=8), OPTIONAL, INTENT(INOUT):: tps_sft
2568 
2569  INTEGER :: m, l , i , ni , j
2570  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, fnl, smb
2571  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:,:,:),SAVE :: W, RotV, ProdI
2572  INTEGER :: N
2573  INTEGER, ALLOCATABLE,DIMENSION(:), SAVE :: j_loc
2574  REAL(KIND=8), ALLOCATABLE,DIMENSION(:,:), SAVE :: dw_loc
2575  REAL(KIND=8), DIMENSION(6,mesh%np) :: V1m, V2m, Vt
2576  REAL(KIND=8), DIMENSION(6,mesh%np,0:m_max) :: Vs
2577  LOGICAL, SAVE :: once = .true.
2578  LOGICAL, SAVE :: once_sft = .true.
2579  REAL(KIND=8):: user_time
2580  !EXTERNAL :: user_time
2581  REAL(KIND=8) :: tps, dummy, tt,tps1
2582 
2583  !fs pour le terme de forcage
2584  !ft pour le terme temporelle
2585  !fp pour le gradient de pression
2586  !fnl pour les termes non lineaires
2587 
2588  !W1, W2 : vitesses dans le format accepte par fft_prod ;
2589  ! W1(1,:)=Vr, W1(2,:)=Vth, W1(3,:)=Vz; 1<->tps n-1, 2<->tps n
2590  !RotV1, RotV2 : pareil pour le rotationnel
2591  !Prod1, Prod2 : pareil pour le produit
2592  !Prod : produits intermediaire pour le calcul du produit vect
2593 
2594  INTEGER, DIMENSION(:,:), POINTER :: jj
2595  INTEGER, POINTER :: me
2596  REAL(KIND=8) :: ray
2597 
2598 
2599  REAL(KIND=8), ALLOCATABLE, DIMENSION(:,:), SAVE :: aps, asp
2600  REAL(KIND=8) :: PI
2601  INTEGER :: li, ki
2602  REAL(KIND=8), DIMENSION(0:2*m_max) :: vr_p, vth_p, vz_p, rotr_p, rotth_p, rotz_p
2603  REAL(KIND=8), DIMENSION(3,0:2*m_max) :: prod
2604 
2605  !--------------------calcul du produit vectoriel----------------------------
2606  !{rot(V) x V}r = {rotV}z*Vth - {rotV}th*Vz = prodr
2607  !{rot(V) x V}th = {rotV}r*Vz - {rotV}z*Vr = prodth
2608  !{rot(V) x V}z = {rotV}th*Vr - {rotV}r*Vth = prodz
2609  !soit 9 sft
2610  !----------------------------------------------------------------------------
2611 
2612  CALL gauss(mesh)
2613  jj => mesh%jj
2614  me => mesh%me
2615 
2616  IF (once) THEN
2617 
2618  once =.false.
2619  !CALCUL DES MATRICES UTILES AUX SFT
2620  pi = acos(-1.d0)
2621  ALLOCATE(asp(0:2*m_max, 0:2*m_max))
2622  ALLOCATE(aps(0:2*m_max, 0:2*m_max))
2623 
2624  !calcul de asp : spectral > physique
2625 
2626  DO li= 0, 2*m_max
2627  DO ki= 0, 2*m_max
2628  IF(ki < m_max+1) THEN
2629  asp(li, ki) = cos(li*2.d0*pi/(2*m_max+1)*ki)
2630  ELSE
2631  asp(li, ki) = sin(li*2.d0*pi/(2*m_max+1)*(ki-m_max))
2632  ENDIF
2633  ENDDO
2634  ENDDO
2635 
2636  !calcul de aps : physique -> spectrale
2637 
2638  DO ki= 0, 2*m_max
2639  DO li= 0, 2*m_max
2640  IF(ki < m_max+1) THEN
2641  IF (ki == 0) THEN
2642  aps(ki,li) = 0.5d0
2643  ELSE
2644  aps(ki, li) = cos(li*2.d0*pi/(2*m_max+1)*ki)
2645  ENDIF
2646  ELSE
2647  aps(ki, li) = sin(li*2.d0*pi/(2*m_max+1)*(ki-m_max))
2648  ENDIF
2649  aps(ki, li) = aps(ki, li) * 2.d0/(float(2*m_max)+1.d0)
2650  ENDDO
2651  ENDDO
2652  WRITE(*,*) 'calcul matrices sft done'
2653 
2654  ALLOCATE(w(3,0:2*m_max,l_g,me))
2655  ALLOCATE(rotv(3,0:2*m_max,l_g,me))
2656  ALLOCATE(prodi(3,0:2*m_max,l_g,me))
2657  ALLOCATE(j_loc(mesh%gauss%n_w))
2658  ALLOCATE(dw_loc(mesh%gauss%k_d,mesh%gauss%n_w))
2659 
2660  ENDIF
2661 
2662 
2663  IF (mode == 0) THEN !pour le calcul des sft
2664  vs = 2*v2-v1
2665  ENDIF
2666 
2667  v1m = v1(:,:,mode)
2668  v2m = v2(:,:,mode)
2669 
2670  u0 = 0.d0
2671 
2672  tps = 0
2673  tps1 = 0
2674  DO m = 1, me
2675  j_loc = jj(:,m)
2676  DO l = 1, l_g
2677 
2678  !tt = user_time(dummy)
2679  dw_loc = dw(:,:,l,m)
2680 
2681  !===Compute radius of Gauss point
2682  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
2683 
2684  DO j=1, 6
2685  fs(j) = sum(ff(j,j_loc) * ww(:,l))
2686  ft(j) = 1/(2*dt)*sum((4*v2m(j,j_loc)-v1m(j,j_loc)) * ww(:,l))
2687  ENDDO
2688 
2689  fp(1) = -sum(p(1,j_loc)*dw_loc(1,:))
2690  fp(2) = -sum(p(2,j_loc)*dw_loc(1,:))
2691  fp(3) = -sum(p(2,j_loc)*ww(:,l))/ray*mode
2692  fp(4) = sum(p(1,j_loc)*ww(:,l))/ray *mode
2693  fp(5) = -sum(p(1,j_loc)*dw_loc(2,:))
2694  fp(6) = -sum(p(2,j_loc)*dw_loc(2,:))
2695 
2696 
2697  IF (mode == 0) THEN
2698 
2699  DO j = 0, m_max
2700 
2701  !tt = user_time(dummy)
2702  !-----------------rotationnel sur les points de Gauss-------------------
2703 
2704  !coeff sur les cosinus
2705  rotv(1,j,l,m) = j/ray*sum(vs(6,j_loc,j)*ww(:,l)) &
2706  -sum(vs(3,j_loc,j)*dw_loc(2,:))
2707  rotv(2,j,l,m) = sum(vs(1,j_loc,j)*dw_loc(2,:)) &
2708  -sum(vs(5,j_loc,j)*dw_loc(1,:))
2709  rotv(3,j,l,m) = 1/ray*sum(vs(3,j_loc,j)*ww(:,l))+sum(vs(3,j_loc,j) &
2710  *dw_loc(1,:))-j/ray*sum(vs(2,j_loc,j)*ww(:,l))
2711  !coeff sur les sinus
2712  IF (j /= 0) THEN
2713  rotv(1,j+m_max,l,m) = -j/ray*sum(vs(5,j_loc,j)*ww(:,l)) &
2714  -sum(vs(4,j_loc,j)*dw_loc(2,:))
2715  rotv(2,j+m_max,l,m) = sum(vs(2,j_loc,j)*dw_loc(2,:)) &
2716  -sum(vs(6,j_loc,j)*dw_loc(1,:))
2717  rotv(3,j+m_max,l,m) = 1/ray*sum(vs(4,j_loc,j)*ww(:,l))+sum(vs(4,j_loc,j) &
2718  *dw_loc(1,:))+j/ray*sum(vs(1,j_loc,j)*ww(:,l))
2719  ENDIF
2720  !-----------------vitesse sur les points de Gauss---------------------------
2721  w(1,j,l,m) = sum(vs(1,j_loc,j)*ww(:,l))
2722  w(2,j,l,m) = sum(vs(3,j_loc,j)*ww(:,l))
2723  w(3,j,l,m) = sum(vs(5,j_loc,j)*ww(:,l))
2724  IF (j /= 0) THEN
2725  w(1,j+m_max,l,m) = sum(vs(2,j_loc,j)*ww(:,l))
2726  w(2,j+m_max,l,m) = sum(vs(4,j_loc,j)*ww(:,l))
2727  w(3,j+m_max,l,m) = sum(vs(6,j_loc,j)*ww(:,l))
2728  ENDIF
2729 
2730  !tps1 = tps1+user_time(dummy) - tt
2731  ENDDO
2732 
2733  tt = user_time(dummy)
2734 
2735  IF (meth == 'sum') THEN
2736  DO li=0, 2*m_max
2737  vr_p(li) = sum(asp(li,:)*w(1,:,l,m))
2738  vth_p(li) = sum(asp(li,:)*w(2,:,l,m))
2739  vz_p(li) = sum(asp(li,:)*w(3,:,l,m))
2740  rotr_p(li) = sum(asp(li,:)*rotv(1,:,l,m))
2741  rotth_p(li) = sum(asp(li,:)*rotv(2,:,l,m))
2742  rotz_p(li) = sum(asp(li,:)*rotv(3,:,l,m))
2743  ENDDO
2744  ELSEIF (meth == 'matmul') THEN
2745  vr_p = matmul(asp,w(1,:,l,m) )
2746  vth_p = matmul(asp,w(2,:,l,m))
2747  vz_p = matmul(asp,w(3,:,l,m))
2748 
2749  rotr_p = matmul(asp,rotv(1,:,l,m))
2750  rotth_p = matmul(asp,rotv(2,:,l,m))
2751  rotz_p = matmul(asp,rotv(3,:,l,m))
2752  ENDIF
2753 
2754  prod(1,:) = rotz_p(:) * vth_p(:) - rotth_p(:) * vz_p(:)
2755  prod(2,:) = rotr_p(:) * vz_p(:) - rotz_p(:) * vr_p(:)
2756  prod(3,:) = rotth_p(:) * vr_p(:) - rotr_p(:) * vth_p(:)
2757 
2758  IF (meth == 'sum') THEN
2759  DO li=0, 2*m_max
2760  prodi(1,li,l,m) = sum(aps(li,:)*prod(1,:))
2761  prodi(2,li,l,m) = sum(aps(li,:)*prod(2,:))
2762  prodi(3,li,l,m) = sum(aps(li,:)*prod(3,:))
2763  ENDDO
2764  ELSEIF (meth == 'matmul') THEN
2765  prodi(1,:,l,m) = matmul(aps, prod(1,:))
2766  prodi(2,:,l,m) = matmul(aps, prod(2,:))
2767  prodi(3,:,l,m) = matmul(aps, prod(3,:))
2768  ENDIF
2769 
2770  IF (PRESENT(tps_sft)) THEN
2771  tps_sft = tps_sft + user_time(dummy) - tt
2772  ENDIF
2773 
2774  ENDIF
2775 
2776  !-----------------assemblage du second membre du terme NL--------------
2777 
2778  fnl(1) = prodi(1,mode,l,m)
2779  fnl(3) = prodi(2,mode,l,m)
2780  fnl(5) = prodi(3,mode,l,m)
2781  IF (mode /= 0) THEN
2782  fnl(2) = prodi(1,mode+m_max,l,m)
2783  fnl(4) = prodi(2,mode+m_max,l,m)
2784  fnl(6) = prodi(3,mode+m_max,l,m)
2785  ENDIF
2786 
2787  smb = (fnl+ft+fp+fs)*ray*rj(l,m)
2788 
2789  !---------- calcul du second membre total
2790 
2791  DO j=1,6
2792  DO ni = 1, n_w
2793  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l)*smb(j)
2794  ENDDO
2795  ENDDO
2796 
2797  ENDDO
2798  ENDDO
2799 
2800  !write(*,*) 'temps d une boucle ss terme nl=',tps
2801  !write(*,*) 'temps du calcul de V et RotV aux points de Gauss=',tps1
2802 
2803  END SUBROUTINE qs_00_ns_inline
2804 
2805 
2806  SUBROUTINE qs_stokes(mesh,mode,ff,V1m,V2m,P,dt,u0)
2807  !=================================
2808  !Calcul du second membre pour le probleme de Stokes
2809  !sans terme non lineaire et source (ff)
2810 
2811  USE gauss_points
2812 
2813  IMPLICIT NONE
2814  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
2815  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
2816  TYPE(mesh_type), TARGET :: mesh
2817  REAL(KIND=8), DIMENSION(6,mesh%np), INTENT(IN) :: V1m, V2m !V(type,noeud)
2818  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
2819  REAL(KIND=8), INTENT(IN) :: dt
2820  INTEGER , INTENT(IN) :: mode
2821 
2822  INTEGER :: m, l , i , ni , j
2823  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, fnl, smb
2824  INTEGER :: N
2825  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
2826  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
2827  LOGICAL, SAVE :: once = .true.
2828  LOGICAL, SAVE :: once_sft = .true.
2829 
2830  !fs pour le terme de forcage
2831  !ft pour le terme temporelle
2832  !fp pour le gradient de pression
2833 
2834  INTEGER, DIMENSION(:,:), POINTER :: jj
2835  INTEGER, POINTER :: me
2836  REAL(KIND=8) :: ray
2837 
2838  CALL gauss(mesh)
2839  jj => mesh%jj
2840  me => mesh%me
2841 
2842  u0 = 0.d0
2843 
2844  DO m = 1, me
2845  j_loc = jj(:,m)
2846  DO l = 1, l_g
2847  dw_loc = dw(:,:,l,m)
2848 
2849  !===Compute radius of Gauss point
2850  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
2851 
2852  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
2853  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
2854 
2855  fs(1) = sum(ff(1,j_loc) * ww(:,l))
2856  ft(1) = 1/(2*dt) * sum((4*v2m(1,j_loc)-v1m(1,j_loc)) * ww(:,l))
2857  fp(1) = -sum(p(1,j_loc)*dw_loc(1,:))
2858 
2859  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
2860 
2861  fs(2) = sum(ff(2,j_loc) * ww(:,l)) !* rj(l,m)
2862  ft(2) = 1/(2*dt) * sum((4*v2m(2,j_loc)-v1m(2,j_loc)) * ww(:,l)) !* rj(l,m)
2863  fp(2) = -sum(p(2,j_loc)*dw_loc(1,:)) !* rj(l,m)
2864 
2865  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
2866 
2867  fs(3) = sum(ff(3,j_loc) * ww(:,l)) !* rj(l,m)
2868  ft(3) = 1/(2*dt) * sum((4*v2m(3,j_loc)-v1m(3,j_loc)) * ww(:,l)) !* rj(l,m)
2869  fp(3) = -sum(p(2,j_loc)*ww(:,l))/ray*mode !* rj(l,m)
2870 
2871  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
2872 
2873  fs(4) = sum(ff(4,j_loc) * ww(:,l)) !* rj(l,m)
2874  ft(4) = 1/(2*dt) * sum((4*v2m(4,j_loc)-v1m(4,j_loc)) * ww(:,l)) !* rj(l,m)
2875  fp(4) = sum(p(1,j_loc)*ww(:,l))/ray *mode !* rj(l,m)
2876 
2877  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
2878 
2879  fs(5) = sum(ff(5,j_loc) * ww(:,l)) !* rj(l,m)
2880  ft(5) = 1/(2*dt) * sum((4*v2m(5,j_loc)-v1m(5,j_loc)) * ww(:,l)) !* rj(l,m)
2881  fp(5) = -sum(p(1,j_loc)*dw_loc(2,:)) !* rj(l,m)
2882 
2883  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
2884 
2885  fs(6) = sum(ff(6,j_loc) * ww(:,l)) !* rj(l,m)
2886  ft(6) = 1/(2*dt) * sum((4*v2m(6,j_loc)-v1m(6,j_loc)) * ww(:,l)) !* rj(l,m)
2887  fp(6) = -sum(p(2,j_loc)*dw_loc(2,:)) !* rj(l,m)
2888 
2889  !-------calcul du second membre pour le terme nonlineaire------------------------
2890  !-------------------------------------------------------------------------------
2891 
2892  smb = (ft+fp+fs)*ray*rj(l,m)
2893  !write(95,*) m,l,ft(3)
2894  !---------- calcul du second membre total
2895  !smb = 1.d0
2896  DO j=1,6
2897  DO ni = 1, n_w
2898  u0(j,jj(ni,m)) = u0(j,jj(ni,m)) + ww(ni,l)*smb(j)
2899  ENDDO
2900  ENDDO
2901 
2902  ENDDO
2903  ENDDO
2904 
2905 
2906  END SUBROUTINE qs_stokes
2907 
2908  SUBROUTINE qs_navier_stokes_2006(mesh,mode,ff,V1m,P,dt,u0,rotv_v)
2909  !=================================
2910  !Calcul du second membre pour le probleme de Stokes
2911  !sans terme non lineaire et source (ff)
2912 
2913  USE gauss_points
2914 
2915  IMPLICIT NONE
2916  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotv_v !analytique (type,noeud)
2917  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT):: u0 !u0(type,noeud)
2918  TYPE(mesh_type), TARGET :: mesh
2919  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1m !V(type,noeud)
2920  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
2921  REAL(KIND=8), INTENT(IN) :: dt
2922  INTEGER, INTENT(IN) :: mode
2923  INTEGER :: m, l , i , ni , j, index
2924  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb
2925  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
2926  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
2927 
2928  !fs pour le terme de forcage
2929  !ft pour le terme temporelle
2930  !fp pour le gradient de pression
2931 
2932  INTEGER, DIMENSION(:,:), POINTER :: jj
2933  INTEGER, POINTER :: me
2934  REAL(KIND=8) :: ray
2935 
2936  CALL gauss(mesh)
2937  jj => mesh%jj
2938  me => mesh%me
2939 
2940  u0 = 0.d0
2941 
2942 
2943  index = 0
2944  DO m = 1, me
2945  j_loc = jj(:,m)
2946  DO l = 1, l_g
2947  index = index +1
2948  dw_loc = dw(:,:,l,m)
2949 
2950  !===Compute radius of Gauss point
2951  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
2952 
2953  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
2954  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
2955 
2956  fs(1) = sum(ff(j_loc,1) * ww(:,l))
2957  ft(1) = sum(v1m(j_loc,1) * ww(:,l))
2958  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
2959 
2960  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
2961 
2962  fs(2) = sum(ff(j_loc,2) * ww(:,l)) !* rj(l,m)
2963  ft(2) = sum(v1m(j_loc,2) * ww(:,l)) !* rj(l,m)
2964  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:)) !* rj(l,m)
2965 
2966  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
2967 
2968  fs(3) = sum(ff(j_loc,3) * ww(:,l)) !* rj(l,m)
2969  ft(3) = sum(v1m(j_loc,3) * ww(:,l)) !* rj(l,m)
2970  fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode !* rj(l,m)
2971 
2972  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
2973 
2974  fs(4) = sum(ff(j_loc,4) * ww(:,l)) !* rj(l,m)
2975  ft(4) = sum(v1m(j_loc,4) * ww(:,l)) !* rj(l,m)
2976  fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode !* rj(l,m)
2977 
2978  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
2979 
2980  fs(5) = sum(ff(j_loc,5) * ww(:,l)) !* rj(l,m)
2981  ft(5) = sum(v1m(j_loc,5) * ww(:,l)) !* rj(l,m)
2982  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:)) !* rj(l,m)
2983 
2984  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
2985 
2986  fs(6) = sum(ff(j_loc,6) * ww(:,l)) !* rj(l,m)
2987  ft(6) = sum(v1m(j_loc,6) * ww(:,l)) !* rj(l,m)
2988  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:)) !* rj(l,m)
2989 
2990  !-------calcul du second membre pour le terme nonlineaire------------------------
2991  !-------------------------------------------------------------------------------
2992 
2993  smb = (ft+fp+fs-rotv_v(:,index))*ray*rj(l,m)
2994 
2995  !---------- calcul du second membre total
2996  DO j=1,6
2997  DO ni = 1, n_w
2998  u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j)
2999  ENDDO
3000  ENDDO
3001 
3002  ENDDO
3003  ENDDO
3004 
3005 
3006  END SUBROUTINE qs_navier_stokes_2006
3007 
3008  SUBROUTINE qs_ns_2006(mesh,mode,ff,V1m,P,dt,u0,rotv_v)
3009  !=================================
3010  !Calcul du second membre pour le probleme de Stokes
3011  !sans terme non lineaire et source (ff)
3012 
3013  USE gauss_points
3014 
3015  IMPLICIT NONE
3016  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotv_v !analytique (type,noeud)
3017  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT):: u0 !u0(type,noeud)
3018  TYPE(mesh_type), TARGET :: mesh
3019  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1m !V(type,noeud)
3020  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
3021  REAL(KIND=8), INTENT(IN) :: dt
3022  INTEGER, INTENT(IN) :: mode
3023  INTEGER :: m, l , i , ni , j, index
3024  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb
3025  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
3026  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3027 
3028  !fs pour le terme de forcage
3029  !ft pour le terme temporelle
3030  !fp pour le gradient de pression
3031 
3032  INTEGER, DIMENSION(:,:), POINTER :: jj
3033  INTEGER, POINTER :: me
3034  REAL(KIND=8) :: ray
3035 
3036  CALL gauss(mesh)
3037  jj => mesh%jj
3038  me => mesh%me
3039 
3040  u0 = 0.d0
3041 
3042 
3043  index = 0
3044  DO m = 1, me
3045  j_loc = jj(:,m)
3046  DO l = 1, l_g
3047  index = index +1
3048  dw_loc = dw(:,:,l,m)
3049 
3050  !===Compute radius of Gauss point
3051  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3052 
3053  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
3054  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
3055 
3056  fs(1) = sum(ff(j_loc,1) * ww(:,l))
3057  ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3058  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3059 
3060  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
3061 
3062  fs(2) = sum(ff(j_loc,2) * ww(:,l)) !* rj(l,m)
3063  ft(2) = sum(v1m(j_loc,2) * ww(:,l)) !* rj(l,m)
3064  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:)) !* rj(l,m)
3065 
3066  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
3067 
3068  fs(3) = sum(ff(j_loc,3) * ww(:,l)) !* rj(l,m)
3069  ft(3) = sum(v1m(j_loc,3) * ww(:,l)) !* rj(l,m)
3070  fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode !* rj(l,m)
3071 
3072  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
3073 
3074  fs(4) = sum(ff(j_loc,4) * ww(:,l)) !* rj(l,m)
3075  ft(4) = sum(v1m(j_loc,4) * ww(:,l)) !* rj(l,m)
3076  fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode !* rj(l,m)
3077 
3078  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
3079 
3080  fs(5) = sum(ff(j_loc,5) * ww(:,l)) !* rj(l,m)
3081  ft(5) = sum(v1m(j_loc,5) * ww(:,l)) !* rj(l,m)
3082  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:)) !* rj(l,m)
3083 
3084  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
3085 
3086  fs(6) = sum(ff(j_loc,6) * ww(:,l)) !* rj(l,m)
3087  ft(6) = sum(v1m(j_loc,6) * ww(:,l)) !* rj(l,m)
3088  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:)) !* rj(l,m)
3089 
3090  !-------calcul du second membre pour le terme nonlineaire------------------------
3091  !-------------------------------------------------------------------------------
3092 
3093  smb = (ft+fp+fs-rotv_v(index,:))*ray*rj(l,m)
3094 
3095  !---------- calcul du second membre total
3096  DO j=1,6
3097  DO ni = 1, n_w
3098  u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j)
3099  ENDDO
3100  ENDDO
3101 
3102  ENDDO
3103  ENDDO
3104 
3105 
3106  END SUBROUTINE qs_ns_2006
3107 
3108  SUBROUTINE qs_ns_stab_new(mesh,mode,ff,vel_tot,V1m,vit,P,dudt,phalf,nlhalf,dt,u0,rotv_v)
3109  !=================================
3110  !Calcul du second membre pour le probleme de Stokes
3111  !sans terme non lineaire et source (ff)
3112 
3113  USE gauss_points
3114  USE chaine_caractere
3115  USE sub_plot
3116  USE input_data
3117  IMPLICIT NONE
3118  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotv_v, dudt, nlhalf
3119  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: vel_tot !(noeud)
3120  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT):: u0 !u0(noeud, type)
3121  TYPE(mesh_type), TARGET :: mesh
3122  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1m, vit !V(noeud, type)
3123  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P, phalf
3124  REAL(KIND=8), INTENT(IN) :: dt
3125  INTEGER, INTENT(IN) :: mode
3126  INTEGER :: m, l , i , ni , j, index, index2, type, k
3127  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb, fv, mult
3128  REAL(KIND=8), DIMENSION(2,6,mesh%gauss%l_G) :: grad
3129  REAL(KIND=8), DIMENSION(6,mesh%gauss%l_G) :: vitloc
3130  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
3131  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3132  REAL(KIND=8), DIMENSION(mesh%me) :: visc_plot
3133  !fs pour le terme de forcage
3134  !ft pour le terme temporelle
3135  !fp pour le gradient de pression
3136 
3137  INTEGER, DIMENSION(:,:), POINTER :: jj
3138  INTEGER, POINTER :: me
3139  REAL(KIND=8), DIMENSION(6) :: visc1, visc2
3140  REAL(KIND=8) :: ray, div1, div2, hloc, cfl, vloc, normal_vit
3141  REAL(KIND=8), SAVE :: coeff_ed_st, R_eff, visc_eff, surf, nu_loc, coeff_visc_ordre_un
3142  LOGICAL, SAVE :: once = .true.
3143  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: h
3144 
3145  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%l_Gs, 2) :: dwni_loc
3146  INTEGER, DIMENSION(mesh%gauss%n_w,2) :: jji_loc
3147  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: u0loci, uloci
3148  REAL(KIND=8), DIMENSION(6,mesh%me) :: viscosity
3149  REAL(KIND=8), DIMENSION(6) :: norm_vit
3150  REAL(KIND=8) :: dul, h2, type_fe, ed_st
3151  INTEGER :: ms, ls, cotei
3152 
3153  CALL gauss(mesh)
3154  jj => mesh%jj
3155  me => mesh%me
3156 
3157  IF (once) THEN
3158  once =.false.
3159 
3160  surf = 0.d0
3161  DO m = 1, me
3162  DO l = 1, l_g
3163  ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
3164  surf = surf + rj(l,m)*ray
3165  END DO
3166  END DO
3167 
3168  IF (mesh%gauss%n_w==3) THEN
3169  type_fe = 1
3170  ELSE
3171  type_fe = 2
3172  END IF
3173  coeff_ed_st = inputs%coeff3*0.02d0/type_fe
3174  coeff_visc_ordre_un = inputs%coeff4*0.25d0
3175  ALLOCATE(h(mesh%mi))
3176  DO ms = 1, mesh%mi
3177  h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
3178  END DO
3179  END IF
3180 
3181  !ATTENTION: JLG Jan 25 2010
3182  !ATTENTION: inputs%coeff1 is assumed to be of the order of the convective velocity
3183  !that simplifies the semi-implicit treatment of the LES viscosity
3184  normal_vit = maxval(vel_tot)
3185  DO type = 1, 6
3186  norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
3187  END DO
3188  !ATTENTION: JLG Jan 25 2010
3189  u0 = 0.d0
3190  r_eff = 0.d0
3191  cfl = 0
3192  index = 0
3193  index2 = 0
3194  IF (inputs%LES) THEN
3195  DO m = 1, me
3196  j_loc = jj(:,m)
3197  hloc = sqrt(sum(rj(:,m)))
3198  vloc = maxval(vel_tot(j_loc))
3199  cfl = max(vloc*dt/hloc,cfl)
3200  visc1 = 0
3201  visc2 = 0
3202  DO l = 1, l_g
3203  index2 = index2 +1
3204  dw_loc = dw(:,:,l,m)
3205 
3206  !===Compute radius of Gauss point
3207  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3208 
3209  !--------calcul de la premiere composante
3210  ft(1) = sum(dudt(j_loc,1) * ww(:,l))
3211  fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
3212  !--------calcul de la seconde composante
3213  ft(2) = sum(dudt(j_loc,2) * ww(:,l))
3214  fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
3215  !--------calcul de la troisieme composante
3216  ft(3) = sum(dudt(j_loc,3) * ww(:,l))
3217  fp(3) = sum(phalf(j_loc,2)*ww(:,l))*mode/ray
3218  !--------calcul de la quatrieme composante
3219  ft(4) = sum(dudt(j_loc,4) * ww(:,l))
3220  fp(4) = -sum(phalf(j_loc,1)*ww(:,l))*mode/ray
3221  !--------calcul de la cinquieme composante
3222  ft(5) = sum(dudt(j_loc,5) * ww(:,l))
3223  fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
3224  !--------calcul de la sixieme composante
3225  ft(6) = sum(dudt(j_loc,6) * ww(:,l))
3226  fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
3227  !-------calcul du second membre pour le terme nonlineaire------------------------
3228 
3229  visc1= max(visc1,abs(ft+fp+nlhalf(index2,:)))
3230 
3231  !--------Calcul du gradient de la vitesse
3232  DO type = 1, 6
3233  DO k = 1 ,2
3234  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3235  END DO
3236  vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3237  END DO
3238 
3239  !--------Calcul de la divergence
3240  div1 = abs(grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray)
3241  div2 = abs(grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray)
3242  visc2(1) = max(visc2(1),div1)
3243  visc2(2) = max(visc2(1),div2)
3244  END DO
3245  visc2(4) = visc2(1); visc2(5) = visc2(1)
3246  visc2(3) = visc2(2); visc2(6) = visc2(2)
3247 
3248  nu_loc = 0.d0
3249  DO type = 1, 6
3250  !visc1(type) = MAX(visc1(type)/normal_vit,visc2(type))
3251  !visc1(type) = MAX(visc1(type),visc2(type)*normal_vit)/MAXVAL(ABS(vitloc(type,:)))
3252  visc1(type) = max(visc1(type),2*visc2(type)*normal_vit)/norm_vit(type)
3253  visc_eff = hloc*min(coeff_visc_ordre_un*normal_vit,inputs%coeff2*hloc*visc1(type))
3254  nu_loc = nu_loc + visc_eff
3255  !======Semi-implicit version==========
3256  viscosity(type,m) = inputs%coeff1*hloc - visc_eff
3257  !======Semi-implicit version==========
3258  END DO
3259  r_eff = r_eff + (nu_loc + dt**2*hloc)*hloc**2*ray
3260  visc_plot(m) = (nu_loc/6)/(coeff_visc_ordre_un*hloc*normal_vit)
3261  END DO
3262  ELSE
3263  viscosity = 0.d0
3264  DO m = 1, me
3265  j_loc = jj(:,m)
3266  hloc = sqrt(sum(rj(:,m)))
3267  vloc = maxval(vel_tot(j_loc))
3268  cfl = max(vloc*dt/hloc,cfl)
3269  END DO
3270  END IF
3271  !DO type = 1, 6
3272  ! CALL average(mesh,viscosity(type,:))
3273  !END DO
3274 
3275  DO m = 1, me
3276  mult(:)= viscosity(:,m)
3277  j_loc = jj(:,m)
3278  DO l = 1, l_g
3279  index = index +1
3280  dw_loc = dw(:,:,l,m)
3281  DO type = 1, 6
3282  DO k = 1 ,2
3283  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3284  END DO
3285  vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3286  END DO
3287 
3288  !===Compute radius of Gauss point
3289  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3290 
3291  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
3292  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
3293  fs(1) = sum(ff(j_loc,1) * ww(:,l))
3294  ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3295  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3296  fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
3297  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
3298  fs(2) = sum(ff(j_loc,2) * ww(:,l))
3299  ft(2) = sum(v1m(j_loc,2) * ww(:,l))
3300  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
3301  fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
3302  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
3303  fs(3) = sum(ff(j_loc,3) * ww(:,l))
3304  ft(3) = sum(v1m(j_loc,3) * ww(:,l))
3305  fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
3306  fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
3307  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
3308  fs(4) = sum(ff(j_loc,4) * ww(:,l))
3309  ft(4) = sum(v1m(j_loc,4) * ww(:,l))
3310  fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
3311  fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
3312  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
3313  fs(5) = sum(ff(j_loc,5) * ww(:,l))
3314  ft(5) = sum(v1m(j_loc,5) * ww(:,l))
3315  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
3316  fv(5) = vitloc(5,l)*(mode/ray)**2
3317  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
3318  fs(6) = sum(ff(j_loc,6) * ww(:,l))
3319  ft(6) = sum(v1m(j_loc,6) * ww(:,l))
3320  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
3321  fv(6) = vitloc(6,l)*(mode/ray)**2
3322  !-------calcul du second membre pour le terme nonlineaire------------------------
3323  !-------------------------------------------------------------------------------
3324 
3325  fv = mult*fv
3326 
3327  smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*rj(l,m)
3328  DO type = 1, 6
3329  grad(:,type,l) = mult(type)*grad(:,type,l)*ray*rj(l,m)
3330  END DO
3331  !---------- calcul du second membre total
3332  DO j=1,6
3333  DO ni = 1, n_w
3334  u0(j_loc(ni),j) = u0(j_loc(ni),j) + ww(ni,l)*smb(j) + sum(dw_loc(:,ni)*grad(:,j,l))
3335  ENDDO
3336  ENDDO
3337 
3338  ENDDO
3339  ENDDO
3340 
3341 
3342  WRITE(*,*) ' CFL = ', cfl, ' inputs%LES ', inputs%LES
3343  IF (inputs%LES) THEN
3344  WRITE(*,*) ' R_eff for mode', mode, normal_vit*6*surf/r_eff
3345  END IF
3346  !IF (mode==0) THEN
3347  !CALL plot_const_p1_label(mesh%jj, mesh%rr, visc_plot, 'tttt_0.plt')
3348  !END IF
3349 
3350  IF (inputs%LES) THEN
3351  ed_st = normal_vit*coeff_ed_st
3352  DO ms = 1, mesh%mi
3353  dwni_loc = mesh%gauss%dwni(:,:,:,ms)
3354  jji_loc = mesh%jji(:,:,ms)
3355  h2 = -ed_st*h(ms)**2
3356 
3357  j_loc(1:n_ws) = mesh%jjsi(1:n_ws,ms)
3358  DO type = 1, 6
3359  DO cotei = 1, 2
3360  uloci(:,cotei) = vit(jji_loc(:,cotei),type)
3361  END DO
3362 
3363  u0loci = 0.d0
3364  DO ls = 1, mesh%gauss%l_Gs
3365  !===Compute radius of Gauss point at face center
3366  ray = mesh%rr(1,j_loc(3))
3367 
3368  !--------Calcul du saut de la derivee normale de la vitesse
3369  dul = sum(dwni_loc(:,ls,:)*uloci)*mesh%gauss%rji(ls,ms)*h2*ray
3370  DO cotei = 1, 2
3371  DO ni = 1, mesh%gauss%n_w
3372  u0loci(ni, cotei) = u0loci(ni, cotei) + dwni_loc(ni,ls,cotei)*dul
3373  END DO
3374  END DO
3375  END DO
3376 
3377  DO cotei = 1, 2
3378  DO ni = 1, mesh%gauss%n_w
3379  u0(jji_loc(ni,cotei),type) = u0(jji_loc(ni,cotei),type) + u0loci(ni,cotei)
3380  END DO
3381  END DO
3382  END DO
3383 
3384  END DO
3385  END IF
3386 
3387  END SUBROUTINE qs_ns_stab_new
3388 
3389  SUBROUTINE qs_ns_stab_2010(mesh, pp_mesh, mode,ff,vel_tot,V1m,vit,P,dudt,phalf,nlhalf,dt,u0,rotv_v)
3390  !=================================
3391  !Calcul du second membre pour le probleme de Stokes
3392  !sans terme non lineaire et source (ff)
3393 
3394  USE gauss_points
3395  USE chaine_caractere
3396  USE sub_plot
3397  IMPLICIT NONE
3398  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotv_v, dudt, nlhalf
3399  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: vel_tot !(noeud)
3400  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT):: u0 !u0(noeud, type)
3401  TYPE(mesh_type), TARGET :: mesh, pp_mesh
3402  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1m, vit !V(noeud, type)
3403  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P, phalf
3404  REAL(KIND=8), INTENT(IN) :: dt
3405  INTEGER, INTENT(IN) :: mode
3406  INTEGER :: m, l , i , ni , j, index, index2, type, k
3407  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb, fv, mult
3408  REAL(KIND=8), DIMENSION(2,6,mesh%gauss%l_G) :: grad
3409  REAL(KIND=8), DIMENSION(6,mesh%gauss%l_G) :: vitloc
3410  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
3411  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3412  REAL(KIND=8), DIMENSION(mesh%me) :: visc_plot
3413  !fs pour le terme de forcage
3414  !ft pour le terme temporelle
3415  !fp pour le gradient de pression
3416 
3417  INTEGER, DIMENSION(:,:), POINTER :: jj
3418  INTEGER, POINTER :: me
3419  REAL(KIND=8), DIMENSION(6) :: visc1, visc2
3420  REAL(KIND=8) :: ray, div1, div2, hloc, cfl, vloc, normal_vit
3421  REAL(KIND=8), SAVE :: coeff_ed_st, R_eff, visc_eff, surf, nu_loc
3422  LOGICAL, SAVE :: once = .true.
3423  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: h
3424 
3425  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%l_Gs, 2) :: dwni_loc
3426  INTEGER, DIMENSION(mesh%gauss%n_w,2) :: jji_loc
3427  REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: u0loci, uloci
3428  REAL(KIND=8), DIMENSION(6,mesh%me) :: viscosity
3429  REAL(KIND=8), DIMENSION(6,pp_mesh%np) :: viscosity_clement
3430  REAL(KIND=8), DIMENSION(mesh%np,6) :: visc_p2
3431  REAL(KIND=8), DIMENSION(6) :: norm_vit
3432  REAL(KIND=8) :: dul, h2, type_fe, ed_st
3433  INTEGER :: ms, ls, cotei
3434 
3435  CALL gauss(mesh)
3436  jj => mesh%jj
3437  me => mesh%me
3438 
3439  IF (once) THEN
3440  once =.false.
3441  surf = sum(rj(:,:))
3442 
3443  IF (mesh%gauss%n_w==3) THEN
3444  type_fe = 1
3445  ELSE
3446  type_fe = 2
3447  END IF
3448  coeff_ed_st = 0.05d0/type_fe
3449 
3450  ALLOCATE(h(mesh%mi))
3451  DO ms = 1, mesh%mi
3452  h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
3453  END DO
3454  END IF
3455 
3456  !ATTENTION: JLG Jan 25 2010
3457  !ATTENTION: inputs%coeff1 is assumed to be of the order of the convective velocity
3458  !that simplifies the semi-implicit treatment of the LES viscosity
3459  normal_vit = maxval(vel_tot)
3460  DO type = 1, 6
3461  norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
3462  END DO
3463  !ATTENTION: JLG Jan 25 2010
3464  u0 = 0.d0
3465  r_eff = 0.d0
3466  cfl = 0
3467  index = 0
3468  index2 = 0
3469  DO m = 1, me
3470  j_loc = jj(:,m)
3471  hloc = sqrt(sum(rj(:,m)))
3472  !vloc = MIN(maxval(vel_tot(j_loc)),1.d0)
3473  vloc = maxval(vel_tot(j_loc))
3474  cfl = max(vloc*dt/hloc,cfl)
3475  visc1 = 0
3476  visc2 = 0
3477  DO l = 1, l_g
3478  index2 = index2 +1
3479  dw_loc = dw(:,:,l,m)
3480 
3481  !===Compute radius of Gauss point
3482  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3483 
3484  !--------calcul de la premiere composante
3485  ft(1) = sum(dudt(j_loc,1) * ww(:,l))
3486  fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
3487  !--------calcul de la seconde composante
3488  ft(2) = sum(dudt(j_loc,2) * ww(:,l))
3489  fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
3490  !--------calcul de la troisieme composante
3491  ft(3) = sum(dudt(j_loc,3) * ww(:,l))
3492  fp(3) = sum(phalf(j_loc,2)*ww(:,l))*mode/ray
3493  !--------calcul de la quatrieme composante
3494  ft(4) = sum(dudt(j_loc,4) * ww(:,l))
3495  fp(4) = -sum(phalf(j_loc,1)*ww(:,l))*mode/ray
3496  !--------calcul de la cinquieme composante
3497  ft(5) = sum(dudt(j_loc,5) * ww(:,l))
3498  fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
3499  !--------calcul de la sixieme composante
3500  ft(6) = sum(dudt(j_loc,6) * ww(:,l))
3501  fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
3502  !-------calcul du second membre pour le terme nonlineaire------------------------
3503 
3504  visc1= max(visc1,abs(ft+fp+nlhalf(index2,:)))
3505 
3506  !--------Calcul du gradient de la vitesse
3507  DO type = 1, 6
3508  DO k = 1 ,2
3509  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3510  END DO
3511  vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3512  END DO
3513 
3514  !--------Calcul de la divergence
3515  div1 = abs(grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray)
3516  div2 = abs(grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray)
3517  visc2(1) = div1; visc2(4) = div1; visc2(5) = div1
3518  visc2(2) = div2; visc2(3) = div2; visc2(6) = div2
3519  END DO
3520 
3521  nu_loc = 0.d0
3522  DO type = 1, 6
3523  !visc1(type) = MAX(visc1(type)/normal_vit,visc2(type))
3524  visc1(type) = max(visc1(type),visc2(type)*normal_vit)/norm_vit(type)
3525  !visc1(type) = visc1(type)/norm_vit(type)
3526  visc_eff = hloc*min(0.25d0*normal_vit,inputs%coeff2*hloc*visc1(type))
3527  nu_loc = nu_loc + visc_eff
3528  !======Semi-implicit version==========
3529  viscosity(type,m) = inputs%coeff1*hloc - visc_eff
3530  !visc1(type) = inputs%coeff1*hloc - visc_eff
3531  !======Semi-implicit version==========
3532  !======Totally explicit version=======
3533  !visc1(type) = -visc_eff
3534  !======Totally explicit version=======
3535  END DO
3536  r_eff = r_eff + (nu_loc + dt**2*hloc)*hloc**2
3537  visc_plot(m) = (nu_loc/6)/(0.25*hloc*normal_vit)
3538  END DO
3539 
3540  DO type = 1, 6
3541  !CALL average(mesh,viscosity(type,:))
3542  CALL clement_c(pp_mesh,viscosity(type,:),viscosity_clement(type,:))
3543  CALL inject_clement(pp_mesh%jj, mesh%jj, viscosity_clement(type,:), visc_p2(:,type))
3544  END DO
3545 !CALL plot_const_p1_label(mesh%jj, mesh%rr, viscosity(1,:), 'tttt_0.plt')
3546 !CALL plot_scalar_field(pp_mesh%jj, pp_mesh%rr, viscosity_clement(1,:), 'v1_p1.plt')
3547 !CALL plot_scalar_field(mesh%jj, mesh%rr, visc_p2(:,1)/(0.25*hloc*normal_vit), 'v1.plt')
3548 !CALL plot_scalar_field(mesh%jj, mesh%rr, visc_p2(:,3)/(0.25*hloc*normal_vit), 'v3.plt')
3549 !CALL plot_scalar_field(mesh%jj, mesh%rr, visc_p2(:,5)/(0.25*hloc*normal_vit), 'v5.plt')
3550 !stop
3551  DO m = 1, me
3552  !mult(:)= viscosity(:,m)
3553  j_loc = jj(:,m)
3554  DO l = 1, l_g
3555 
3556  index = index +1
3557  dw_loc = dw(:,:,l,m)
3558  DO type = 1, 6
3559  mult(type)= sum(visc_p2(j_loc(:),type)*ww(:,l))
3560  DO k = 1 ,2
3561  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3562  END DO
3563  vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3564  END DO
3565 
3566  !===Compute radius of Gauss point
3567  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3568 
3569  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
3570  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
3571  fs(1) = sum(ff(j_loc,1) * ww(:,l))
3572  ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3573  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3574  fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
3575  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
3576  fs(2) = sum(ff(j_loc,2) * ww(:,l))
3577  ft(2) = sum(v1m(j_loc,2) * ww(:,l))
3578  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
3579  fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
3580  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
3581  fs(3) = sum(ff(j_loc,3) * ww(:,l))
3582  ft(3) = sum(v1m(j_loc,3) * ww(:,l))
3583  fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode
3584  fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
3585  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
3586  fs(4) = sum(ff(j_loc,4) * ww(:,l))
3587  ft(4) = sum(v1m(j_loc,4) * ww(:,l))
3588  fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode
3589  fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
3590  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
3591  fs(5) = sum(ff(j_loc,5) * ww(:,l))
3592  ft(5) = sum(v1m(j_loc,5) * ww(:,l))
3593  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
3594  fv(5) = vitloc(5,l)*(mode/ray)**2
3595  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
3596  fs(6) = sum(ff(j_loc,6) * ww(:,l))
3597  ft(6) = sum(v1m(j_loc,6) * ww(:,l))
3598  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
3599  fv(6) = vitloc(6,l)*(mode/ray)**2
3600  !-------calcul du second membre pour le terme nonlineaire------------------------
3601  !-------------------------------------------------------------------------------
3602 
3603  !mult = visc1
3604  fv = mult*fv
3605 
3606  smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*rj(l,m)
3607  DO type = 1, 6
3608  grad(:,type,l) = mult(type)*grad(:,type,l)*ray*rj(l,m)
3609  END DO
3610  !---------- calcul du second membre total
3611  DO j=1,6
3612  DO ni = 1, n_w
3613  u0(j_loc(ni),j) = u0(j_loc(ni),j) + ww(ni,l)*smb(j) + sum(dw_loc(:,ni)*grad(:,j,l))
3614  ENDDO
3615  ENDDO
3616 
3617  ENDDO
3618  ENDDO
3619 
3620  WRITE(*,*) ' CFL = ', cfl, ' R_eff for mode', mode, 6*surf/r_eff
3621  !IF (mode==0) THEN
3622  ! CALL plot_const_p1_label(mesh%jj, mesh%rr, visc_plot, 'tttt_0.plt')
3623  !END IF
3624 
3625  ed_st = normal_vit*coeff_ed_st
3626  DO ms = 1, mesh%mi
3627  dwni_loc = mesh%gauss%dwni(:,:,:,ms)
3628  jji_loc = mesh%jji(:,:,ms)
3629  h2 = -ed_st*h(ms)**2
3630 
3631  j_loc(1:n_ws) = mesh%jjsi(1:n_ws,ms)
3632  DO type = 1, 6
3633  DO cotei = 1, 2
3634  uloci(:,cotei) = vit(jji_loc(:,cotei),type)
3635  END DO
3636 
3637  u0loci = 0.d0
3638  DO ls = 1, mesh%gauss%l_Gs
3639  !===Compute radius of Gauss point
3640  !ray = SUM(mesh%rr(1,j_loc(1:n_ws))*wws(:,ls))
3641  ray = mesh%rr(1,j_loc(3))
3642 
3643  !--------Calcul du saut de la derivee normale de la vitesse
3644  dul = sum(dwni_loc(:,ls,:)*uloci)*mesh%gauss%rji(ls,ms)*h2*ray
3645  DO cotei = 1, 2
3646  DO ni = 1, mesh%gauss%n_w
3647  u0loci(ni, cotei) = u0loci(ni, cotei) + dwni_loc(ni,ls,cotei)*dul
3648  END DO
3649  END DO
3650  END DO
3651 
3652  DO cotei = 1, 2
3653  DO ni = 1, mesh%gauss%n_w
3654  u0(jji_loc(ni,cotei),type) = u0(jji_loc(ni,cotei),type) + u0loci(ni,cotei)
3655  END DO
3656  END DO
3657  END DO
3658 
3659  END DO
3660 
3661  END SUBROUTINE qs_ns_stab_2010
3662 
3663 
3664  SUBROUTINE qs_ns_stab_2008(mesh,mode,ff,vel_tot,V1m,vit,P,dt,u0,rotv_v)
3665  !=================================
3666  !Calcul du second membre pour le probleme de Stokes
3667  !sans terme non lineaire et source (ff)
3668 
3669  USE gauss_points
3670  USE chaine_caractere
3671  IMPLICIT NONE
3672  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotv_v !analytique (type,noeud)
3673  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: vel_tot !(noeud)
3674  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT):: u0 !u0(noeud, type)
3675  TYPE(mesh_type), TARGET :: mesh
3676  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1m, vit !V(noeud, type)
3677  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
3678  REAL(KIND=8), INTENT(IN) :: dt
3679  INTEGER, INTENT(IN) :: mode
3680  INTEGER :: m, l , i , ni , j, index, type, k
3681  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb, fv, mult
3682  REAL(KIND=8), DIMENSION(2,6,mesh%gauss%l_G) :: grad
3683  REAL(KIND=8), DIMENSION(6,mesh%gauss%l_G) :: vitloc
3684  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
3685  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3686 
3687  !fs pour le terme de forcage
3688  !ft pour le terme temporelle
3689  !fp pour le gradient de pression
3690 
3691  INTEGER, DIMENSION(:,:), POINTER :: jj
3692  INTEGER, POINTER :: me
3693  REAL(KIND=8) :: ray, div1, div2, vit1=0.005d0, vit2=.075d0, visc1, visc2, hloc, cfl, vloc
3694  REAL(KIND=8), SAVE :: R_eff, visc_eff
3695  LOGICAL, SAVE :: once = .true.
3696 
3697  CALL gauss(mesh)
3698  jj => mesh%jj
3699  me => mesh%me
3700 
3701  IF (once) THEN
3702  once =.false.
3703  END IF
3704 
3705  u0 = 0.d0
3706  r_eff = 0.d0
3707  cfl = 0
3708  index = 0
3709  DO m = 1, me
3710  j_loc = jj(:,m)
3711  hloc = sqrt(sum(rj(:,m)))
3712  vloc = min(maxval(vel_tot(j_loc)),1.d0)
3713  !vloc = maxval(vel_tot(j_loc))
3714  cfl = max(vloc*dt/hloc,cfl)
3715  visc1 = 0
3716  visc2 = 0
3717  DO l = 1, l_g
3718  dw_loc = dw(:,:,l,m)
3719 
3720  !===Compute radius of Gauss point
3721  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3722 
3723  !--------Calcul du gradient de la vitesse
3724  DO type = 1, 6
3725  DO k = 1 ,2
3726  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
3727  END DO
3728  vitloc(type,l) = sum(vit(j_loc,type)*ww(:,l))
3729  END DO
3730 
3731  !--------Calcul de la divergence
3732  div1 = grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray
3733  div2 = grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray
3734  visc1 = max(visc1,abs(div1))
3735  visc2 = max(visc2,abs(div2))
3736  END DO
3737  visc1 = max(visc1,visc2)
3738  !visc1 = 0.d0
3739  !visc_eff = hloc*MIN(inputs%coeff1,inputs%coeff2*hloc*visc1) - inputs%coeff1*hloc
3740  visc_eff = hloc*min(inputs%coeff1,inputs%coeff2*visc1)
3741  !======Semi-implicit version==========
3742  !visc1 = visc_eff - inputs%coeff1*hloc
3743  !======Semi-implicit version==========
3744  !======Totally explicit version=======
3745  visc1 = visc_eff
3746  !======Totally explicit version=======
3747  r_eff = r_eff + visc_eff*sum(rj(:,m))
3748  visc2 = visc1
3749 
3750  DO l = 1, l_g
3751  index = index +1
3752  dw_loc = dw(:,:,l,m)
3753 
3754  !===Compute radius of Gauss point
3755  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3756 
3757  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
3758  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
3759 
3760  fs(1) = sum(ff(j_loc,1) * ww(:,l))
3761  ft(1) = sum(v1m(j_loc,1) * ww(:,l))
3762  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3763  fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
3764  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
3765 
3766  fs(2) = sum(ff(j_loc,2) * ww(:,l)) !* rj(l,m)
3767  ft(2) = sum(v1m(j_loc,2) * ww(:,l)) !* rj(l,m)
3768  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:)) !* rj(l,m)
3769  fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
3770  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
3771 
3772  fs(3) = sum(ff(j_loc,3) * ww(:,l)) !* rj(l,m)
3773  ft(3) = sum(v1m(j_loc,3) * ww(:,l)) !* rj(l,m)
3774  fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode !* rj(l,m)
3775  fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
3776  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
3777 
3778  fs(4) = sum(ff(j_loc,4) * ww(:,l)) !* rj(l,m)
3779  ft(4) = sum(v1m(j_loc,4) * ww(:,l)) !* rj(l,m)
3780  fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode !* rj(l,m)
3781  fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
3782  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
3783 
3784  fs(5) = sum(ff(j_loc,5) * ww(:,l)) !* rj(l,m)
3785  ft(5) = sum(v1m(j_loc,5) * ww(:,l)) !* rj(l,m)
3786  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:)) !* rj(l,m)
3787  fv(5) = vitloc(5,l)*(mode/ray)**2
3788  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
3789 
3790  fs(6) = sum(ff(j_loc,6) * ww(:,l)) !* rj(l,m)
3791  ft(6) = sum(v1m(j_loc,6) * ww(:,l)) !* rj(l,m)
3792  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:)) !* rj(l,m)
3793  fv(6) = vitloc(6,l)*(mode/ray)**2
3794  !-------calcul du second membre pour le terme nonlineaire------------------------
3795  !-------------------------------------------------------------------------------
3796 
3797  mult(1) = visc1; mult(4) = visc1; mult(5) = visc1
3798  mult(2) = visc2; mult(3) = visc2; mult(6) = visc2
3799  fv = mult*fv
3800 
3801  smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*rj(l,m)
3802  DO type = 1, 6
3803  grad(:,type,l) = mult(type)*grad(:,type,l)*ray*rj(l,m)
3804  END DO
3805  !---------- calcul du second membre total
3806  DO j=1,6
3807  DO ni = 1, n_w
3808  u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j) + sum(dw_loc(:,ni)*grad(:,j,l))
3809  ENDDO
3810  ENDDO
3811 
3812  ENDDO
3813  ENDDO
3814 
3815  WRITE(*,*) ' CFL = ', cfl, ' R_eff ', 1/r_eff
3816  !IF (cfl > 10.d0) THEN
3817  ! WRITE(*,*) ' CFL too large '
3818  ! STOP
3819  !END IF
3820 
3821 
3822  END SUBROUTINE qs_ns_stab_2008
3823 
3824  SUBROUTINE qs_stokes_2006(mesh,mode,ff,V1m,P,dt,u0)
3825  !=================================
3826  !Calcul du second membre pour le probleme de Stokes
3827  !sans terme non lineaire et source (ff)
3828 
3829  USE gauss_points
3830 
3831  IMPLICIT NONE
3832  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff !analytique (type,noeud)
3833  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0 !u0(type,noeud)
3834  TYPE(mesh_type), TARGET :: mesh
3835  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: V1m !V(type,noeud)
3836  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: P
3837  REAL(KIND=8), INTENT(IN) :: dt
3838  INTEGER , INTENT(IN) :: mode
3839 
3840  INTEGER :: m, l , i , ni , j
3841  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, fnl, smb
3842  INTEGER :: N
3843  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
3844  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
3845  LOGICAL, SAVE :: once = .true.
3846  LOGICAL, SAVE :: once_sft = .true.
3847  REAL(KIND=8):: user_time
3848  !EXTERNAL :: user_time
3849  REAL(KIND=8) :: tps, dummy, tt
3850 
3851  !fs pour le terme de forcage
3852  !ft pour le terme temporelle
3853  !fp pour le gradient de pression
3854 
3855  INTEGER, DIMENSION(:,:), POINTER :: jj
3856  INTEGER, POINTER :: me
3857  REAL(KIND=8) :: ray
3858 
3859  CALL gauss(mesh)
3860  jj => mesh%jj
3861  me => mesh%me
3862 
3863  tps = 0.d0
3864  tt = 0.d0
3865 
3866  u0 = 0.d0
3867  fs = 0.d0
3868  ft = 0.d0
3869  fp = 0.d0
3870  fnl = 0.d0
3871 
3872  DO m = 1, me
3873  j_loc = jj(:,m)
3874  DO l = 1, l_g
3875  dw_loc = dw(:,:,l,m)
3876 
3877  !===Compute radius of Gauss point
3878  ray = sum(mesh%rr(1,j_loc)*ww(:,l))
3879 
3880  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
3881  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
3882 
3883  fs(1) = sum(ff(j_loc,1) * ww(:,l))
3884  ft(1) = 1/(2*dt) * sum(v1m(j_loc,1) * ww(:,l))
3885  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
3886 
3887  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
3888 
3889  fs(2) = sum(ff(j_loc,2) * ww(:,l)) !* rj(l,m)
3890  ft(2) = 1/(2*dt) * sum(v1m(j_loc,2) * ww(:,l)) !* rj(l,m)
3891  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:)) !* rj(l,m)
3892 
3893  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
3894 
3895  fs(3) = sum(ff(j_loc,3) * ww(:,l)) !* rj(l,m)
3896  ft(3) = 1/(2*dt) * sum(v1m(j_loc,3) * ww(:,l)) !* rj(l,m)
3897  fp(3) = -sum(p(j_loc,2)*ww(:,l))/ray*mode !* rj(l,m)
3898 
3899  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
3900 
3901  fs(4) = sum(ff(j_loc,4) * ww(:,l)) !* rj(l,m)
3902  ft(4) = 1/(2*dt) * sum(v1m(j_loc,4) * ww(:,l)) !* rj(l,m)
3903  fp(4) = sum(p(j_loc,1)*ww(:,l))/ray *mode !* rj(l,m)
3904 
3905  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
3906 
3907  fs(5) = sum(ff(j_loc,5) * ww(:,l)) !* rj(l,m)
3908  ft(5) = 1/(2*dt) * sum(v1m(j_loc,5) * ww(:,l)) !* rj(l,m)
3909  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:)) !* rj(l,m)
3910 
3911  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
3912 
3913  fs(6) = sum(ff(j_loc,6) * ww(:,l)) !* rj(l,m)
3914  ft(6) = 1/(2*dt) * sum(v1m(j_loc,6) * ww(:,l)) !* rj(l,m)
3915  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:)) !* rj(l,m)
3916 
3917  !-------calcul du second membre pour le terme nonlineaire------------------------
3918  !-------------------------------------------------------------------------------
3919 
3920  smb = (ft+fp+fs)*ray*rj(l,m)
3921  !write(95,*) m,l,ft(3)
3922  !---------- calcul du second membre total
3923  !smb = 1.d0
3924  DO j=1,6
3925  DO ni = 1, n_w
3926  u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l)*smb(j)
3927  ENDDO
3928  ENDDO
3929 
3930  ENDDO
3931  ENDDO
3932 
3933 
3934  END SUBROUTINE qs_stokes_2006
3935 
3936 
3937 
3938  SUBROUTINE qs_01_rot(mesh,mode,V,u0)
3939  !=================================
3940  !sans le terme de couplage
3941 
3942 
3943  USE gauss_points
3944 
3945  TYPE(mesh_type), TARGET :: mesh
3946  INTEGER, INTENT(IN) :: mode
3947  REAL(KIND=8), DIMENSION(mesh%np,6), INTENT(IN) :: V
3948  REAL(KIND=8), DIMENSION(:,:), INTENT(OUT) :: u0
3949 
3950  INTEGER :: m, l , i , ni , j
3951  REAL(KIND=8), DIMENSION(6) :: frot
3952  REAL(KIND=8):: user_time
3953  !EXTERNAL :: user_time
3954  REAL(KIND=8) :: tps, dummy, tt
3955  INTEGER, DIMENSION(:,:), POINTER :: jj
3956  INTEGER, POINTER :: me
3957  REAL(KIND=8) :: ray
3958 
3959  CALL gauss(mesh)
3960  jj => mesh%jj
3961  me => mesh%me
3962 
3963  u0 = 0.d0
3964  frot = 0.d0
3965 
3966  DO m = 1, me
3967 
3968  DO l = 1, l_g
3969 
3970  !===Compute radius of Gauss point
3971  ray = 0
3972  DO ni = 1, n_w; i = jj(ni,m)
3973  ray = ray + mesh%rr(1,i)*ww(ni,l)
3974  ENDDO
3975  !-----------------rotationnel ------------------------------
3976  !coeff sur les cosinus
3977  !ancien format
3978  !frot(1) = mode/ray*SUM(V(6,jj(:,m))*ww(:,l)) &
3979  ! -SUM(V(3,jj(:,m))*dw(2,:,l,m))
3980  !frot(3) = SUM(V(1,jj(:,m))*dw(2,:,l,m)) &
3981  ! -SUM(V(5,jj(:,m))*dw(1,:,l,m))
3982  !frot(5) = 1/ray*SUM(V(3,jj(:,m))*ww(:,l))+SUM(V(3,jj(:,m)) &
3983  ! *dw(1,:,l,m))-mode/ray*SUM(V(2,jj(:,m))*ww(:,l))
3984  ! !coeff sur les sinus
3985  !IF (mode /= 0) THEN
3986  ! frot(2) = -mode/ray*SUM(V(5,jj(:,m))*ww(:,l)) &
3987  ! -SUM(V(4,jj(:,m))*dw(2,:,l,m))
3988  ! frot(4) = SUM(V(2,jj(:,m))*dw(2,:,l,m)) &
3989  ! -SUM(V(6,jj(:,m))*dw(1,:,l,m))
3990  ! frot(6) = 1/ray*SUM(V(4,jj(:,m))*ww(:,l))+SUM(V(4,jj(:,m)) &
3991  ! *dw(1,:,l,m))+mode/ray*SUM(V(1,jj(:,m))*ww(:,l))
3992  !ENDIF
3993 
3994 
3995  frot(1) = mode/ray*sum(v(jj(:,m),6)*ww(:,l)) &
3996  -sum(v(jj(:,m),3)*dw(2,:,l,m))
3997  frot(3) = sum(v(jj(:,m),1)*dw(2,:,l,m)) &
3998  -sum(v(jj(:,m),5)*dw(1,:,l,m))
3999  frot(5) = 1/ray*sum(v(jj(:,m),3)*ww(:,l))+sum(v(jj(:,m),3) &
4000  *dw(1,:,l,m))-mode/ray*sum(v(jj(:,m),2)*ww(:,l))
4001  !coeff sur les sinus
4002  IF (mode /= 0) THEN
4003  frot(2) = -mode/ray*sum(v(jj(:,m),5)*ww(:,l)) &
4004  -sum(v(jj(:,m),4)*dw(2,:,l,m))
4005  frot(4) = sum(v(jj(:,m),2)*dw(2,:,l,m)) &
4006  -sum(v(jj(:,m),6)*dw(1,:,l,m))
4007  frot(6) = 1/ray*sum(v(jj(:,m),4)*ww(:,l))+sum(v(jj(:,m),4) &
4008  *dw(1,:,l,m))+mode/ray*sum(v(jj(:,m),1)*ww(:,l))
4009  ENDIF
4010  frot = frot * rj(l,m)
4011 
4012  !---------- calcul du second membre total
4013 
4014  DO j=1,6
4015  DO ni = 1, n_w
4016  u0(jj(ni,m),j) = u0(jj(ni,m),j) + ww(ni,l) * ray * frot(j)
4017  ENDDO
4018  ENDDO
4019 
4020  ENDDO
4021  ENDDO
4022 
4023  END SUBROUTINE qs_01_rot
4024 
4025 
4026  SUBROUTINE average(mesh,visc)
4027  USE gauss_points
4028  USE sub_plot
4029  IMPLICIT NONE
4030  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: visc
4031  TYPE(mesh_type) :: mesh
4032  REAL(KIND=8), DIMENSION(mesh%me) :: vint
4033  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: surf, sloc
4034  LOGICAL, SAVE :: once=.true.
4035  INTEGER :: m, n, mp, it
4036 
4037  IF (once) THEN
4038  once =.false.
4039  ALLOCATE(surf(mesh%me),sloc(mesh%me))
4040  DO m = 1, mesh%me
4041  surf(m) = sum(mesh%gauss%rj(:,m))
4042  END DO
4043  DO m = 1, mesh%me
4044  sloc(m) = 0.d0
4045  DO n = 1, 3
4046  IF (mesh%neigh(n,m)==0) cycle
4047  mp = mesh%neigh(n,m)
4048  sloc(m) = sloc(m) + surf(m)+surf(mp)
4049  END DO
4050  END DO
4051 
4052  END IF
4053  vint = 0.d0
4054  DO m = 1, mesh%me
4055  DO n = 1, 3
4056  IF (mesh%neigh(n,m)==0) cycle
4057  mp = mesh%neigh(n,m)
4058  vint(m) = vint(m)+ visc(m)*surf(m) + visc(mp)*surf(mp)
4059  END DO
4060  vint(m) = vint(m)/sloc(m)
4061  END DO
4062  visc = vint
4063 
4064  END SUBROUTINE average
4065 
4066  SUBROUTINE clement_c(mesh,phi,phi_s)
4067  !USE GAUSS_POINTS
4068  USE def_type_mesh
4069  IMPLICIT NONE
4070  TYPE(mesh_type) :: mesh
4071  REAL(KIND=8), DIMENSION(:) :: phi, phi_s
4072 
4073  INTEGER :: m, l, k, i
4074  REAL(KIND=8), DIMENSION(mesh%np) :: a0, b1, b2
4075  REAL(KIND=8), DIMENSION(:), ALLOCATABLE, SAVE :: s, a11, a22, a12, oneoverdet
4076  REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE, SAVE :: rg
4077  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
4078  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: phi_loc
4079  REAL(KIND=8), DIMENSION(mesh%gauss%k_d) :: r_loc
4080  REAL(KIND=8) :: x1l, x2l, phil, a1, a2, det
4081  LOGICAL :: const
4082  LOGICAL, SAVE :: once=.true.
4083 
4084  IF (mesh%gauss%n_w/=3) THEN
4085  WRITE(*,*) ' BUG CLEMENT: Only P1 programmed yet'
4086  phi_s= phi
4087  RETURN
4088  END IF
4089  const=.false.
4090  IF (SIZE(phi)==mesh%me) const=.true.
4091 
4092  IF (once) THEN
4093  once = .false.
4094  ALLOCATE(rg(mesh%gauss%k_d,mesh%np), s(mesh%np))
4095  s = 0.d0
4096  rg = 0.d0
4097  DO m = 1, mesh%me
4098  j_loc = mesh%jj(:,m)
4099  s(j_loc(:)) = s(j_loc(:)) + sum(mesh%gauss%rj(:,m))
4100  DO k = 1, mesh%gauss%k_d
4101  r_loc(k) = 0.d0
4102  DO l = 1, mesh%gauss%l_G
4103  r_loc(k) = r_loc(k) + sum(mesh%rr(k,j_loc(:))*mesh%gauss%ww(:,l))*mesh%gauss%rj(l,m)
4104  END DO
4105  rg(k,j_loc(:)) = rg(k,j_loc(:)) + r_loc(k)
4106  END DO
4107  END DO
4108  DO i = 1, mesh%np
4109  rg(:,i) = rg(:,i)/s(i)
4110  END DO
4111  ALLOCATE(a11(mesh%np), a22(mesh%np), a12(mesh%np), oneoverdet(mesh%np))
4112  a11 = 0.d0
4113  a12 = 0.d0
4114  a22 = 0.d0
4115  DO m = 1, mesh%me
4116  j_loc = mesh%jj(:,m)
4117  DO l = 1, mesh%gauss%l_G
4118  x1l = sum(mesh%rr(1,j_loc(:))*mesh%gauss%ww(:,l))
4119  x2l = sum(mesh%rr(2,j_loc(:))*mesh%gauss%ww(:,l))
4120  a11(j_loc(:)) = a11(j_loc(:)) + (x1l-rg(1,j_loc(:)))**2*mesh%gauss%rj(l,m)
4121  a22(j_loc(:)) = a22(j_loc(:)) + (x2l-rg(2,j_loc(:)))**2*mesh%gauss%rj(l,m)
4122  a12(j_loc(:)) = a12(j_loc(:)) + (x1l-rg(1,j_loc(:)))*(x2l-rg(2,j_loc(:)))*mesh%gauss%rj(l,m)
4123  END DO
4124  END DO
4125  oneoverdet = 1.d0/(a11*a22-a12*a12)
4126  END IF
4127 
4128 
4129  a0 = 0.d0
4130  b1 = 0.d0
4131  b2 = 0.d0
4132  DO m = 1, mesh%me
4133  j_loc = mesh%jj(:,m)
4134  !IF (const) THEN
4135  phi_loc = phi(m)
4136  !ELSE
4137  ! phi_loc = phi(j_loc)
4138  !END IF
4139  DO l = 1, mesh%gauss%l_G
4140  x1l = sum(mesh%rr(1,j_loc(:))*mesh%gauss%ww(:,l))
4141  x2l = sum(mesh%rr(2,j_loc(:))*mesh%gauss%ww(:,l))
4142  phil = sum(phi_loc*mesh%gauss%ww(:,l))
4143  b1(j_loc(:)) = b1(j_loc(:)) + (x1l-rg(1,j_loc(:)))*phil*mesh%gauss%rj(l,m)
4144  b2(j_loc(:)) = b2(j_loc(:)) + (x2l-rg(2,j_loc(:)))*phil*mesh%gauss%rj(l,m)
4145  a0(j_loc) = a0(j_loc) + phil*mesh%gauss%rj(l,m)
4146  END DO
4147  END DO
4148 
4149  DO i = 1, mesh%np
4150  a0(i)=a0(i)/s(i)
4151  !det = a11(i)*a22(i)-a12(i)*a12(i)
4152  a1 =(b1(i)*a22(i)-b2(i)*a12(i))*oneoverdet(i)
4153  a2 =(b2(i)*a11(i)-b1(i)*a12(i))*oneoverdet(i)
4154  phi_s(i) = a0(i)+a1*(mesh%rr(1,i)-rg(1,i))+a2*(mesh%rr(2,i)-rg(2,i))
4155  END DO
4156 
4157  END SUBROUTINE clement_c
4158 
4159 
4160  SUBROUTINE inject_clement(jj_c, jj_f, pp_c, pp_f)
4161 
4162  IMPLICIT NONE
4163 
4164  INTEGER, DIMENSION(:,:), INTENT(IN) :: jj_c, jj_f
4165  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: pp_c
4166  REAL(KIND=8), DIMENSION(:), INTENT(OUT) :: pp_f
4167  REAL(KIND=8) :: half = 0.5
4168  INTEGER:: m
4169 
4170  DO m = 1, SIZE(jj_f,2)
4171  pp_f(jj_f(1:3,m)) = pp_c(jj_c(:,m))
4172  pp_f(jj_f(4,m)) = (pp_c(jj_c(2,m)) + pp_c(jj_c(3,m)))*half
4173  pp_f(jj_f(5,m)) = (pp_c(jj_c(3,m)) + pp_c(jj_c(1,m)))*half
4174  pp_f(jj_f(6,m)) = (pp_c(jj_c(1,m)) + pp_c(jj_c(2,m)))*half
4175  END DO
4176 
4177  END SUBROUTINE inject_clement
4178 END MODULE fem_s_axi
fem_s_axi::average
subroutine average(mesh, visc)
Definition: fem_s_direct_axi.f90:4027
fem_s_axi::clement_c
subroutine clement_c(mesh, phi, phi_s)
Definition: fem_s_direct_axi.f90:4067
fem_s_axi::qs_00_adv_diff_vect_3d_init_ssr
subroutine qs_00_adv_diff_vect_3d_init_ssr(mesh, alpha, mode, gg, ff, V, dt, u0)
Definition: fem_s_direct_axi.f90:1711
fem_s_axi::qs_00_inst_vect_3d_init
subroutine qs_00_inst_vect_3d_init(mesh, alpha, mode, ff, V2, dt, u0)
Definition: fem_s_direct_axi.f90:1512
sub_plot
Definition: sub_plot.f90:4
fem_s_axi::inject_clement
subroutine inject_clement(jj_c, jj_f, pp_c, pp_f)
Definition: fem_s_direct_axi.f90:4161
gauss_points
Definition: associate_gauss.f90:4
fem_s_axi::qs_00_inst_init_3d_ssr
subroutine qs_00_inst_init_3d_ssr(mesh, ff, T1, dt, u0)
Definition: fem_s_direct_axi.f90:323
fem_s_axi::qs_00_inst_vect_3d_init_ssr
subroutine qs_00_inst_vect_3d_init_ssr(mesh, alpha, mode, ff, V, dt, u0)
Definition: fem_s_direct_axi.f90:1609
fem_s_axi::qs_00_rot
subroutine qs_00_rot(mesh, mode, V, Rot)
Definition: fem_s_direct_axi.f90:746
gauss_points::l_g
integer, public l_g
Definition: associate_gauss.f90:9
gauss_points::n_ws
integer, public n_ws
Definition: associate_gauss.f90:9
fem_s_axi::qs_ns_stab_new
subroutine qs_ns_stab_new(mesh, mode, ff, vel_tot, V1m, vit, P, dudt, phalf, nlhalf, dt, u0, rotv_v)
Definition: fem_s_direct_axi.f90:3109
fem_s_axi::qs_01_grad
subroutine qs_01_grad(mesh, mode, pp, u0)
Definition: fem_s_direct_axi.f90:1184
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:74
gauss_points::rj
real(kind=8), dimension(:,:), pointer rj
Definition: associate_gauss.f90:14
fem_s_axi::qs_01_grad_gl
subroutine qs_01_grad_gl(mesh, mod_max, pp, u0)
Definition: fem_s_direct_axi.f90:1240
fem_s_axi::qs_ns_2006
subroutine qs_ns_2006(mesh, mode, ff, V1m, P, dt, u0, rotv_v)
Definition: fem_s_direct_axi.f90:3009
fem_s_axi::qs_00_ns_inline
subroutine qs_00_ns_inline(mesh, mode, m_max, ff, V1, V2, P, dt, u0, meth, tps_sft)
Definition: fem_s_direct_axi.f90:2546
input_data::inputs
type(my_data), public inputs
Definition: read_my_data.f90:237
fem_s_axi::qs_01_div_hybrid_old
subroutine qs_01_div_hybrid_old(uu_mesh, pp_mesh, mode, gg, u0_c)
Definition: fem_s_direct_axi.f90:884
fem_s_axi::moy
subroutine moy(mesh, p, RESULT)
Definition: fem_s_direct_axi.f90:35
fem_s_axi::qs_00_adv_diff_vect_3d_ssr
subroutine qs_00_adv_diff_vect_3d_ssr(mesh, alpha, mode, gg, ff, V1, V2, dt, u0)
Definition: fem_s_direct_axi.f90:1916
fem_s_axi::qs_00_lap
subroutine qs_00_lap(mesh, alpha, mode, ff, V2, V1, dt, u0)
Definition: fem_s_direct_axi.f90:615
fem_s_axi::qs_00_lap_init
subroutine qs_00_lap_init(mesh, alpha, mode, ff, V, dt, u0)
Definition: fem_s_direct_axi.f90:681
fem_s_axi::qs_00_adv_diff_init_3d
subroutine qs_00_adv_diff_init_3d(eps, mode, mesh, ff, T1, dt, gg, u0)
Definition: fem_s_direct_axi.f90:500
fem_s_axi::qs_00_adv_diff_init_3d_ssr
subroutine qs_00_adv_diff_init_3d_ssr(eps, mode, mesh, ff, T1, dt, gg, u0)
Definition: fem_s_direct_axi.f90:553
gauss_points::n_w
integer, public n_w
Definition: associate_gauss.f90:9
fem_s_axi::qs_00_inst_init_3d
subroutine qs_00_inst_init_3d(mesh, ff, T1, dt, u0)
Definition: fem_s_direct_axi.f90:274
fem_s_axi::qs_00_adv_diff_3d
subroutine qs_00_adv_diff_3d(eps, mode, mesh, ff, T1, T2, gg, dt, u0)
Definition: fem_s_direct_axi.f90:373
fem_s_axi::dirichlet
subroutine dirichlet(jjs, us, ff)
Definition: fem_s_direct_axi.f90:18
fem_s_axi::qs_00_div
subroutine qs_00_div(mesh, mode, V, u0)
Definition: fem_s_direct_axi.f90:804
fem_s_axi::qs_00_inst_vect_3d_ssr
subroutine qs_00_inst_vect_3d_ssr(mesh, alpha, mode, ff, V1, V2, dt, u0)
Definition: fem_s_direct_axi.f90:1412
fem_s_axi::qs_stokes
subroutine qs_stokes(mesh, mode, ff, V1m, V2m, P, dt, u0)
Definition: fem_s_direct_axi.f90:2807
fem_s_axi::qs_00_inst_3d_ssr
subroutine qs_00_inst_3d_ssr(mesh, ff, T1, T2, dt, u0)
Definition: fem_s_direct_axi.f90:174
fem_s_axi::qs_navier_stokes_2006
subroutine qs_navier_stokes_2006(mesh, mode, ff, V1m, P, dt, u0, rotv_v)
Definition: fem_s_direct_axi.f90:2909
fem_s_axi::qs_01_rot
subroutine qs_01_rot(mesh, mode, V, u0)
Definition: fem_s_direct_axi.f90:3939
fem_s_axi::qs_ns_stab_2008
subroutine qs_ns_stab_2008(mesh, mode, ff, vel_tot, V1m, vit, P, dt, u0, rotv_v)
Definition: fem_s_direct_axi.f90:3665
fem_s_axi::qs_00_adv_diff_vect_3d_init
subroutine qs_00_adv_diff_vect_3d_init(mesh, alpha, mode, gg, ff, V, dt, u0)
Definition: fem_s_direct_axi.f90:1814
fem_s_axi::qs_00_stokes_3d_init
subroutine qs_00_stokes_3d_init(mesh, alpha, mode, ff, P, V, dt, u0)
Definition: fem_s_direct_axi.f90:2125
fem_s_axi::qs_01_div_hybrid_2006
subroutine qs_01_div_hybrid_2006(uu_mesh, pp_mesh, mode, gg, u0_c)
Definition: fem_s_direct_axi.f90:1101
fem_s_axi::qs_stokes_2006
subroutine qs_stokes_2006(mesh, mode, ff, V1m, P, dt, u0)
Definition: fem_s_direct_axi.f90:3825
fem_s_axi
Definition: fem_s_direct_axi.f90:6
input_data
Definition: read_my_data.f90:233
fem_s_axi::qs_00_inst_vect_3d
subroutine qs_00_inst_vect_3d(mesh, alpha, mode, ff, V1, V2, dt, u0)
Definition: fem_s_direct_axi.f90:1309
fem_s_axi::qs_01_div_hybrid
subroutine qs_01_div_hybrid(uu_mesh, pp_mesh, mode, gg, u0_c)
Definition: fem_s_direct_axi.f90:1012
gauss_points::gauss
subroutine gauss(mesh)
Definition: associate_gauss.f90:24
chaine_caractere
Definition: chaine_caractere.f90:4
fem_s_axi::qs_00_adv_diff_3d_ssr
subroutine qs_00_adv_diff_3d_ssr(eps, mode, mesh, ff, T1, T2, gg, dt, u0)
Definition: fem_s_direct_axi.f90:431
def_type_mesh
Definition: def_type_mesh.f90:4
fem_s_axi::qs_00_ssr
subroutine qs_00_ssr(mesh, ff, u0)
Definition: fem_s_direct_axi.f90:82
fem_s_axi::qs_00_stokes_3d_new
subroutine qs_00_stokes_3d_new(mesh, alpha, mode, ff, V1, V2, P, dt, u0)
Definition: fem_s_direct_axi.f90:2446
gauss_points::dw
real(kind=8), dimension(:,:,:,:), pointer dw
Definition: associate_gauss.f90:12
fem_s_axi::qs_00
subroutine qs_00(mesh, ff, u0)
Definition: fem_s_direct_axi.f90:128
fem_s_axi::qs_00_stokes_3d
subroutine qs_00_stokes_3d(mesh, alpha, mode, ff, V1, V2, P, dt, u0)
Definition: fem_s_direct_axi.f90:2248
gauss_points::ww
real(kind=8), dimension(:,:), pointer ww
Definition: associate_gauss.f90:10
fem_s_axi::qs_00_stokes_3d_new_ssr
subroutine qs_00_stokes_3d_new_ssr(mesh, alpha, mode, ff, V1, V2, P, dt, u0)
Definition: fem_s_direct_axi.f90:2350
fem_s_axi::qs_00_adv_diff_vect_3d
subroutine qs_00_adv_diff_vect_3d(mesh, alpha, mode, gg, ff, V1, V2, dt, u0)
Definition: fem_s_direct_axi.f90:2020
fem_s_axi::qs_00_inst_3d
subroutine qs_00_inst_3d(mesh, ff, T1, T2, dt, u0)
Definition: fem_s_direct_axi.f90:223
fem_s_axi::qs_ns_stab_2010
subroutine qs_ns_stab_2010(mesh, pp_mesh, mode, ff, vel_tot, V1m, vit, P, dudt, phalf, nlhalf, dt, u0, rotv_v)
Definition: fem_s_direct_axi.f90:3390