SFEMaNS  version 4.1 (work in progress)
Reference documentation for SFEMaNS
 All Classes Files Functions Variables Groups Pages
fem_rhs_axi.f90
Go to the documentation of this file.
1 MODULE fem_rhs_axi
2 
3 CONTAINS
4 
5  SUBROUTINE qs_ns_stab_new(mesh,vv_1_LA,vv_2_LA,mode,ff,vel_tot,V1m,vit,P,dudt,phalf,nlhalf,dt,&
6  vb_2_14,vb_2_23,vb_1_5,vb_1_6,rotv_v, vel_tot_max)
7  !=================================
8  !RHS for Navier-Stokes
9  USE def_type_mesh
10  USE chaine_caractere
11  USE my_util
12  USE input_data
13  !USE sub_plot
14  IMPLICIT NONE
15  TYPE(petsc_csr_la) :: vv_1_la,vv_2_la
16  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotv_v, dudt, nlhalf
17  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: vel_tot !(noeud)
18  TYPE(mesh_type), TARGET :: mesh
19  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: v1m, vit !V(noeud, type)
20  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: p, phalf
21  REAL(KIND=8), INTENT(IN) :: dt
22  INTEGER, INTENT(IN) :: mode
23 !TEST LES LC October 28, 2014
24  REAL(KIND=8) :: vel_tot_max
25 !TEST LES LC October 28, 2014
26  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb, fv, mult
27  REAL(KIND=8), DIMENSION(2,6,mesh%gauss%l_G) :: grad
28  REAL(KIND=8), DIMENSION(6,mesh%gauss%l_G) :: vitloc
29  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
30  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
31  REAL(KIND=8), DIMENSION(mesh%dom_me) :: visc_plot
32  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_w) :: v14_loc, v23_loc
33  INTEGER, DIMENSION(2*mesh%gauss%n_w) :: idxm_2
34  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v5_loc, v6_loc
35  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm_1
36  REAL(KIND=8), DIMENSION(6) :: visc1, visc2
37  REAL(KIND=8) :: ray, div1, div2, cfl, vloc, normal_vit
38  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: h
39  REAL(KIND=8), SAVE :: coeff_ed_st, r_eff, &
40  visc_eff, surf, nu_loc, coeff_visc_ordre_un
41  LOGICAL, SAVE :: once = .true.
42  REAL(KIND=8), DIMENSION(6,mesh%dom_me) :: viscosity
43  REAL(KIND=8), DIMENSION(6) :: norm_vit
44  REAL(KIND=8) :: type_fe
45  INTEGER :: m, l , i , ni , index, index2, type, k
46  INTEGER :: ms, nw, ix, ki, iglob
47 !!$ WARNING FL Variables removed (used for edge_stab)
48 !!$ INTEGER :: cotei, ls
49 !!$ REAL(KIND=8) :: dul, ed_st, h2
50 !!$ INTEGER, DIMENSION(mesh%gauss%n_w,2) :: jji_loc
51 !!$ REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: u0loci, uloci
52 !!$ REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%l_Gs, 2) :: dwni_loc
53 !!$ WARNING FL Variables removed (used for edge_stab)
54 #include "petsc/finclude/petsc.h"
55  vec :: vb_2_14,vb_2_23,vb_1_5,vb_1_6
56  petscerrorcode :: ierr
57 
58  CALL vecset(vb_2_14, 0.d0, ierr)
59  CALL vecset(vb_2_23, 0.d0, ierr)
60  CALL vecset(vb_1_5, 0.d0, ierr)
61  CALL vecset(vb_1_6, 0.d0, ierr)
62 
63  IF (once) THEN
64  once =.false.
65 
66  IF (.NOT.inputs%LES) THEN
67  inputs%LES_coeff1=0.d0
68  inputs%LES_coeff2=0.d0
69  inputs%LES_coeff3=0.d0
70  inputs%LES_coeff4=0.d0
71  END IF
72 
73  surf = 0.d0
74  DO m = 1, mesh%dom_me
75  DO l = 1, mesh%gauss%l_G
76  ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
77  surf = surf + mesh%gauss%rj(l,m)*ray
78  END DO
79  END DO
80 
81  IF (mesh%gauss%n_w==3) THEN
82  type_fe = 1
83  ELSE
84  type_fe = 2
85  END IF
86  coeff_ed_st = inputs%LES_coeff3*0.02d0/type_fe
87  coeff_visc_ordre_un = inputs%LES_coeff4
88  IF (mesh%edge_stab) THEN
89  ALLOCATE(h(mesh%mi))
90  DO ms = 1, mesh%mi
91  h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
92  END DO
93  END IF
94  END IF
95 
96  !ATTENTION: JLG Jan 25 2010
97  !ATTENTION: inputs%LES_coeff1 is assumed to be of the order of the convective velocity
98  !that simplifies the semi-implicit treatment of the LES viscosity
99 !!$ normal_vit = MAXVAL(vel_tot)
100 !TEST LES LC October 28, 2014
101  normal_vit = vel_tot_max
102 !TEST LES LC October 28, 2014
103  DO TYPE = 1, 6
104  norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
105  END DO
106  !ATTENTION: JLG Jan 25 2010
107  r_eff = 0.d0
108  cfl = 0
109  index = 0
110  index2 = 0
111  IF (inputs%LES) THEN
112  DO m = 1, mesh%dom_me
113  j_loc = mesh%jj(:,m)
114  vloc = maxval(vel_tot(j_loc))
115  cfl = max(vloc*dt/mesh%hloc(m),cfl)
116  visc1 = 0
117  visc2 = 0
118  DO l = 1, mesh%gauss%l_G
119  index2 = index2 +1
120  dw_loc = mesh%gauss%dw(:,:,l,m)
121 
122  !--------radius of gauss point
123  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
124  !fs Source term
125  !ft Time derivative
126  !fp Pressure gradient
127  !--------calcul de la premiere composante
128  ft(1) = sum(dudt(j_loc,1) *mesh%gauss%ww(:,l))
129  fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
130  !--------calcul de la seconde composante
131  ft(2) = sum(dudt(j_loc,2) * mesh%gauss%ww(:,l))
132  fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
133  !--------calcul de la troisieme composante
134  ft(3) = sum(dudt(j_loc,3) *mesh%gauss%ww(:,l))
135  fp(3) = sum(phalf(j_loc,2)*mesh%gauss%ww(:,l))*mode/ray
136  !--------calcul de la quatrieme composante
137  ft(4) = sum(dudt(j_loc,4) *mesh%gauss%ww(:,l))
138  fp(4) = -sum(phalf(j_loc,1)*mesh%gauss%ww(:,l))*mode/ray
139  !--------calcul de la cinquieme composante
140  ft(5) = sum(dudt(j_loc,5) *mesh%gauss%ww(:,l))
141  fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
142  !--------calcul de la sixieme composante
143  ft(6) = sum(dudt(j_loc,6) *mesh%gauss%ww(:,l))
144  fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
145  !-------calcul du second membre pour le terme nonlineaire------------------------
146 
147  visc1 = max(visc1,abs(ft+fp+nlhalf(index2,:)))
148 
149  !--------Calcul du gradient de la vitesse
150  DO TYPE = 1, 6
151  DO k = 1 ,2
152  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
153  END DO
154  vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
155  END DO
156 
157  !--------Calcul de la divergence
158  div1 = abs(grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray)
159  div2 = abs(grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray)
160  visc2(1) = max(visc2(1),div1)
161  visc2(2) = max(visc2(2),div2)
162  END DO
163  visc2(4) = visc2(1); visc2(5) = visc2(1)
164  visc2(3) = visc2(2); visc2(6) = visc2(2)
165 
166  nu_loc = 0.d0
167  DO TYPE = 1, 6
168  !visc1(type) = MAX(visc1(type)/normal_vit,visc2(type))
169  !visc1(type) = MAX(visc1(type),visc2(type)*normal_vit)/MAXVAL(ABS(vitloc(type,:)))
170  visc1(type) = max(visc1(type),2*visc2(type)*normal_vit)/norm_vit(type)
171  visc_eff = mesh%hloc(m)*min(coeff_visc_ordre_un*normal_vit,inputs%LES_coeff2*mesh%hloc(m)*visc1(type))
172  nu_loc = nu_loc + visc_eff
173  !======Semi-implicit version==========
174  viscosity(type,m) = inputs%LES_coeff1*mesh%hloc(m) - visc_eff
175  !======Semi-implicit version==========
176  END DO
177  r_eff = r_eff + (nu_loc + dt**2*mesh%hloc(m))*mesh%hloc(m)**2*ray
178  visc_plot(m) = (nu_loc/6)/(coeff_visc_ordre_un*mesh%hloc(m)*normal_vit)
179  END DO
180  ELSE
181  cfl = 0.d0
182  viscosity = 0.d0
183  DO m = 1, mesh%dom_me
184  j_loc = mesh%jj(:,m)
185  vloc = maxval(vel_tot(j_loc))
186  cfl = max(vloc*dt/mesh%hloc(m),cfl)
187  END DO
188  END IF
189 
190 
191  !DO type = 1, 6
192  ! CALL average(mesh,viscosity(type,:))
193  !END DO
194 
195  nw = mesh%gauss%n_w
196  DO m = 1, mesh%dom_me
197  mult(:)= viscosity(:,m)
198  j_loc = mesh%jj(:,m)
199 
200  DO ni = 1, nw
201  i = j_loc(ni)
202  iglob = vv_1_la%loc_to_glob(1,i)
203  idxm_1(ni) = iglob-1
204  END DO
205 
206  DO ki = 1, 2
207  DO ni = 1, nw
208  i = j_loc(ni)
209  iglob = vv_2_la%loc_to_glob(ki,i)
210  ix = (ki-1)*nw+ni
211  idxm_2(ix) = iglob-1
212  END DO
213  END DO
214 
215  v14_loc = 0.d0
216  v23_loc = 0.d0
217  v5_loc = 0.d0
218  v6_loc = 0.d0
219  DO l = 1, mesh%gauss%l_G
220  index = index +1
221  dw_loc = mesh%gauss%dw(:,:,l,m)
222  DO TYPE = 1, 6
223  DO k = 1 ,2
224  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
225  END DO
226  vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
227  END DO
228 
229  !===Compute radius of Gauss point
230  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
231 
232  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
233  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
234  fs(1) = sum(ff(j_loc,1) * mesh%gauss%ww(:,l))
235  ft(1) = sum(v1m(j_loc,1) * mesh%gauss%ww(:,l))
236  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
237  fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
238  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
239  fs(2) = sum(ff(j_loc,2) * mesh%gauss%ww(:,l))
240  ft(2) = sum(v1m(j_loc,2) * mesh%gauss%ww(:,l))
241  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
242  fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
243  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
244  fs(3) = sum(ff(j_loc,3) * mesh%gauss%ww(:,l))
245  ft(3) = sum(v1m(j_loc,3) * mesh%gauss%ww(:,l))
246  fp(3) = -sum(p(j_loc,2)*mesh%gauss%ww(:,l))/ray*mode
247  fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
248  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
249  fs(4) = sum(ff(j_loc,4) * mesh%gauss%ww(:,l))
250  ft(4) = sum(v1m(j_loc,4) * mesh%gauss%ww(:,l))
251  fp(4) = sum(p(j_loc,1)*mesh%gauss%ww(:,l))/ray *mode
252  fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
253  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
254  fs(5) = sum(ff(j_loc,5) * mesh%gauss%ww(:,l))
255  ft(5) = sum(v1m(j_loc,5) * mesh%gauss%ww(:,l))
256  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
257  fv(5) = vitloc(5,l)*(mode/ray)**2
258  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
259  fs(6) = sum(ff(j_loc,6) * mesh%gauss%ww(:,l))
260  ft(6) = sum(v1m(j_loc,6) * mesh%gauss%ww(:,l))
261  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
262  fv(6) = vitloc(6,l)*(mode/ray)**2
263  !-------calcul du second membre pour le terme nonlineaire------------------------
264  !-------------------------------------------------------------------------------
265 
266  fv = mult*fv
267 
268  smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*mesh%gauss%rj(l,m)
269  DO TYPE = 1, 6
270  grad(:,type,l) = mult(type)*grad(:,type,l)*ray*mesh%gauss%rj(l,m)
271  END DO
272 
273 !!$ DO j=1,6
274 !!$ DO ni = 1, mesh%gauss%n_w
275 !!$ u0(j_loc(ni),j) = u0(j_loc(ni),j) + mesh%gauss%ww(ni,l)*smb(j) + SUM(dw_loc(:,ni)*grad(:,j,l))
276 !!$ ENDDO
277 !!$ ENDDO
278 
279  ki = 1
280  DO ni = 1, nw
281  ix = (ki-1)*nw+ni
282  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(1) + sum(dw_loc(:,ni)*grad(:,1,l))
283  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(2) + sum(dw_loc(:,ni)*grad(:,2,l))
284  v5_loc(ix) = v5_loc(ix) + mesh%gauss%ww(ni,l)*smb(5) + sum(dw_loc(:,ni)*grad(:,5,l))
285  v6_loc(ix) = v6_loc(ix) + mesh%gauss%ww(ni,l)*smb(6) + sum(dw_loc(:,ni)*grad(:,6,l))
286  END DO
287 
288  ki = 2
289  DO ni = 1, nw
290  ix = (ki-1)*nw+ni
291  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(4) + sum(dw_loc(:,ni)*grad(:,4,l))
292  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(3) + sum(dw_loc(:,ni)*grad(:,3,l))
293  END DO
294 
295  ENDDO
296 
297  CALL vecsetvalues(vb_2_14, 2*nw, idxm_2, v14_loc, add_values, ierr)
298  CALL vecsetvalues(vb_2_23, 2*nw, idxm_2, v23_loc, add_values, ierr)
299  CALL vecsetvalues(vb_1_5, nw, idxm_1, v5_loc, add_values, ierr)
300  CALL vecsetvalues(vb_1_6, nw, idxm_1, v6_loc, add_values, ierr)
301  ENDDO
302 
303  !WRITE(*,*) ' CFL = ', cfl, ' R_eff for mode', mode, normal_vit*6*surf/R_eff
304  !IF (mode==0) THEN
305  !CALL plot_const_p1_label(mesh%jj, mesh%rr, visc_plot, 'tttt_0.plt')
306  !END IF
307 
308  IF (inputs%LES) THEN
309  IF (mesh%edge_stab) THEN
310  CALL error_petsc('BUG in qs_ns_stab_new: terms with edge_stab not yet assembled')
311 !!$ ed_st = normal_vit*coeff_ed_st
312 !!$ DO ms = 1, mesh%mi
313 !!$ dwni_loc = mesh%gauss%dwni(:,:,:,ms)
314 !!$ jji_loc = mesh%jji(:,:,ms)
315 !!$ h2 = -ed_st*h(ms)**2
316 !!$
317 !!$ j_loc(1:mesh%gauss%n_ws) = mesh%jjsi(1:mesh%gauss%n_ws,ms)
318 !!$ DO TYPE = 1, 6
319 !!$ DO cotei = 1, 2
320 !!$ uloci(:,cotei) = vit(jji_loc(:,cotei),TYPE)
321 !!$ END DO
322 !!$
323 !!$ u0loci = 0.d0
324 !!$ DO ls = 1, mesh%gauss%l_Gs
325 !!$ !===Compute radius of Gauss point at face center
326 !!$ ray = mesh%rr(1,j_loc(3))
327 !!$
328 !!$ !--------Calcul du saut de la derivee normale de la vitesse
329 !!$ dul = SUM(dwni_loc(:,ls,:)*uloci)*mesh%gauss%rji(ls,ms)*h2*ray
330 !!$ DO cotei = 1, 2
331 !!$ DO ni = 1, mesh%gauss%n_w
332 !!$ u0loci(ni, cotei) = u0loci(ni, cotei) + dwni_loc(ni,ls,cotei)*dul
333 !!$ END DO
334 !!$ END DO
335 !!$ END DO
336 !!$
337 !!$ DO cotei = 1, 2
338 !!$ DO ni = 1, mesh%gauss%n_w
339 !!$ u0(jji_loc(ni,cotei),TYPE) = u0(jji_loc(ni,cotei),TYPE) + u0loci(ni,cotei)
340 !!$ END DO
341 !!$ END DO
342 !!$ END DO
343 
344  END IF
345  END IF
346 
347  CALL vecassemblybegin(vb_2_14,ierr)
348  CALL vecassemblyend(vb_2_14,ierr)
349  CALL vecassemblybegin(vb_2_23,ierr)
350  CALL vecassemblyend(vb_2_23,ierr)
351  CALL vecassemblybegin(vb_1_5,ierr)
352  CALL vecassemblyend(vb_1_5,ierr)
353  CALL vecassemblybegin(vb_1_6,ierr)
354  CALL vecassemblyend(vb_1_6,ierr)
355 
356  END SUBROUTINE qs_ns_stab_new
357 
358  SUBROUTINE qs_ns_stab_3x3(mesh,vv_3_LA,mode,ff,vel_tot,V1m,vit,P,dudt,phalf,nlhalf,dt,&
359  vb_145, vb_236,rotv_v, vel_tot_max)
360  !=================================
361  !RHS for Navier-Stokes
362  USE def_type_mesh
363  USE chaine_caractere
364  USE my_util
365  USE input_data
366  USE rhs_para_assembling
367  !USE sub_plot
368  IMPLICIT NONE
369  TYPE(petsc_csr_la) :: vv_3_la
370  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotv_v, dudt, nlhalf
371  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: vel_tot !(noeud)
372  TYPE(mesh_type), TARGET :: mesh
373  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: v1m, vit !V(noeud, type)
374  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: p, phalf
375  REAL(KIND=8), INTENT(IN) :: dt
376  INTEGER, INTENT(IN) :: mode
377 !TEST LES LC October 28, 2014
378  REAL(KIND=8) :: vel_tot_max
379 !TEST LES LC October 28, 2014
380  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, fv, mult
381  REAL(KIND=8), DIMENSION(2,6,mesh%gauss%l_G) :: grad
382  REAL(KIND=8), DIMENSION(6,mesh%gauss%l_G) :: vitloc
383  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
384  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
385  REAL(KIND=8), DIMENSION(mesh%dom_me) :: visc_plot
386  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_w) :: v14_loc, v23_loc
387  INTEGER, DIMENSION(2*mesh%gauss%n_w) :: idxm_2
388  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v5_loc, v6_loc
389  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm_1
390  REAL(KIND=8), DIMENSION(mesh%dom_me*mesh%gauss%l_G,6) :: rhs_gauss
391  REAL(KIND=8), DIMENSION(6) :: visc1, visc2
392  REAL(KIND=8) :: ray, div1, div2, cfl, vloc, normal_vit
393 
394  REAL(KIND=8), ALLOCATABLE, DIMENSION(:), SAVE :: h
395  REAL(KIND=8), SAVE :: coeff_ed_st, r_eff, &
396  visc_eff, surf, nu_loc, coeff_visc_ordre_un
397  LOGICAL, SAVE :: once = .true.
398  REAL(KIND=8), DIMENSION(6,mesh%dom_me) :: viscosity
399  REAL(KIND=8), DIMENSION(6) :: norm_vit
400  REAL(KIND=8) :: type_fe
401  INTEGER :: m, l , i , ni , index, index2, type, k
402  INTEGER :: ms, nw, ix, ki, iglob
403 !!$ WARNING FL Variables removed (used for edge_stab)
404 !!$ INTEGER :: cotei, ls
405 !!$ REAL(KIND=8) :: dul, ed_st, h2
406 !!$ INTEGER, DIMENSION(mesh%gauss%n_w,2) :: jji_loc
407 !!$ REAL(KIND=8), DIMENSION(mesh%gauss%n_w,2) :: u0loci, uloci
408 !!$ REAL(KIND=8), DIMENSION(mesh%gauss%n_w,mesh%gauss%l_Gs, 2) :: dwni_loc
409 !!$ WARNING FL Variables removed (used for edge_stab)
410 #include "petsc/finclude/petsc.h"
411  vec :: vb_145, vb_236
412  !PetscErrorCode :: ierr
413 
414  !CALL VecSet(vb_3_145, 0.d0, ierr)
415  !CALL VecSet(vb_3_236, 0.d0, ierr)
416 
417  IF (once) THEN
418  once =.false.
419 
420  IF (.NOT.inputs%LES) THEN
421  inputs%LES_coeff1=0.d0
422  inputs%LES_coeff2=0.d0
423  inputs%LES_coeff3=0.d0
424  inputs%LES_coeff4=0.d0
425  END IF
426 
427  surf = 0.d0
428  DO m = 1, mesh%dom_me
429  DO l = 1, mesh%gauss%l_G
430  ray = sum(mesh%rr(1,mesh%jj(:,m))*mesh%gauss%ww(:,l))
431  surf = surf + mesh%gauss%rj(l,m)*ray
432  END DO
433  END DO
434 
435  IF (mesh%gauss%n_w==3) THEN
436  type_fe = 1
437  ELSE
438  type_fe = 2
439  END IF
440  coeff_ed_st = inputs%LES_coeff3*0.02d0/type_fe
441  coeff_visc_ordre_un = inputs%LES_coeff4
442  IF (mesh%edge_stab) THEN
443  ALLOCATE(h(mesh%mi))
444  DO ms = 1, mesh%mi
445  h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
446  END DO
447  END IF
448  END IF
449 
450  !ATTENTION: JLG Jan 25 2010
451  !ATTENTION: inputs%LES_coeff1 is assumed to be of the order of the convective velocity
452  !that simplifies the semi-implicit treatment of the LES viscosity
453 !!$ normal_vit = MAXVAL(vel_tot)
454 !TEST LES LC October 28, 2014
455  normal_vit = vel_tot_max
456 !TEST LES LC October 28, 2014
457  DO TYPE = 1, 6
458  norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
459  END DO
460  !ATTENTION: JLG Jan 25 2010
461  r_eff = 0.d0
462  cfl = 0
463  index = 0
464  index2 = 0
465  IF (inputs%LES) THEN
466  DO m = 1, mesh%dom_me
467  j_loc = mesh%jj(:,m)
468  vloc = maxval(vel_tot(j_loc))
469  cfl = max(vloc*dt/mesh%hloc(m),cfl)
470  visc1 = 0
471  visc2 = 0
472  DO l = 1, mesh%gauss%l_G
473  index2 = index2 +1
474  dw_loc = mesh%gauss%dw(:,:,l,m)
475 
476  !--------radius of gauss point
477  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
478  !fs Source term
479  !ft Time derivative
480  !fp Pressure gradient
481  !--------calcul de la premiere composante
482  ft(1) = sum(dudt(j_loc,1) *mesh%gauss%ww(:,l))
483  fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
484  !--------calcul de la seconde composante
485  ft(2) = sum(dudt(j_loc,2) * mesh%gauss%ww(:,l))
486  fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
487  !--------calcul de la troisieme composante
488  ft(3) = sum(dudt(j_loc,3) *mesh%gauss%ww(:,l))
489  fp(3) = sum(phalf(j_loc,2)*mesh%gauss%ww(:,l))*mode/ray
490  !--------calcul de la quatrieme composante
491  ft(4) = sum(dudt(j_loc,4) *mesh%gauss%ww(:,l))
492  fp(4) = -sum(phalf(j_loc,1)*mesh%gauss%ww(:,l))*mode/ray
493  !--------calcul de la cinquieme composante
494  ft(5) = sum(dudt(j_loc,5) *mesh%gauss%ww(:,l))
495  fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
496  !--------calcul de la sixieme composante
497  ft(6) = sum(dudt(j_loc,6) *mesh%gauss%ww(:,l))
498  fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
499  !-------calcul du second membre pour le terme nonlineaire------------------------
500 
501  visc1= max(visc1,abs(ft+fp+nlhalf(index2,:)))
502 
503  !--------Calcul du gradient de la vitesse
504  DO TYPE = 1, 6
505  DO k = 1 ,2
506  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
507  END DO
508  vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
509  END DO
510 
511  !--------Calcul de la divergence
512  div1 = abs(grad(1,1,l) + vitloc(1,l)/ray + grad(2,5,l) + vitloc(4,l)*mode/ray)
513  div2 = abs(grad(1,2,l) + vitloc(2,l)/ray + grad(2,6,l) - vitloc(3,l)*mode/ray)
514  visc2(1) = max(visc2(1),div1)
515  visc2(2) = max(visc2(2),div2)
516  END DO
517  visc2(4) = visc2(1); visc2(5) = visc2(1)
518  visc2(3) = visc2(2); visc2(6) = visc2(2)
519 
520  nu_loc = 0.d0
521  DO TYPE = 1, 6
522  !visc1(type) = MAX(visc1(type)/normal_vit,visc2(type))
523  !visc1(type) = MAX(visc1(type),visc2(type)*normal_vit)/MAXVAL(ABS(vitloc(type,:)))
524  visc1(type) = max(visc1(type),2*visc2(type)*normal_vit)/norm_vit(type)
525  visc_eff = mesh%hloc(m)*min(coeff_visc_ordre_un*normal_vit,inputs%LES_coeff2*mesh%hloc(m)*visc1(type))
526  nu_loc = nu_loc + visc_eff
527  !======Semi-implicit version==========
528  viscosity(type,m) = inputs%LES_coeff1*mesh%hloc(m) - visc_eff
529  !======Semi-implicit version==========
530  END DO
531  r_eff = r_eff + (nu_loc + dt**2*mesh%hloc(m))*mesh%hloc(m)**2*ray
532  visc_plot(m) = (nu_loc/6)/(coeff_visc_ordre_un*mesh%hloc(m)*normal_vit)
533  END DO
534  ELSE
535  cfl = 0.d0
536  viscosity = 0.d0
537  DO m = 1, mesh%dom_me
538  j_loc = mesh%jj(:,m)
539  vloc = maxval(vel_tot(j_loc))
540  cfl = max(vloc*dt/mesh%hloc(m),cfl)
541  END DO
542  END IF
543 
544 
545  !DO type = 1, 6
546  ! CALL average(mesh,viscosity(type,:))
547  !END DO
548 
549  nw = mesh%gauss%n_w
550  DO m = 1, mesh%dom_me
551  mult(:)= viscosity(:,m)
552  j_loc = mesh%jj(:,m)
553 
554  DO ni = 1, nw
555  i = j_loc(ni)
556  iglob = vv_3_la%loc_to_glob(3,i)
557  idxm_1(ni) = iglob-1
558  END DO
559 
560  DO ki = 1, 2
561  DO ni = 1, nw
562  i = j_loc(ni)
563  iglob = vv_3_la%loc_to_glob(ki,i)
564  ix = (ki-1)*nw+ni
565  idxm_2(ix) = iglob-1
566  END DO
567  END DO
568 
569  v14_loc = 0.d0
570  v23_loc = 0.d0
571  v5_loc = 0.d0
572  v6_loc = 0.d0
573  DO l = 1, mesh%gauss%l_G
574  index = index +1
575  dw_loc = mesh%gauss%dw(:,:,l,m)
576  DO TYPE = 1, 6
577  DO k = 1 ,2
578  grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
579  END DO
580  vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
581  END DO
582 
583  !===Compute radius of Gauss point
584  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
585 
586  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
587  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
588  fs(1) = sum(ff(j_loc,1) * mesh%gauss%ww(:,l))
589  ft(1) = sum(v1m(j_loc,1) * mesh%gauss%ww(:,l))
590  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
591  fv(1) = ((mode*vitloc(1,l)+vitloc(4,l))*mode +mode*vitloc(4,l)+vitloc(1,l))/ray**2
592  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
593  fs(2) = sum(ff(j_loc,2) * mesh%gauss%ww(:,l))
594  ft(2) = sum(v1m(j_loc,2) * mesh%gauss%ww(:,l))
595  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
596  fv(2) = ((mode*vitloc(2,l)-vitloc(3,l))*mode -mode*vitloc(3,l)+vitloc(2,l))/ray**2
597  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
598  fs(3) = sum(ff(j_loc,3) * mesh%gauss%ww(:,l))
599  ft(3) = sum(v1m(j_loc,3) * mesh%gauss%ww(:,l))
600  fp(3) = -sum(p(j_loc,2)*mesh%gauss%ww(:,l))/ray*mode
601  fv(3) = (-mode*vitloc(2,l)+vitloc(3,l) +(mode*vitloc(3,l)-vitloc(2,l))*mode)/ray**2
602  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
603  fs(4) = sum(ff(j_loc,4) * mesh%gauss%ww(:,l))
604  ft(4) = sum(v1m(j_loc,4) * mesh%gauss%ww(:,l))
605  fp(4) = sum(p(j_loc,1)*mesh%gauss%ww(:,l))/ray *mode
606  fv(4) = (mode*vitloc(1,l)+vitloc(4,l) +(mode*vitloc(4,l)+vitloc(1,l))*mode)/ray**2
607  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
608  fs(5) = sum(ff(j_loc,5) * mesh%gauss%ww(:,l))
609  ft(5) = sum(v1m(j_loc,5) * mesh%gauss%ww(:,l))
610  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
611  fv(5) = vitloc(5,l)*(mode/ray)**2
612  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
613  fs(6) = sum(ff(j_loc,6) * mesh%gauss%ww(:,l))
614  ft(6) = sum(v1m(j_loc,6) * mesh%gauss%ww(:,l))
615  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
616  fv(6) = vitloc(6,l)*(mode/ray)**2
617  !-------calcul du second membre pour le terme nonlineaire------------------------
618  !-------------------------------------------------------------------------------
619 
620  fv = mult*fv
621 
622  rhs_gauss(index, :) = (ft+fp+fs+fv-rotv_v(index,:))
623 
624  ENDDO
625 
626  ENDDO
627 
628  CALL rhs_3x3(mesh, vv_3_la, mode, rhs_gauss, vb_145, vb_236)
629 
630  END SUBROUTINE qs_ns_stab_3x3
631 
632  SUBROUTINE qs_01_div_hybrid_2006 (vv_mesh, pp_mesh, pp_LA, mode, gg, pb_1, pb_2)
633  !================================================
634  USE def_type_mesh
635  IMPLICIT NONE
636  TYPE(mesh_type) :: vv_mesh, pp_mesh
637  TYPE(petsc_csr_la) :: pp_la
638  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: gg
639  INTEGER , INTENT(IN) :: mode
640  REAL(KIND=8) :: x
641  REAL(KIND=8), DIMENSION(pp_mesh%gauss%n_w, vv_mesh%gauss%l_G) :: w_c
642  REAL(KIND=8), DIMENSION(pp_mesh%gauss%n_w) :: v1_loc, v2_loc
643  INTEGER, DIMENSION(pp_mesh%gauss%n_w) :: idxm
644  INTEGER :: m, l, n, i, nw, nwc, ni, iglob
645  REAL(KIND=8), DIMENSION(3,2) :: f
646  REAL(KIND=8) :: ray
647  INTEGER, DIMENSION(vv_mesh%gauss%n_w) :: j_loc
648  INTEGER, DIMENSION(pp_mesh%gauss%n_w) :: jc_loc
649 #include "petsc/finclude/petsc.h"
650  vec :: pb_1, pb_2
651  petscerrorcode :: ierr
652 
653  CALL vecset(pb_1, 0.d0, ierr)
654  CALL vecset(pb_2, 0.d0, ierr)
655 
656  nw = vv_mesh%gauss%n_w
657  DO l = 1, vv_mesh%gauss%l_G
658  !P1/P2
659  w_c(1,l) = vv_mesh%gauss%ww(1,l) + 0.5*(vv_mesh%gauss%ww(nw-1,l) + vv_mesh%gauss%ww(nw,l))
660  w_c(2,l) = vv_mesh%gauss%ww(2,l) + 0.5*(vv_mesh%gauss%ww(nw,l) + vv_mesh%gauss%ww(4,l))
661  w_c(3,l) = vv_mesh%gauss%ww(3,l) + 0.5*(vv_mesh%gauss%ww(4,l) + vv_mesh%gauss%ww(nw-1,l))
662  END DO
663 
664  nwc = pp_mesh%gauss%n_w
665  DO m = 1, vv_mesh%dom_me
666  j_loc(:) = vv_mesh%jj(:,m)
667  jc_loc(:)= pp_mesh%jj(:,m)
668 
669  DO ni = 1, nwc
670  i = jc_loc(ni)
671  iglob = pp_la%loc_to_glob(1,i)
672  idxm(ni) = iglob-1
673  END DO
674 
675  v1_loc = 0.d0
676  v2_loc = 0.d0
677  DO l = 1, vv_mesh%gauss%l_G
678 
679  !--------radius of P2 Gauss point
680  ray = 0
681  DO n = 1, nw; i = vv_mesh%jj(n,m)
682  ray = ray + vv_mesh%rr(1,i)*vv_mesh%gauss%ww(n,l)
683  END DO
684  !----------------------
685 
686  !calcul de la divergence sur les cos
687  f(1,1) = (ray*sum(gg(j_loc,1)*vv_mesh%gauss%dw(1,:,l,m)) &
688  + sum(gg(j_loc,1)*vv_mesh%gauss%ww(:,l)))
689  f(2,1) = mode*sum(gg(j_loc,4)*vv_mesh%gauss%ww(:,l))
690  f(3,1) = sum(gg(j_loc,5)*vv_mesh%gauss%dw(2,:,l,m)) * ray
691 
692  !calcul de la divergence sur les sin
693  f(1,2) = (ray*sum(gg(j_loc,2)*vv_mesh%gauss%dw(1,:,l,m)) &
694  + sum(gg(j_loc,2)*vv_mesh%gauss%ww(:,l)))
695  f(2,2) =-mode*sum(gg(j_loc,3)*vv_mesh%gauss%ww(:,l))
696  f(3,2) = sum(gg(j_loc,6)*vv_mesh%gauss%dw(2,:,l,m)) * ray
697 
698  f = f *vv_mesh%gauss%rj(l,m)
699 
700 !!$ DO j=1,2
701 !!$ x = f(1,j)+f(2,j)+f(3,j)
702 !!$ DO n=1, nwc
703 !!$ u0_c(jj_c(n,m),j) = u0_c(jj_c(n,m),j) + w_c(n,l)*x
704 !!$ ENDDO
705 !!$ ENDDO
706 
707  x = f(1,1)+f(2,1)+f(3,1)
708  DO ni = 1, nwc
709  v1_loc(ni) = v1_loc(ni) + w_c(ni,l)*x
710  END DO
711 
712  x = f(1,2)+f(2,2)+f(3,2)
713  DO ni = 1, nwc
714  v2_loc(ni) = v2_loc(ni) + w_c(ni,l)*x
715  END DO
716 
717  ENDDO
718  CALL vecsetvalues(pb_1, nwc, idxm, v1_loc, add_values, ierr)
719  CALL vecsetvalues(pb_2, nwc, idxm, v2_loc, add_values, ierr)
720  ENDDO
721  CALL vecassemblybegin(pb_1,ierr)
722  CALL vecassemblyend(pb_1,ierr)
723 
724  CALL vecassemblybegin(pb_2,ierr)
725  CALL vecassemblyend(pb_2,ierr)
726 
727  END SUBROUTINE qs_01_div_hybrid_2006
728 
729  SUBROUTINE qs_00 (mesh, LA, ff, vect)
730  !=================================
731  USE def_type_mesh
732  IMPLICIT NONE
733  TYPE(mesh_type), TARGET :: mesh
734  type(petsc_csr_la) :: la
735  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff
736  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: ff_loc
737  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
738  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v_loc
739  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm
740  INTEGER :: i, m, l, ni, iglob
741  REAL(KIND=8) :: fl, ray
742 #include "petsc/finclude/petsc.h"
743  vec :: vect
744  petscerrorcode :: ierr
745 
746  CALL vecset(vect, 0.d0, ierr)
747 
748  DO m = 1, mesh%dom_me
749  jj_loc = mesh%jj(:,m)
750  ff_loc = ff(jj_loc)
751  DO ni = 1, mesh%gauss%n_w
752  i = mesh%jj(ni,m)
753  iglob = la%loc_to_glob(1,i)
754  idxm(ni) = iglob-1
755  END DO
756 
757  v_loc = 0.d0
758  DO l = 1, mesh%gauss%l_G
759  ray = 0
760  DO ni = 1, mesh%gauss%n_w
761  i = jj_loc(ni)
762  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
763  END DO
764  fl = sum(ff_loc*mesh%gauss%ww(:,l))*mesh%gauss%rj(l,m)*ray
765  DO ni = 1, mesh%gauss%n_w
766  v_loc(ni) = v_loc(ni) + mesh%gauss%ww(ni,l) * fl
767  END DO
768  ENDDO
769  CALL vecsetvalues(vect, mesh%gauss%n_w, idxm, v_loc, add_values, ierr)
770  ENDDO
771  CALL vecassemblybegin(vect,ierr)
772  CALL vecassemblyend(vect,ierr)
773  END SUBROUTINE qs_00
774 
775  SUBROUTINE qs_00_gauss (mesh, LA, heat_capa, ff, ff_gauss, vect)
776  !=================================
777  USE def_type_mesh
778  IMPLICIT NONE
779  TYPE(mesh_type), TARGET :: mesh
780  type(petsc_csr_la) :: la
781  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: heat_capa
782  REAL(KIND=8), DIMENSION(:), INTENT(IN) :: ff, ff_gauss
783  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: ff_loc
784  INTEGER, DIMENSION(mesh%gauss%n_w) :: jj_loc
785  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v_loc
786  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm
787  INTEGER :: i, m, l, ni, iglob, index
788  REAL(KIND=8) :: fl, ray
789 #include "petsc/finclude/petsc.h"
790  vec :: vect
791  petscerrorcode :: ierr
792 
793  CALL vecset(vect, 0.d0, ierr)
794 
795  index = 0
796  DO m = 1, mesh%dom_me
797  jj_loc = mesh%jj(:,m)
798  ff_loc = ff(jj_loc)
799  DO ni = 1, mesh%gauss%n_w
800  i = mesh%jj(ni,m)
801  iglob = la%loc_to_glob(1,i)
802  idxm(ni) = iglob-1
803  END DO
804 
805  v_loc = 0.d0
806  DO l = 1, mesh%gauss%l_G
807  index =index + 1
808  ray = 0
809  DO ni = 1, mesh%gauss%n_w
810  i = jj_loc(ni)
811  ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
812  END DO
813  fl = (heat_capa(m) * sum(ff_loc*mesh%gauss%ww(:,l)) + ff_gauss(index))*mesh%gauss%rj(l,m)*ray
814  DO ni = 1, mesh%gauss%n_w
815  v_loc(ni) = v_loc(ni) + mesh%gauss%ww(ni,l) * fl
816  END DO
817  ENDDO
818  CALL vecsetvalues(vect, mesh%gauss%n_w, idxm, v_loc, add_values, ierr)
819  ENDDO
820  CALL vecassemblybegin(vect,ierr)
821  CALL vecassemblyend(vect,ierr)
822  END SUBROUTINE qs_00_gauss
823 
824  SUBROUTINE qs_ns_momentum_stab_new(mesh,vv_1_LA,vv_2_LA,mode,ff,V1m,P,&
825  vb_2_14,vb_2_23,vb_1_5,vb_1_6, rotb_b, tensor, tensor_surface_gauss,&
826  stab, momentum, momentum_les, visc_grad_vel, visc_entro)
827  !=================================
828  !RHS for Navier-Stokes with momentum
829  USE def_type_mesh
830  USE chaine_caractere
831  USE my_util
832  USE input_data
833  !USE sub_plot
834  IMPLICIT NONE
835  TYPE(petsc_csr_la) :: vv_1_la,vv_2_la
836  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotb_b
837  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tensor !(node)
838  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tensor_surface_gauss !(gauss)
839  TYPE(mesh_type), TARGET :: mesh
840  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: v1m, momentum, momentum_les !V(node, type)
841  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: p
842  REAL(KIND=8), INTENT(IN) :: stab
843  INTEGER, INTENT(IN) :: mode
844  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: visc_grad_vel !V(r/th/z, gauss, type)
845  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: visc_entro
846  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb, mult
847  REAL(KIND=8), DIMENSION(6) :: fnl, fvgm, fvgm_les, fvgu
848  REAL(KIND=8), DIMENSION(2,6,mesh%gauss%l_G) :: grad_mom
849  REAL(KIND=8), DIMENSION(6,mesh%gauss%l_G) :: momloc, momloc_les
850  REAL(KIND=8), DIMENSION(3,6,mesh%gauss%l_G) :: tensor_loc
851  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
852  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
853  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_w) :: v14_loc, v23_loc
854  INTEGER, DIMENSION(2*mesh%gauss%n_w) :: idxm_2
855  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v5_loc, v6_loc
856  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm_1
857  REAL(KIND=8) :: ray
858  REAL(KIND=8), DIMENSION(mesh%dom_me) :: stab_loc
859  LOGICAL, SAVE :: once = .true.
860  REAL(KIND=8), SAVE :: type_fe
861  INTEGER :: m, l , i , ni , index, type, k
862  INTEGER :: nw, ix, ki, iglob
863 #include "petsc/finclude/petsc.h"
864  vec :: vb_2_14,vb_2_23,vb_1_5,vb_1_6
865  petscerrorcode :: ierr
866 
867  CALL vecset(vb_2_14, 0.d0, ierr)
868  CALL vecset(vb_2_23, 0.d0, ierr)
869  CALL vecset(vb_1_5, 0.d0, ierr)
870  CALL vecset(vb_1_6, 0.d0, ierr)
871 
872  IF (once) THEN
873  once =.false.
874 
875  IF (.NOT.inputs%LES) THEN
876  inputs%LES_coeff1=0.d0
877  END IF
878 
879  IF (mesh%gauss%n_w==3) THEN
880  type_fe = 1
881  ELSE
882  type_fe = 2
883  END IF
884  END IF ! end once
885 
886  IF (inputs%LES) THEN
887  DO m = 1, mesh%dom_me
888  stab_loc(m) = stab + inputs%LES_coeff1*mesh%hloc(m)
889  END DO
890  ELSE
891  stab_loc(:) = stab
892  END IF
893 
894  index = 0
895  nw = mesh%gauss%n_w
896 
897  DO m = 1, mesh%dom_me
898  mult(:)= - visc_entro(m,1) !visc_entro is the same in type 1 or 2
899 
900  j_loc = mesh%jj(:,m)
901 
902  DO ni = 1, nw
903  i = j_loc(ni)
904  iglob = vv_1_la%loc_to_glob(1,i)
905  idxm_1(ni) = iglob-1
906  END DO
907 
908  DO ki = 1, 2
909  DO ni = 1, nw
910  i = j_loc(ni)
911  iglob = vv_2_la%loc_to_glob(ki,i)
912  ix = (ki-1)*nw+ni
913  idxm_2(ix) = iglob-1
914  END DO
915  END DO
916 
917  v14_loc = 0.d0
918  v23_loc = 0.d0
919  v5_loc = 0.d0
920  v6_loc = 0.d0
921  DO l = 1, mesh%gauss%l_G
922  index = index +1
923  dw_loc = mesh%gauss%dw(:,:,l,m)
924  DO TYPE = 1, 6
925  DO k = 1, 3
926  tensor_loc(k,type,l) = sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
927  + tensor_surface_gauss(k,index,type)
928  END DO
929  END DO
930 
931  !===Compute radius of Gauss point
932  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
933 
934  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
935  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
936  fs(1) = sum(ff(j_loc,1) * mesh%gauss%ww(:,l))
937  ft(1) = sum(v1m(j_loc,1) * mesh%gauss%ww(:,l))
938  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
939  fnl(1) = (-mode*tensor_loc(1,4,l) + tensor_loc(2,3,l))/ray
940  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
941  fs(2) = sum(ff(j_loc,2) * mesh%gauss%ww(:,l))
942  ft(2) = sum(v1m(j_loc,2) * mesh%gauss%ww(:,l))
943  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
944  fnl(2) = (mode*tensor_loc(1,3,l) + tensor_loc(2,4,l))/ray
945  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
946  fs(3) = sum(ff(j_loc,3) * mesh%gauss%ww(:,l))
947  ft(3) = sum(v1m(j_loc,3) * mesh%gauss%ww(:,l))
948  fp(3) = -sum(p(j_loc,2)*mesh%gauss%ww(:,l))/ray*mode
949  fnl(3) = (-tensor_loc(1,3,l) - mode*tensor_loc(2,4,l))/ray
950  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
951  fs(4) = sum(ff(j_loc,4) * mesh%gauss%ww(:,l))
952  ft(4) = sum(v1m(j_loc,4) * mesh%gauss%ww(:,l))
953  fp(4) = sum(p(j_loc,1)*mesh%gauss%ww(:,l))/ray *mode
954  fnl(4) = (-tensor_loc(1,4,l) + mode*tensor_loc(2,3,l))/ray
955  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
956  fs(5) = sum(ff(j_loc,5) * mesh%gauss%ww(:,l))
957  ft(5) = sum(v1m(j_loc,5) * mesh%gauss%ww(:,l))
958  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
959  fnl(5) = -tensor_loc(3,4,l)*mode/ray
960  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
961  fs(6) = sum(ff(j_loc,6) * mesh%gauss%ww(:,l))
962  ft(6) = sum(v1m(j_loc,6) * mesh%gauss%ww(:,l))
963  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
964  fnl(6) = tensor_loc(3,3,l)*mode/ray
965  !-------calcul du second membre pour le terme nonlineaire------------------------
966  !-------------------------------------------------------------------------------
967 
968  IF (inputs%if_level_set) THEN
969  DO TYPE = 1, 6
970 
971  DO k = 1 ,2
972  grad_mom(k,type,l) = stab_loc(m)* sum(momentum(j_loc,type)*dw_loc(k,:)) &
973  + mult(type)*sum(momentum_les(j_loc,type)*dw_loc(k,:))
974  END DO
975 
976  momloc(type,l) = sum(momentum(j_loc,type)*mesh%gauss%ww(:,l))
977  momloc_les(type,l) = sum(momentum_les(j_loc,type)*mesh%gauss%ww(:,l))
978 
979  END DO
980  fvgm(1) = ((mode*momloc(1,l)+momloc(4,l))*mode + mode*momloc(4,l)+momloc(1,l))/ray**2
981  fvgm(2) = ((mode*momloc(2,l)-momloc(3,l))*mode - mode*momloc(3,l)+momloc(2,l))/ray**2
982  fvgm(3) = (-mode*momloc(2,l)+momloc(3,l) + (mode*momloc(3,l)-momloc(2,l))*mode)/ray**2
983  fvgm(4) = (mode*momloc(1,l)+momloc(4,l) + (mode*momloc(4,l)+momloc(1,l))*mode)/ray**2
984  fvgm(5) = momloc(5,l)*(mode/ray)**2
985  fvgm(6) = momloc(6,l)*(mode/ray)**2
986 
987  fvgm_les(1) = ((mode*momloc_les(1,l)+momloc_les(4,l))*mode + mode*momloc_les(4,l)+momloc_les(1,l))/ray**2
988  fvgm_les(2) = ((mode*momloc_les(2,l)-momloc_les(3,l))*mode - mode*momloc_les(3,l)+momloc_les(2,l))/ray**2
989  fvgm_les(3) = (-mode*momloc_les(2,l)+momloc_les(3,l) + (mode*momloc_les(3,l)-momloc_les(2,l))*mode)/ray**2
990  fvgm_les(4) = (mode*momloc_les(1,l)+momloc_les(4,l) + (mode*momloc_les(4,l)+momloc_les(1,l))*mode)/ray**2
991  fvgm_les(5) = momloc_les(5,l)*(mode/ray)**2
992  fvgm_les(6) = momloc_les(6,l)*(mode/ray)**2
993 
994  fvgm = stab_loc(m)*fvgm + mult*fvgm_les
995 
996  fvgu(1) = (-visc_grad_vel(1,index,4)*mode + visc_grad_vel(2,index,3))/ray
997  fvgu(2) = (visc_grad_vel(1,index,3)*mode + visc_grad_vel(2,index,4))/ray
998  fvgu(3) = (-visc_grad_vel(1,index,3) - visc_grad_vel(2,index,4)*mode)/ray
999  fvgu(4) = (-visc_grad_vel(1,index,4) + visc_grad_vel(2,index,3)*mode)/ray
1000  fvgu(5) = -visc_grad_vel(3,index,4)*mode/ray
1001  fvgu(6) = visc_grad_vel(3,index,3)*mode/ray
1002 
1003  DO TYPE = 1, 6
1004  grad_mom(:,type,l) = grad_mom(:,type,l)*ray*mesh%gauss%rj(l,m)
1005  END DO
1006  tensor_loc(:,:,l) = tensor_loc(:,:,l) + visc_grad_vel(:,index,:)
1007 
1008  ELSE
1009  grad_mom = 0.d0
1010  fvgm = 0.d0
1011  fvgu = 0.d0
1012  END IF
1013 
1014  ! if NOT level_set then :
1015  ! rhs = (BDF2(n&n_m1 terms) - Grad_pressure + source_in_ns + precession + lorentz)*test_function
1016  ! + tensor_explicit:Grad(test_function) (tensor by FFT)
1017 
1018  ! if level_set then :
1019  !rhs = (BDF2(n&n_m1 terms) - Grad_pressure + source_in_ns + precession + lorentz)*test_function
1020  ! + (tensor_explicit + visc_grad_vel):Grad(test_function) (tensor and visc_grad_vel by FFT)
1021  ! + (LES + stab*grad_mom):GRAD(test_function)
1022 
1023  smb = (ft+fp+fs+rotb_b(index,:)+fnl+fvgm+fvgu)*ray*mesh%gauss%rj(l,m)
1024 
1025  ki = 1
1026  DO ni = 1, nw
1027  ix = (ki-1)*nw+ni
1028  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(1) + sum(dw_loc(:,ni)*grad_mom(:,1,l)) &
1029  + (dw_loc(1,ni)*tensor_loc(1,1,l) + dw_loc(2,ni)*tensor_loc(1,5,l))*ray*mesh%gauss%rj(l,m)
1030  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(2) + sum(dw_loc(:,ni)*grad_mom(:,2,l)) &
1031  + (dw_loc(1,ni)*tensor_loc(1,2,l) + dw_loc(2,ni)*tensor_loc(1,6,l))*ray*mesh%gauss%rj(l,m)
1032  v5_loc(ix) = v5_loc(ix) + mesh%gauss%ww(ni,l)*smb(5) + sum(dw_loc(:,ni)*grad_mom(:,5,l)) &
1033  + (dw_loc(1,ni)*tensor_loc(3,1,l) + dw_loc(2,ni)*tensor_loc(3,5,l))*ray*mesh%gauss%rj(l,m)
1034  v6_loc(ix) = v6_loc(ix) + mesh%gauss%ww(ni,l)*smb(6) + sum(dw_loc(:,ni)*grad_mom(:,6,l)) &
1035  + (dw_loc(1,ni)*tensor_loc(3,2,l) + dw_loc(2,ni)*tensor_loc(3,6,l))*ray*mesh%gauss%rj(l,m)
1036  END DO
1037 
1038  ki = 2
1039  DO ni = 1, nw
1040  ix = (ki-1)*nw+ni
1041  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(4) + sum(dw_loc(:,ni)*grad_mom(:,4,l)) &
1042  + (dw_loc(1,ni)*tensor_loc(2,2,l) + dw_loc(2,ni)*tensor_loc(2,6,l))*ray*mesh%gauss%rj(l,m)
1043  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(3) + sum(dw_loc(:,ni)*grad_mom(:,3,l)) &
1044  + (dw_loc(1,ni)*tensor_loc(2,1,l) + dw_loc(2,ni)*tensor_loc(2,5,l))*ray*mesh%gauss%rj(l,m)
1045  END DO
1046  ENDDO
1047 
1048  CALL vecsetvalues(vb_2_14, 2*nw, idxm_2, v14_loc, add_values, ierr)
1049  CALL vecsetvalues(vb_2_23, 2*nw, idxm_2, v23_loc, add_values, ierr)
1050  CALL vecsetvalues(vb_1_5, nw, idxm_1, v5_loc, add_values, ierr)
1051  CALL vecsetvalues(vb_1_6, nw, idxm_1, v6_loc, add_values, ierr)
1052  ENDDO
1053 
1054  IF (inputs%LES) THEN
1055  IF (mesh%edge_stab) THEN
1056  CALL error_petsc('BUG in qs_ns_stab_new: terms with edge_stab not yet assembled')
1057  END IF
1058  END IF
1059 
1060  CALL vecassemblybegin(vb_2_14,ierr)
1061  CALL vecassemblyend(vb_2_14,ierr)
1062  CALL vecassemblybegin(vb_2_23,ierr)
1063  CALL vecassemblyend(vb_2_23,ierr)
1064  CALL vecassemblybegin(vb_1_5,ierr)
1065  CALL vecassemblyend(vb_1_5,ierr)
1066  CALL vecassemblybegin(vb_1_6,ierr)
1067  CALL vecassemblyend(vb_1_6,ierr)
1068 
1069  END SUBROUTINE qs_ns_momentum_stab_new
1070 
1071  SUBROUTINE qs_ns_momentum_stab_3x3(mesh,vv_3_LA,mode,ff,V1m,P,&
1072  vb_3_145,vb_3_236, rotb_b, tensor, tensor_surface_gauss,&
1073  stab, momentum, visc_grad_vel)
1074  !=================================
1075  !RHS for Navier-Stokes
1076  USE def_type_mesh
1077  USE chaine_caractere
1078  USE my_util
1079  USE input_data
1080  !USE sub_plot
1081  IMPLICIT NONE
1082  TYPE(petsc_csr_la) :: vv_3_la
1083  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff, rotb_b
1084  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tensor !(node)
1085  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tensor_surface_gauss !(gauss)
1086  TYPE(mesh_type), TARGET :: mesh
1087  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: v1m, momentum !V(noeud, type)
1088  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: p
1089  REAL(KIND=8), INTENT(IN) :: stab
1090  INTEGER, INTENT(IN) :: mode
1091  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: visc_grad_vel !V(r/th/z, gauss, type)
1092  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, smb, fnl
1093  REAL(KIND=8), DIMENSION(6) :: fvgm, fvgu, fvgmt
1094  REAL(KIND=8), DIMENSION(2,6,mesh%gauss%l_G) :: grad_mom, grad_t_mom
1095  REAL(KIND=8), DIMENSION(6,mesh%gauss%l_G) :: momloc
1096  REAL(KIND=8), DIMENSION(3,6,mesh%gauss%l_G) :: tensor_loc
1097  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1098  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
1099  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_w) :: v14_loc, v23_loc
1100  INTEGER, DIMENSION(2*mesh%gauss%n_w) :: idxm_2
1101  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v5_loc, v6_loc
1102  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm_1
1103  REAL(KIND=8) :: ray
1104  REAL(KIND=8), DIMENSION(mesh%dom_me) :: stab_loc
1105  LOGICAL, SAVE :: once = .true.
1106  REAL(KIND=8), SAVE :: type_fe
1107  INTEGER :: m, l , i , ni , index, type, k
1108  INTEGER :: nw, ix, ki, iglob
1109 #include "petsc/finclude/petsc.h"
1110  vec :: vb_3_145, vb_3_236
1111  petscerrorcode :: ierr
1112 
1113  CALL vecset(vb_3_145, 0.d0, ierr)
1114  CALL vecset(vb_3_236, 0.d0, ierr)
1115 
1116  IF (once) THEN
1117  once =.false.
1118 
1119  IF (.NOT.inputs%LES) THEN
1120  inputs%LES_coeff1=0.d0
1121  END IF
1122 
1123  IF (mesh%gauss%n_w==3) THEN
1124  type_fe = 1
1125  ELSE
1126  type_fe = 2
1127  END IF
1128  END IF !end once
1129 
1130  stab_loc(:) = stab
1131 
1132  nw = mesh%gauss%n_w
1133  index = 0
1134 
1135  DO m = 1, mesh%dom_me
1136 
1137  j_loc = mesh%jj(:,m)
1138 
1139  DO ni = 1, nw
1140  i = j_loc(ni)
1141  iglob = vv_3_la%loc_to_glob(3,i)
1142  idxm_1(ni) = iglob-1
1143  END DO
1144 
1145  DO ki = 1, 2
1146  DO ni = 1, nw
1147  i = j_loc(ni)
1148  iglob = vv_3_la%loc_to_glob(ki,i)
1149  ix = (ki-1)*nw+ni
1150  idxm_2(ix) = iglob-1
1151  END DO
1152  END DO
1153 
1154  v14_loc = 0.d0
1155  v23_loc = 0.d0
1156  v5_loc = 0.d0
1157  v6_loc = 0.d0
1158  DO l = 1, mesh%gauss%l_G
1159  index = index +1
1160  dw_loc = mesh%gauss%dw(:,:,l,m)
1161  DO TYPE = 1, 6
1162  DO k = 1, 3
1163  tensor_loc(k,type,l) = sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
1164  + tensor_surface_gauss(k,index,type)
1165  END DO
1166  END DO
1167 
1168  !===Compute radius of Gauss point
1169  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
1170 
1171  !-calcul des seconds membres pr les termes de forçage, temporels et de gradP
1172  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1173  fs(1) = sum(ff(j_loc,1) * mesh%gauss%ww(:,l))
1174  ft(1) = sum(v1m(j_loc,1) * mesh%gauss%ww(:,l))
1175  fp(1) = -sum(p(j_loc,1)*dw_loc(1,:))
1176  fnl(1) = (-mode*tensor_loc(1,4,l) + tensor_loc(2,3,l))/ray
1177  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1178  fs(2) = sum(ff(j_loc,2) * mesh%gauss%ww(:,l))
1179  ft(2) = sum(v1m(j_loc,2) * mesh%gauss%ww(:,l))
1180  fp(2) = -sum(p(j_loc,2)*dw_loc(1,:))
1181  fnl(2) = (mode*tensor_loc(1,3,l) + tensor_loc(2,4,l))/ray
1182  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1183  fs(3) = sum(ff(j_loc,3) * mesh%gauss%ww(:,l))
1184  ft(3) = sum(v1m(j_loc,3) * mesh%gauss%ww(:,l))
1185  fp(3) = -sum(p(j_loc,2)*mesh%gauss%ww(:,l))/ray*mode
1186  fnl(3) = (-tensor_loc(1,3,l) - mode*tensor_loc(2,4,l))/ray
1187  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1188  fs(4) = sum(ff(j_loc,4) * mesh%gauss%ww(:,l))
1189  ft(4) = sum(v1m(j_loc,4) * mesh%gauss%ww(:,l))
1190  fp(4) = sum(p(j_loc,1)*mesh%gauss%ww(:,l))/ray *mode
1191  fnl(4) = (-tensor_loc(1,4,l) + mode*tensor_loc(2,3,l))/ray
1192  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1193  fs(5) = sum(ff(j_loc,5) * mesh%gauss%ww(:,l))
1194  ft(5) = sum(v1m(j_loc,5) * mesh%gauss%ww(:,l))
1195  fp(5) = -sum(p(j_loc,1)*dw_loc(2,:))
1196  fnl(5) = -tensor_loc(3,4,l)*mode/ray
1197  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1198  fs(6) = sum(ff(j_loc,6) * mesh%gauss%ww(:,l))
1199  ft(6) = sum(v1m(j_loc,6) * mesh%gauss%ww(:,l))
1200  fp(6) = -sum(p(j_loc,2)*dw_loc(2,:))
1201  fnl(6) = tensor_loc(3,3,l)*mode/ray
1202  !-------calcul du second membre pour le terme nonlineaire------------------------
1203  !-------------------------------------------------------------------------------
1204 
1205  IF (inputs%if_level_set) THEN
1206  DO TYPE = 1, 6
1207  DO k = 1 ,2
1208  grad_mom(k,type,l) = stab_loc(m)*sum(momentum(j_loc,type)*dw_loc(k,:))
1209  END DO
1210 
1211  momloc(type,l) = sum(momentum(j_loc,type)*mesh%gauss%ww(:,l))
1212  END DO
1213 
1214  fvgm(1) = ((mode*momloc(1,l)+momloc(4,l))*mode + mode*momloc(4,l)+momloc(1,l))/ray**2
1215  fvgm(2) = ((mode*momloc(2,l)-momloc(3,l))*mode - mode*momloc(3,l)+momloc(2,l))/ray**2
1216  fvgm(3) = (-mode*momloc(2,l)+momloc(3,l) + (mode*momloc(3,l)-momloc(2,l))*mode)/ray**2
1217  fvgm(4) = (mode*momloc(1,l)+momloc(4,l) + (mode*momloc(4,l)+momloc(1,l))*mode)/ray**2
1218  fvgm(5) = momloc(5,l)*(mode/ray)**2
1219  fvgm(6) = momloc(6,l)*(mode/ray)**2
1220 
1221  !===Compute Grad_T_mom . r or z derivative of test function
1222  grad_t_mom(1,1,l) = sum(momentum(j_loc,1)*dw_loc(1,:))
1223  grad_t_mom(2,1,l) = sum(momentum(j_loc,5)*dw_loc(1,:))
1224  grad_t_mom(1,2,l) = sum(momentum(j_loc,2)*dw_loc(1,:))
1225  grad_t_mom(2,2,l) = sum(momentum(j_loc,6)*dw_loc(1,:))
1226  grad_t_mom(1,3,l) = (mode*momloc(2,l) - momloc(3,l))/ray
1227  grad_t_mom(2,3,l) = mode*momloc(6,l)/ray
1228  grad_t_mom(1,4,l) = (-mode*momloc(1,l) - momloc(4,l))/ray
1229  grad_t_mom(2,4,l) = -mode*momloc(5,l)/ray
1230  grad_t_mom(1,5,l) = sum(momentum(j_loc,1)*dw_loc(2,:))
1231  grad_t_mom(2,5,l) = sum(momentum(j_loc,5)*dw_loc(2,:))
1232  grad_t_mom(1,6,l) = sum(momentum(j_loc,2)*dw_loc(2,:))
1233  grad_t_mom(2,6,l) = sum(momentum(j_loc,6)*dw_loc(2,:))
1234 
1235  !===Compute Grad_T_mom . NO (r or z derivative) of test function
1236  fvgmt(1) = -mode*sum(momentum(j_loc,4)*dw_loc(1,:))/ray &
1237  + (mode*momloc(4,l)+momloc(1,l))/ray**2
1238  fvgmt(2) = mode*sum(momentum(j_loc,3)*dw_loc(1,:))/ray &
1239  + (-mode*momloc(3,l)+momloc(2,l))/ray**2
1240  fvgmt(3) = -sum(momentum(j_loc,3)*dw_loc(1,:))/ray &
1241  + mode*(mode*momloc(3,l)-momloc(2,l))/ray**2
1242  fvgmt(4) = -sum(momentum(j_loc,4)*dw_loc(1,:))/ray &
1243  + mode*(mode*momloc(4,l)+momloc(1,l))/ray**2
1244  fvgmt(5) = -mode*sum(momentum(j_loc,4)*dw_loc(2,:))/ray
1245  fvgmt(6) = mode*sum(momentum(j_loc,3)*dw_loc(2,:))/ray
1246 
1247  !===update grad_mom
1248  grad_mom(:,:,l) = grad_mom(:,:,l) + stab*grad_t_mom(:,:,l)
1249  fvgm = stab_loc(m)*fvgm + stab*fvgmt
1250 
1251  fvgu(1) = (-visc_grad_vel(1,index,4)*mode + visc_grad_vel(2,index,3))/ray
1252  fvgu(2) = (visc_grad_vel(1,index,3)*mode + visc_grad_vel(2,index,4))/ray
1253  fvgu(3) = (-visc_grad_vel(1,index,3) - visc_grad_vel(2,index,4)*mode)/ray
1254  fvgu(4) = (-visc_grad_vel(1,index,4) + visc_grad_vel(2,index,3)*mode)/ray
1255  fvgu(5) = -visc_grad_vel(3,index,4)*mode/ray
1256  fvgu(6) = visc_grad_vel(3,index,3)*mode/ray
1257 
1258  DO TYPE = 1, 6
1259  grad_mom(:,type,l) = grad_mom(:,type,l)*ray*mesh%gauss%rj(l,m)
1260  END DO
1261  tensor_loc(:,:,l) = tensor_loc(:,:,l) + visc_grad_vel(:,index,:)
1262 
1263  ELSE
1264  grad_mom = 0.d0
1265  fvgm = 0.d0
1266  fvgu = 0.d0
1267  END IF
1268 
1269  ! if NOT level_set then :
1270  ! rhs = (BDF2(n&n_m1 terms) - Grad_pressure + source_in_ns + precession + lorentz)*test_function
1271  ! + tensor_explicit:Grad(test_function) (tensor by FFT)
1272 
1273  ! if level_set then :
1274  !rhs = (BDF2(n&n_m1 terms) - Grad_pressure + source_in_ns + precession + lorentz)*test_function
1275  ! + (tensor_explicit + visc_grad_vel):Grad(test_function) (tensor and visc_grad_vel by FFT)
1276  ! + (LES + stab*grad_mom):GRAD(test_function)
1277 
1278  smb = (ft+fp+fs+rotb_b(index,:)+fnl+fvgm+fvgu)*ray*mesh%gauss%rj(l,m)
1279 
1280  ki = 1
1281  DO ni = 1, nw
1282  ix = (ki-1)*nw+ni
1283  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(1) + sum(dw_loc(:,ni)*grad_mom(:,1,l)) &
1284  + (dw_loc(1,ni)*tensor_loc(1,1,l) + dw_loc(2,ni)*tensor_loc(1,5,l))*ray*mesh%gauss%rj(l,m)
1285  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(2) + sum(dw_loc(:,ni)*grad_mom(:,2,l)) &
1286  + (dw_loc(1,ni)*tensor_loc(1,2,l) + dw_loc(2,ni)*tensor_loc(1,6,l))*ray*mesh%gauss%rj(l,m)
1287  v5_loc(ix) = v5_loc(ix) + mesh%gauss%ww(ni,l)*smb(5) + sum(dw_loc(:,ni)*grad_mom(:,5,l)) &
1288  + (dw_loc(1,ni)*tensor_loc(3,1,l) + dw_loc(2,ni)*tensor_loc(3,5,l))*ray*mesh%gauss%rj(l,m)
1289  v6_loc(ix) = v6_loc(ix) + mesh%gauss%ww(ni,l)*smb(6) + sum(dw_loc(:,ni)*grad_mom(:,6,l)) &
1290  + (dw_loc(1,ni)*tensor_loc(3,2,l) + dw_loc(2,ni)*tensor_loc(3,6,l))*ray*mesh%gauss%rj(l,m)
1291  END DO
1292 
1293  ki = 2
1294  DO ni = 1, nw
1295  ix = (ki-1)*nw+ni
1296  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(4) + sum(dw_loc(:,ni)*grad_mom(:,4,l)) &
1297  + (dw_loc(1,ni)*tensor_loc(2,2,l) + dw_loc(2,ni)*tensor_loc(2,6,l))*ray*mesh%gauss%rj(l,m)
1298  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(3) + sum(dw_loc(:,ni)*grad_mom(:,3,l)) &
1299  + (dw_loc(1,ni)*tensor_loc(2,1,l) + dw_loc(2,ni)*tensor_loc(2,5,l))*ray*mesh%gauss%rj(l,m)
1300  END DO
1301 
1302  ENDDO
1303 
1304  CALL vecsetvalues(vb_3_145, 2*nw, idxm_2, v14_loc, add_values, ierr)
1305  CALL vecsetvalues(vb_3_236, 2*nw, idxm_2, v23_loc, add_values, ierr)
1306  CALL vecsetvalues(vb_3_145, nw, idxm_1, v5_loc, add_values, ierr)
1307  CALL vecsetvalues(vb_3_236, nw, idxm_1, v6_loc, add_values, ierr)
1308  ENDDO
1309 
1310  IF (inputs%LES) THEN
1311  IF (mesh%edge_stab) THEN
1312  CALL error_petsc('BUG in qs_ns_stab_new: terms with edge_stab not yet assembled')
1313  END IF
1314  END IF
1315 
1316  CALL vecassemblybegin(vb_3_145,ierr)
1317  CALL vecassemblyend(vb_3_145,ierr)
1318  CALL vecassemblybegin(vb_3_236,ierr)
1319  CALL vecassemblyend(vb_3_236,ierr)
1320 
1321  END SUBROUTINE qs_ns_momentum_stab_3x3
1322 
1323  SUBROUTINE qs_ns_mom_compute_residual_les(mesh,vv_1_LA,vv_2_LA,mode,ff,V1m,P,rotb_b, &
1324  tensor,tensor_gauss,vb_2_14,vb_2_23,vb_1_5,vb_1_6)
1325  !=================================
1326  USE def_type_mesh
1327  USE chaine_caractere
1328  USE my_util
1329  USE input_data
1330  IMPLICIT NONE
1331  TYPE(petsc_csr_la) :: vv_1_la,vv_2_la
1332  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff
1333  TYPE(mesh_type), TARGET :: mesh
1334  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: v1m !V(noeud, type)
1335  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tensor, tensor_gauss
1336  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: p
1337  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rotb_b
1338  INTEGER, INTENT(IN) :: mode
1339  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, ftensor, smb
1340  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1341  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
1342  REAL(KIND=8), DIMENSION(3,6,mesh%gauss%l_G) :: tensor_loc
1343  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_w) :: v14_loc, v23_loc
1344  INTEGER, DIMENSION(2*mesh%gauss%n_w) :: idxm_2
1345  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v5_loc, v6_loc
1346  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm_1
1347  REAL(KIND=8) :: ray
1348  INTEGER :: m, l , i , ni , index, type, k
1349  INTEGER :: nw, ix, ki, iglob
1350 #include "petsc/finclude/petsc.h"
1351  vec :: vb_2_14,vb_2_23,vb_1_5,vb_1_6
1352  petscerrorcode :: ierr
1353 
1354  CALL vecset(vb_2_14, 0.d0, ierr)
1355  CALL vecset(vb_2_23, 0.d0, ierr)
1356  CALL vecset(vb_1_5, 0.d0, ierr)
1357  CALL vecset(vb_1_6, 0.d0, ierr)
1358 
1359  index = 0
1360  nw = mesh%gauss%n_w
1361  DO m = 1, mesh%dom_me
1362  j_loc = mesh%jj(:,m)
1363 
1364  DO ni = 1, nw
1365  i = j_loc(ni)
1366  iglob = vv_1_la%loc_to_glob(1,i)
1367  idxm_1(ni) = iglob-1
1368  END DO
1369 
1370  DO ki = 1, 2
1371  DO ni = 1, nw
1372  i = j_loc(ni)
1373  iglob = vv_2_la%loc_to_glob(ki,i)
1374  ix = (ki-1)*nw+ni
1375  idxm_2(ix) = iglob-1
1376  END DO
1377  END DO
1378 
1379  v14_loc = 0.d0
1380  v23_loc = 0.d0
1381  v5_loc = 0.d0
1382  v6_loc = 0.d0
1383  DO l = 1, mesh%gauss%l_G
1384  index = index +1
1385  dw_loc = mesh%gauss%dw(:,:,l,m)
1386 
1387  !===Compute radius of Gauss point
1388  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
1389 
1390  !===Compute tensor on gauss points:
1391  DO TYPE = 1, 6
1392  DO k = 1, 3
1393  tensor_loc(k,type,l) = -sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
1394  + tensor_gauss(k,index,type)
1395  END DO
1396  END DO
1397 
1398  !===Compute residual part of source term, time derivative, pressure gradient without r&z derivatives
1399  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1400  fs(1) = sum(ff(j_loc,1) * mesh%gauss%ww(:,l))
1401  ft(1) = sum(v1m(j_loc,1) * mesh%gauss%ww(:,l))
1402  fp(1) = sum(p(j_loc,1)*dw_loc(1,:))
1403  ftensor(1) = -mode/ray*tensor_loc(1,4,l) + tensor_loc(2,3,l)/ray
1404  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1405  fs(2) = sum(ff(j_loc,2) * mesh%gauss%ww(:,l))
1406  ft(2) = sum(v1m(j_loc,2) * mesh%gauss%ww(:,l))
1407  fp(2) = sum(p(j_loc,2)*dw_loc(1,:))
1408  ftensor(2) = mode/ray*tensor_loc(1,3,l) + tensor_loc(2,4,l)/ray
1409  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1410  fs(3) = sum(ff(j_loc,3) * mesh%gauss%ww(:,l))
1411  ft(3) = sum(v1m(j_loc,3) * mesh%gauss%ww(:,l))
1412  fp(3) = sum(p(j_loc,2)*mesh%gauss%ww(:,l))/ray*mode
1413  ftensor(3) = -mode/ray*tensor_loc(2,4,l) - tensor_loc(1,3,l)/ray
1414  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1415  fs(4) = sum(ff(j_loc,4) * mesh%gauss%ww(:,l))
1416  ft(4) = sum(v1m(j_loc,4) * mesh%gauss%ww(:,l))
1417  fp(4) = -sum(p(j_loc,1)*mesh%gauss%ww(:,l))/ray *mode
1418  ftensor(4) = mode/ray*tensor_loc(2,3,l) - tensor_loc(1,4,l)/ray
1419  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1420  fs(5) = sum(ff(j_loc,5) * mesh%gauss%ww(:,l))
1421  ft(5) = sum(v1m(j_loc,5) * mesh%gauss%ww(:,l))
1422  fp(5) = sum(p(j_loc,1)*dw_loc(2,:))
1423  ftensor(5) = -mode/ray*tensor_loc(3,4,l)
1424  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1425  fs(6) = sum(ff(j_loc,6) * mesh%gauss%ww(:,l))
1426  ft(6) = sum(v1m(j_loc,6) * mesh%gauss%ww(:,l))
1427  fp(6) = sum(p(j_loc,2)*dw_loc(2,:))
1428  ftensor(6) = mode/ray*tensor_loc(3,3,l)
1429  !-------calcul du second membre pour le terme nonlineaire------------------------
1430  !-------------------------------------------------------------------------------
1431 
1432  smb = (ft+fp-fs+ftensor-rotb_b(index,:))*ray*mesh%gauss%rj(l,m)
1433 
1434  ki = 1
1435  DO ni = 1, nw
1436  ix = (ki-1)*nw+ni
1437  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(1) + &
1438  (dw_loc(1,ni)*tensor_loc(1,1,l) + dw_loc(2,ni)*tensor_loc(1,5,l))*ray*mesh%gauss%rj(l,m)
1439  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(2) + &
1440  (dw_loc(1,ni)*tensor_loc(1,2,l) + dw_loc(2,ni)*tensor_loc(1,6,l))*ray*mesh%gauss%rj(l,m)
1441  v5_loc(ix) = v5_loc(ix) + mesh%gauss%ww(ni,l)*smb(5) + &
1442  (dw_loc(1,ni)*tensor_loc(3,1,l) + dw_loc(2,ni)*tensor_loc(3,5,l))*ray*mesh%gauss%rj(l,m)
1443  v6_loc(ix) = v6_loc(ix) + mesh%gauss%ww(ni,l)*smb(6) + &
1444  (dw_loc(1,ni)*tensor_loc(3,2,l) + dw_loc(2,ni)*tensor_loc(3,6,l))*ray*mesh%gauss%rj(l,m)
1445  END DO
1446 
1447  ki = 2
1448  DO ni = 1, nw
1449  ix = (ki-1)*nw+ni
1450  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(4) + &
1451  (dw_loc(1,ni)*tensor_loc(2,2,l) + dw_loc(2,ni)*tensor_loc(2,6,l))*ray*mesh%gauss%rj(l,m)
1452  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(3) + &
1453  (dw_loc(1,ni)*tensor_loc(2,1,l) + dw_loc(2,ni)*tensor_loc(2,5,l))*ray*mesh%gauss%rj(l,m)
1454  END DO
1455  ENDDO
1456 
1457  CALL vecsetvalues(vb_2_14, 2*nw, idxm_2, v14_loc, add_values, ierr)
1458  CALL vecsetvalues(vb_2_23, 2*nw, idxm_2, v23_loc, add_values, ierr)
1459  CALL vecsetvalues(vb_1_5, nw, idxm_1, v5_loc, add_values, ierr)
1460  CALL vecsetvalues(vb_1_6, nw, idxm_1, v6_loc, add_values, ierr)
1461  ENDDO
1462 
1463  CALL vecassemblybegin(vb_2_14,ierr)
1464  CALL vecassemblyend(vb_2_14,ierr)
1465  CALL vecassemblybegin(vb_2_23,ierr)
1466  CALL vecassemblyend(vb_2_23,ierr)
1467  CALL vecassemblybegin(vb_1_5,ierr)
1468  CALL vecassemblyend(vb_1_5,ierr)
1469  CALL vecassemblybegin(vb_1_6,ierr)
1470  CALL vecassemblyend(vb_1_6,ierr)
1471 
1472  END SUBROUTINE qs_ns_mom_compute_residual_les
1473 
1474  SUBROUTINE qs_ns_mom_compute_residual_les_3x3(mesh,vv_3_LA,mode,ff,V1m,P,rotb_b, &
1475  tensor,tensor_gauss,vb_3_145,vb_3_236)
1476  !=================================
1477  USE def_type_mesh
1478  USE chaine_caractere
1479  USE my_util
1480  USE input_data
1481  IMPLICIT NONE
1482  TYPE(petsc_csr_la) :: vv_3_la
1483  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: ff
1484  TYPE(mesh_type), TARGET :: mesh
1485  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: v1m !V(noeud, type)
1486  REAL(KIND=8), DIMENSION(:,:,:), INTENT(IN) :: tensor, tensor_gauss
1487  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: p
1488  REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: rotb_b
1489  INTEGER, INTENT(IN) :: mode
1490  REAL(KIND=8), DIMENSION(6) :: fs , ft , fp, ftensor, smb
1491  INTEGER, DIMENSION(mesh%gauss%n_w) :: j_loc
1492  REAL(KIND=8), DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
1493  REAL(KIND=8), DIMENSION(3,6,mesh%gauss%l_G) :: tensor_loc
1494  REAL(KIND=8), DIMENSION(2*mesh%gauss%n_w) :: v14_loc, v23_loc
1495  INTEGER, DIMENSION(2*mesh%gauss%n_w) :: idxm_2
1496  REAL(KIND=8), DIMENSION(mesh%gauss%n_w) :: v5_loc, v6_loc
1497  INTEGER, DIMENSION(mesh%gauss%n_w) :: idxm_1
1498  REAL(KIND=8) :: ray
1499  INTEGER :: m, l , i , ni , index, type, k
1500  INTEGER :: nw, ix, ki, iglob
1501 #include "petsc/finclude/petsc.h"
1502  vec :: vb_3_145,vb_3_236
1503  petscerrorcode :: ierr
1504 
1505  CALL vecset(vb_3_145, 0.d0, ierr)
1506  CALL vecset(vb_3_236, 0.d0, ierr)
1507 
1508  index = 0
1509  nw = mesh%gauss%n_w
1510  DO m = 1, mesh%dom_me
1511  j_loc = mesh%jj(:,m)
1512 
1513  DO ni = 1, nw
1514  i = j_loc(ni)
1515  iglob = vv_3_la%loc_to_glob(3,i)
1516  idxm_1(ni) = iglob-1
1517  END DO
1518 
1519  DO ki = 1, 2
1520  DO ni = 1, nw
1521  i = j_loc(ni)
1522  iglob = vv_3_la%loc_to_glob(ki,i)
1523  ix = (ki-1)*nw+ni
1524  idxm_2(ix) = iglob-1
1525  END DO
1526  END DO
1527 
1528  v14_loc = 0.d0
1529  v23_loc = 0.d0
1530  v5_loc = 0.d0
1531  v6_loc = 0.d0
1532  DO l = 1, mesh%gauss%l_G
1533  index = index +1
1534  dw_loc = mesh%gauss%dw(:,:,l,m)
1535 
1536  !===Compute radius of Gauss point
1537  ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
1538 
1539  !===Compute tensor on gauss points:
1540  DO TYPE = 1, 6
1541  DO k = 1, 3
1542  tensor_loc(k,type,l) = -sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
1543  + tensor_gauss(k,index,type)
1544  END DO
1545  END DO
1546 
1547  !===Compute residual part of source term, time derivative, pressure gradient without r&z derivatives
1548  !--------calcul de la premiere composante : u0(1,:) <--> f(r,m,c)
1549  fs(1) = sum(ff(j_loc,1) * mesh%gauss%ww(:,l))
1550  ft(1) = sum(v1m(j_loc,1) * mesh%gauss%ww(:,l))
1551  fp(1) = sum(p(j_loc,1)*dw_loc(1,:))
1552  ftensor(1) = -mode/ray*tensor_loc(1,4,l) + tensor_loc(2,3,l)/ray
1553  !--------calcul de la seconde composante : u0(2,:) <--> f(r,m,s)
1554  fs(2) = sum(ff(j_loc,2) * mesh%gauss%ww(:,l))
1555  ft(2) = sum(v1m(j_loc,2) * mesh%gauss%ww(:,l))
1556  fp(2) = sum(p(j_loc,2)*dw_loc(1,:))
1557  ftensor(2) = mode/ray*tensor_loc(1,3,l) + tensor_loc(2,4,l)/ray
1558  !--------calcul de la troisieme composante : u0(3,:) <--> f(th,m,c)
1559  fs(3) = sum(ff(j_loc,3) * mesh%gauss%ww(:,l))
1560  ft(3) = sum(v1m(j_loc,3) * mesh%gauss%ww(:,l))
1561  fp(3) = sum(p(j_loc,2)*mesh%gauss%ww(:,l))/ray*mode
1562  ftensor(3) = -mode/ray*tensor_loc(2,4,l) - tensor_loc(1,3,l)/ray
1563  !--------calcul de la quatrieme composante : u0(4,:) <--> f(th,m,s)
1564  fs(4) = sum(ff(j_loc,4) * mesh%gauss%ww(:,l))
1565  ft(4) = sum(v1m(j_loc,4) * mesh%gauss%ww(:,l))
1566  fp(4) = -sum(p(j_loc,1)*mesh%gauss%ww(:,l))/ray *mode
1567  ftensor(4) = mode/ray*tensor_loc(2,3,l) - tensor_loc(1,4,l)/ray
1568  !--------calcul de la cinquieme composante : u0(5,:) <--> f(z,m,c)
1569  fs(5) = sum(ff(j_loc,5) * mesh%gauss%ww(:,l))
1570  ft(5) = sum(v1m(j_loc,5) * mesh%gauss%ww(:,l))
1571  fp(5) = sum(p(j_loc,1)*dw_loc(2,:))
1572  ftensor(5) = -mode/ray*tensor_loc(3,4,l)
1573  !--------calcul de la sixieme composante : u0(6,:) <--> f(z,m,s)
1574  fs(6) = sum(ff(j_loc,6) * mesh%gauss%ww(:,l))
1575  ft(6) = sum(v1m(j_loc,6) * mesh%gauss%ww(:,l))
1576  fp(6) = sum(p(j_loc,2)*dw_loc(2,:))
1577  ftensor(6) = mode/ray*tensor_loc(3,3,l)
1578  !-------calcul du second membre pour le terme nonlineaire------------------------
1579  !-------------------------------------------------------------------------------
1580 
1581  smb = (ft+fp-fs+ftensor-rotb_b(index,:))*ray*mesh%gauss%rj(l,m)
1582 
1583  ki = 1
1584  DO ni = 1, nw
1585  ix = (ki-1)*nw+ni
1586  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(1) + &
1587  (dw_loc(1,ni)*tensor_loc(1,1,l) + dw_loc(2,ni)*tensor_loc(1,5,l))*ray*mesh%gauss%rj(l,m)
1588  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(2) + &
1589  (dw_loc(1,ni)*tensor_loc(1,2,l) + dw_loc(2,ni)*tensor_loc(1,6,l))*ray*mesh%gauss%rj(l,m)
1590  v5_loc(ix) = v5_loc(ix) + mesh%gauss%ww(ni,l)*smb(5) + &
1591  (dw_loc(1,ni)*tensor_loc(3,1,l) + dw_loc(2,ni)*tensor_loc(3,5,l))*ray*mesh%gauss%rj(l,m)
1592  v6_loc(ix) = v6_loc(ix) + mesh%gauss%ww(ni,l)*smb(6) + &
1593  (dw_loc(1,ni)*tensor_loc(3,2,l) + dw_loc(2,ni)*tensor_loc(3,6,l))*ray*mesh%gauss%rj(l,m)
1594  END DO
1595 
1596  ki = 2
1597  DO ni = 1, nw
1598  ix = (ki-1)*nw+ni
1599  v14_loc(ix) = v14_loc(ix) + mesh%gauss%ww(ni,l)*smb(4) + &
1600  (dw_loc(1,ni)*tensor_loc(2,2,l) + dw_loc(2,ni)*tensor_loc(2,6,l))*ray*mesh%gauss%rj(l,m)
1601  v23_loc(ix) = v23_loc(ix) + mesh%gauss%ww(ni,l)*smb(3) + &
1602  (dw_loc(1,ni)*tensor_loc(2,1,l) + dw_loc(2,ni)*tensor_loc(2,5,l))*ray*mesh%gauss%rj(l,m)
1603  END DO
1604  ENDDO
1605 
1606  CALL vecsetvalues(vb_3_145, 2*nw, idxm_2, v14_loc, add_values, ierr)
1607  CALL vecsetvalues(vb_3_236, 2*nw, idxm_2, v23_loc, add_values, ierr)
1608  CALL vecsetvalues(vb_3_145, nw, idxm_1, v5_loc, add_values, ierr)
1609  CALL vecsetvalues(vb_3_236, nw, idxm_1, v6_loc, add_values, ierr)
1610  ENDDO
1611 
1612  CALL vecassemblybegin(vb_3_145,ierr)
1613  CALL vecassemblyend(vb_3_145,ierr)
1614  CALL vecassemblybegin(vb_3_236,ierr)
1615  CALL vecassemblyend(vb_3_236,ierr)
1616 
1617  END SUBROUTINE qs_ns_mom_compute_residual_les_3x3
1618 
1619 END MODULE fem_rhs_axi
fem_rhs_axi::qs_00_gauss
subroutine qs_00_gauss(mesh, LA, heat_capa, ff, ff_gauss, vect)
Definition: fem_rhs_axi.f90:775
fem_rhs_axi
Definition: fem_rhs_axi.f90:1
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:3108
fem_rhs_axi::qs_ns_stab_3x3
subroutine qs_ns_stab_3x3(mesh, vv_3_LA, mode, ff, vel_tot, V1m, vit, P, dudt, phalf, nlhalf, dt, vb_145, vb_236, rotv_v, vel_tot_max)
Definition: fem_rhs_axi.f90:358
def_type_mesh::petsc_csr_la
Definition: def_type_mesh.f90:18
def_type_mesh::mesh_type
Definition: def_type_mesh.f90:73
fem_rhs_axi::qs_ns_momentum_stab_3x3
subroutine qs_ns_momentum_stab_3x3(mesh, vv_3_LA, mode, ff, V1m, P, vb_3_145, vb_3_236, rotb_b, tensor, tensor_surface_gauss, stab, momentum, visc_grad_vel)
Definition: fem_rhs_axi.f90:1071
rhs_para_assembling::rhs_3x3
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
Definition: assembling_rhs.f90:5
f
section doc_intro_frame_work_num_app Numerical approximation subsection doc_intro_fram_work_num_app_Fourier_FEM Fourier Finite element representation The SFEMaNS code uses a hybrid Fourier Finite element formulation The Fourier decomposition allows to approximate the problem’s solutions for each Fourier mode modulo nonlinear terms that are made explicit The variables are then approximated on a meridian section of the domain with a finite element method The numerical approximation of a function f $f f is written in the following generic f with f(r,\theta, z)\f $the cylindrical coordinates
def_type_mesh
Definition: def_type_mesh.f90:4
my_util
Definition: my_util.f90:1
fem_rhs_axi::qs_ns_mom_compute_residual_les
subroutine qs_ns_mom_compute_residual_les(mesh, vv_1_LA, vv_2_LA, mode, ff, V1m, P, rotb_b, tensor, tensor_gauss, vb_2_14, vb_2_23, vb_1_5, vb_1_6)
Definition: fem_rhs_axi.f90:1323
input_data
Definition: read_my_data.f90:200
fem_rhs_axi::qs_ns_momentum_stab_new
subroutine qs_ns_momentum_stab_new(mesh, vv_1_LA, vv_2_LA, mode, ff, V1m, P, vb_2_14, vb_2_23, vb_1_5, vb_1_6, rotb_b, tensor, tensor_surface_gauss, stab, momentum, momentum_LES, visc_grad_vel, visc_entro)
Definition: fem_rhs_axi.f90:824
chaine_caractere
Definition: chaine_caractere_par.f90:4
fem_s_axi::qs_00
subroutine qs_00(mesh, ff, u0)
Definition: fem_s_direct_axi.f90:127
fem_rhs_axi::qs_ns_mom_compute_residual_les_3x3
subroutine qs_ns_mom_compute_residual_les_3x3(mesh, vv_3_LA, mode, ff, V1m, P, rotb_b, tensor, tensor_gauss, vb_3_145, vb_3_236)
Definition: fem_rhs_axi.f90:1474
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:1100
rhs_para_assembling
Definition: assembling_rhs.f90:1
my_util::error_petsc
subroutine error_petsc(string)
Definition: my_util.f90:15