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)
16 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotv_v, dudt, nlhalf
17 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: vel_tot
19 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m, vit
20 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p, phalf
21 REAL(KIND=8),
INTENT(IN) :: dt
22 INTEGER,
INTENT(IN) :: mode
24 REAL(KIND=8) :: vel_tot_max
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
54 #include "petsc/finclude/petsc.h"
55 vec :: vb_2_14,vb_2_23,vb_1_5,vb_1_6
56 petscerrorcode :: ierr
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)
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
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
81 IF (mesh%gauss%n_w==3)
THEN
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
91 h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
101 normal_vit = vel_tot_max
104 norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
112 DO m = 1, mesh%dom_me
114 vloc = maxval(vel_tot(j_loc))
115 cfl = max(vloc*dt/mesh%hloc(m),cfl)
118 DO l = 1, mesh%gauss%l_G
120 dw_loc = mesh%gauss%dw(:,:,l,m)
123 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
128 ft(1) = sum(dudt(j_loc,1) *mesh%gauss%ww(:,l))
129 fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
131 ft(2) = sum(dudt(j_loc,2) * mesh%gauss%ww(:,l))
132 fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
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
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
140 ft(5) = sum(dudt(j_loc,5) *mesh%gauss%ww(:,l))
141 fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
143 ft(6) = sum(dudt(j_loc,6) *mesh%gauss%ww(:,l))
144 fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
147 visc1 = max(visc1,abs(ft+fp+nlhalf(index2,:)))
152 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
154 vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
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)
163 visc2(4) = visc2(1); visc2(5) = visc2(1)
164 visc2(3) = visc2(2); visc2(6) = visc2(2)
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
174 viscosity(type,m) = inputs%LES_coeff1*mesh%hloc(m) - visc_eff
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)
183 DO m = 1, mesh%dom_me
185 vloc = maxval(vel_tot(j_loc))
186 cfl = max(vloc*dt/mesh%hloc(m),cfl)
196 DO m = 1, mesh%dom_me
197 mult(:)= viscosity(:,m)
202 iglob = vv_1_la%loc_to_glob(1,i)
209 iglob = vv_2_la%loc_to_glob(ki,i)
219 DO l = 1, mesh%gauss%l_G
221 dw_loc = mesh%gauss%dw(:,:,l,m)
224 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
226 vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
230 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
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
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
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
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
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
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
268 smb = (ft+fp+fs+fv-rotv_v(index,:))*ray*mesh%gauss%rj(l,m)
270 grad(:,type,l) = mult(type)*grad(:,type,l)*ray*mesh%gauss%rj(l,m)
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))
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))
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)
309 IF (mesh%edge_stab)
THEN
310 CALL
error_petsc(
'BUG in qs_ns_stab_new: terms with edge_stab not yet assembled')
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)
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)
370 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotv_v, dudt, nlhalf
371 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: vel_tot
373 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m, vit
374 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: p, phalf
375 REAL(KIND=8),
INTENT(IN) :: dt
376 INTEGER,
INTENT(IN) :: mode
378 REAL(KIND=8) :: vel_tot_max
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
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
410 #include "petsc/finclude/petsc.h"
411 vec :: vb_145, vb_236
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
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
435 IF (mesh%gauss%n_w==3)
THEN
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
445 h(ms) = sum(mesh%gauss%rji(:,ms))/type_fe
455 normal_vit = vel_tot_max
458 norm_vit(type) = sum(abs(vit(:,type)))/mesh%np + 1.d-14
466 DO m = 1, mesh%dom_me
468 vloc = maxval(vel_tot(j_loc))
469 cfl = max(vloc*dt/mesh%hloc(m),cfl)
472 DO l = 1, mesh%gauss%l_G
474 dw_loc = mesh%gauss%dw(:,:,l,m)
477 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
482 ft(1) = sum(dudt(j_loc,1) *mesh%gauss%ww(:,l))
483 fp(1) = sum(phalf(j_loc,1)*dw_loc(1,:))
485 ft(2) = sum(dudt(j_loc,2) * mesh%gauss%ww(:,l))
486 fp(2) = sum(phalf(j_loc,2)*dw_loc(1,:))
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
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
494 ft(5) = sum(dudt(j_loc,5) *mesh%gauss%ww(:,l))
495 fp(5) = sum(phalf(j_loc,1)*dw_loc(2,:))
497 ft(6) = sum(dudt(j_loc,6) *mesh%gauss%ww(:,l))
498 fp(6) = sum(phalf(j_loc,2)*dw_loc(2,:))
501 visc1= max(visc1,abs(ft+fp+nlhalf(index2,:)))
506 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
508 vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
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)
517 visc2(4) = visc2(1); visc2(5) = visc2(1)
518 visc2(3) = visc2(2); visc2(6) = visc2(2)
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
528 viscosity(type,m) = inputs%LES_coeff1*mesh%hloc(m) - visc_eff
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)
537 DO m = 1, mesh%dom_me
539 vloc = maxval(vel_tot(j_loc))
540 cfl = max(vloc*dt/mesh%hloc(m),cfl)
550 DO m = 1, mesh%dom_me
551 mult(:)= viscosity(:,m)
556 iglob = vv_3_la%loc_to_glob(3,i)
563 iglob = vv_3_la%loc_to_glob(ki,i)
573 DO l = 1, mesh%gauss%l_G
575 dw_loc = mesh%gauss%dw(:,:,l,m)
578 grad(k,type,l) = sum(vit(j_loc,type)*dw_loc(k,:))
580 vitloc(type,l) = sum(vit(j_loc,type)*mesh%gauss%ww(:,l))
584 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
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
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
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
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
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
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
622 rhs_gauss(index, :) = (ft+fp+fs+fv-rotv_v(index,:))
628 CALL
rhs_3x3(mesh, vv_3_la, mode, rhs_gauss, vb_145, vb_236)
638 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: gg
639 INTEGER ,
INTENT(IN) :: mode
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
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"
651 petscerrorcode :: ierr
653 CALL vecset(pb_1, 0.d0, ierr)
654 CALL vecset(pb_2, 0.d0, ierr)
656 nw = vv_mesh%gauss%n_w
657 DO l = 1, vv_mesh%gauss%l_G
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))
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)
671 iglob = pp_la%loc_to_glob(1,i)
677 DO l = 1, vv_mesh%gauss%l_G
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)
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
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
698 f =
f *vv_mesh%gauss%rj(l,m)
707 x =
f(1,1)+
f(2,1)+
f(3,1)
709 v1_loc(ni) = v1_loc(ni) + w_c(ni,l)*x
712 x =
f(1,2)+
f(2,2)+
f(3,2)
714 v2_loc(ni) = v2_loc(ni) + w_c(ni,l)*x
718 CALL vecsetvalues(pb_1, nwc, idxm, v1_loc, add_values, ierr)
719 CALL vecsetvalues(pb_2, nwc, idxm, v2_loc, add_values, ierr)
721 CALL vecassemblybegin(pb_1,ierr)
722 CALL vecassemblyend(pb_1,ierr)
724 CALL vecassemblybegin(pb_2,ierr)
725 CALL vecassemblyend(pb_2,ierr)
729 SUBROUTINE qs_00 (mesh, LA, ff, vect)
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"
744 petscerrorcode :: ierr
746 CALL vecset(vect, 0.d0, ierr)
748 DO m = 1, mesh%dom_me
749 jj_loc = mesh%jj(:,m)
751 DO ni = 1, mesh%gauss%n_w
753 iglob = la%loc_to_glob(1,i)
758 DO l = 1, mesh%gauss%l_G
760 DO ni = 1, mesh%gauss%n_w
762 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
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
769 CALL vecsetvalues(vect, mesh%gauss%n_w, idxm, v_loc, add_values, ierr)
771 CALL vecassemblybegin(vect,ierr)
772 CALL vecassemblyend(vect,ierr)
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"
791 petscerrorcode :: ierr
793 CALL vecset(vect, 0.d0, ierr)
796 DO m = 1, mesh%dom_me
797 jj_loc = mesh%jj(:,m)
799 DO ni = 1, mesh%gauss%n_w
801 iglob = la%loc_to_glob(1,i)
806 DO l = 1, mesh%gauss%l_G
809 DO ni = 1, mesh%gauss%n_w
811 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
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
818 CALL vecsetvalues(vect, mesh%gauss%n_w, idxm, v_loc, add_values, ierr)
820 CALL vecassemblybegin(vect,ierr)
821 CALL vecassemblyend(vect,ierr)
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)
836 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotb_b
837 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: tensor
838 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: tensor_surface_gauss
840 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m, momentum, momentum_les
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
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
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
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)
875 IF (.NOT.inputs%LES)
THEN
876 inputs%LES_coeff1=0.d0
879 IF (mesh%gauss%n_w==3)
THEN
887 DO m = 1, mesh%dom_me
888 stab_loc(m) = stab + inputs%LES_coeff1*mesh%hloc(m)
897 DO m = 1, mesh%dom_me
898 mult(:)= - visc_entro(m,1)
904 iglob = vv_1_la%loc_to_glob(1,i)
911 iglob = vv_2_la%loc_to_glob(ki,i)
921 DO l = 1, mesh%gauss%l_G
923 dw_loc = mesh%gauss%dw(:,:,l,m)
926 tensor_loc(k,type,l) = sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
927 + tensor_surface_gauss(k,index,type)
932 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
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
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
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
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
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
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
968 IF (inputs%if_level_set)
THEN
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,:))
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))
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
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
994 fvgm = stab_loc(m)*fvgm + mult*fvgm_les
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
1004 grad_mom(:,type,l) = grad_mom(:,type,l)*ray*mesh%gauss%rj(l,m)
1006 tensor_loc(:,:,l) = tensor_loc(:,:,l) + visc_grad_vel(:,index,:)
1023 smb = (ft+fp+fs+rotb_b(index,:)+fnl+fvgm+fvgu)*ray*mesh%gauss%rj(l,m)
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)
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)
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)
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')
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)
1072 vb_3_145,vb_3_236, rotb_b, tensor, tensor_surface_gauss,&
1073 stab, momentum, visc_grad_vel)
1083 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff, rotb_b
1084 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: tensor
1085 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: tensor_surface_gauss
1087 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m, momentum
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
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
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
1113 CALL vecset(vb_3_145, 0.d0, ierr)
1114 CALL vecset(vb_3_236, 0.d0, ierr)
1119 IF (.NOT.inputs%LES)
THEN
1120 inputs%LES_coeff1=0.d0
1123 IF (mesh%gauss%n_w==3)
THEN
1135 DO m = 1, mesh%dom_me
1137 j_loc = mesh%jj(:,m)
1141 iglob = vv_3_la%loc_to_glob(3,i)
1142 idxm_1(ni) = iglob-1
1148 iglob = vv_3_la%loc_to_glob(ki,i)
1150 idxm_2(ix) = iglob-1
1158 DO l = 1, mesh%gauss%l_G
1160 dw_loc = mesh%gauss%dw(:,:,l,m)
1163 tensor_loc(k,type,l) = sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
1164 + tensor_surface_gauss(k,index,type)
1169 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
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
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
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
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
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
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
1205 IF (inputs%if_level_set)
THEN
1208 grad_mom(k,type,l) = stab_loc(m)*sum(momentum(j_loc,type)*dw_loc(k,:))
1211 momloc(type,l) = sum(momentum(j_loc,type)*mesh%gauss%ww(:,l))
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
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,:))
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
1248 grad_mom(:,:,l) = grad_mom(:,:,l) + stab*grad_t_mom(:,:,l)
1249 fvgm = stab_loc(m)*fvgm + stab*fvgmt
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
1259 grad_mom(:,type,l) = grad_mom(:,type,l)*ray*mesh%gauss%rj(l,m)
1261 tensor_loc(:,:,l) = tensor_loc(:,:,l) + visc_grad_vel(:,index,:)
1278 smb = (ft+fp+fs+rotb_b(index,:)+fnl+fvgm+fvgu)*ray*mesh%gauss%rj(l,m)
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)
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)
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)
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')
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)
1324 tensor,tensor_gauss,vb_2_14,vb_2_23,vb_1_5,vb_1_6)
1332 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1334 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m
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
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
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)
1361 DO m = 1, mesh%dom_me
1362 j_loc = mesh%jj(:,m)
1366 iglob = vv_1_la%loc_to_glob(1,i)
1367 idxm_1(ni) = iglob-1
1373 iglob = vv_2_la%loc_to_glob(ki,i)
1375 idxm_2(ix) = iglob-1
1383 DO l = 1, mesh%gauss%l_G
1385 dw_loc = mesh%gauss%dw(:,:,l,m)
1388 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
1393 tensor_loc(k,type,l) = -sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
1394 + tensor_gauss(k,index,type)
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
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
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
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
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)
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)
1432 smb = (ft+fp-fs+ftensor-rotb_b(index,:))*ray*mesh%gauss%rj(l,m)
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)
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)
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)
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)
1475 tensor,tensor_gauss,vb_3_145,vb_3_236)
1483 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: ff
1485 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: v1m
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
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
1505 CALL vecset(vb_3_145, 0.d0, ierr)
1506 CALL vecset(vb_3_236, 0.d0, ierr)
1510 DO m = 1, mesh%dom_me
1511 j_loc = mesh%jj(:,m)
1515 iglob = vv_3_la%loc_to_glob(3,i)
1516 idxm_1(ni) = iglob-1
1522 iglob = vv_3_la%loc_to_glob(ki,i)
1524 idxm_2(ix) = iglob-1
1532 DO l = 1, mesh%gauss%l_G
1534 dw_loc = mesh%gauss%dw(:,:,l,m)
1537 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
1542 tensor_loc(k,type,l) = -sum(tensor(k,j_loc,type)*mesh%gauss%ww(:,l)) &
1543 + tensor_gauss(k,index,type)
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
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
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
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
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)
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)
1581 smb = (ft+fp-fs+ftensor-rotb_b(index,:))*ray*mesh%gauss%rj(l,m)
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)
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)
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)
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)
subroutine qs_00_gauss(mesh, LA, heat_capa, ff, ff_gauss, vect)
subroutine qs_ns_stab_new(mesh, mode, ff, vel_tot, V1m, vit, P, dudt, phalf, nlhalf, dt, u0, rotv_v)
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)
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)
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
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
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)
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)
subroutine qs_00(mesh, ff, u0)
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)
subroutine qs_01_div_hybrid_2006(uu_mesh, pp_mesh, mode, gg, u0_c)
subroutine error_petsc(string)