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