16 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab
17 INTEGER,
INTENT(IN) :: type_op, mode
19 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
20 INTEGER,
DIMENSION(:),
ALLOCATABLE :: idxm, idxn
21 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
22 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc, b_loc
23 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: mat_loc
24 REAL(KIND=8) :: ray, eps1, eps2
25 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, k_max, n_w, ix, jx, ki, kj
26 REAL(KIND=8) :: viscolm, xij
27 #include "petsc/finclude/petsc.h"
29 petscerrorcode :: ierr
30 CALL matzeroentries(matrix, ierr)
31 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
34 IF (type_op == 1)
THEN
38 ELSEIF (type_op == 2)
THEN
42 ELSEIF (type_op == 3)
THEN
50 CALL
error_petsc(
'BUG in qs_diff_mass_vect_M, type_op not correct')
53 IF (k_max/=
SIZE(la%loc_to_glob,1))
THEN
54 CALL
error_petsc(
'BUG in qs_diff_mass_vect_petsc_M, k_max/=SIZE(LA%loc_to_glob,1)')
58 DO l = 1, mesh%gauss%l_G
61 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
66 ALLOCATE(mat_loc(k_max*n_w,k_max*n_w),idxm(k_max*n_w),idxn(k_max*n_w))
72 DO l = 1, mesh%gauss%l_G
75 DO ni = 1, n_w; i = jj_loc(ni)
76 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
78 viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
83 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
86 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
87 + viscolm*eps1*wwprod(ni,nj,l)/ray &
88 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
89 + viscolm*mode**2*wwprod(ni,nj,l)/ray
91 b_loc(ni,nj) = b_loc(ni,nj) + eps2*viscolm*2*mode*wwprod(ni,nj,l)/ray
100 iglob = la%loc_to_glob(ki,i)
106 jglob = la%loc_to_glob(kj,j)
110 mat_loc(ix,jx) = mat_loc(ix,jx) + a_loc(ni,nj)
112 mat_loc(ix,jx) = mat_loc(ix,jx) + b_loc(ni,nj)
118 CALL matsetvalues(matrix, k_max*n_w, idxm, k_max*n_w, idxn, mat_loc, add_values, ierr)
120 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
121 CALL matassemblyend(matrix,mat_final_assembly,ierr)
123 DEALLOCATE(mat_loc,idxm,idxn)
132 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab
133 INTEGER,
INTENT(IN) :: mode
134 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
135 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
136 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
137 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
139 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
140 REAL(KIND=8) :: viscolm, xij
141 #include "petsc/finclude/petsc.h"
143 petscerrorcode :: ierr
144 CALL matzeroentries(matrix, ierr)
145 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
147 DO l = 1, mesh%gauss%l_G
148 DO ni = 1, mesh%gauss%n_w
149 DO nj = 1, mesh%gauss%n_w
150 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
156 DO m = 1, mesh%dom_me
157 jj_loc = mesh%jj(:,m)
160 DO nj = 1, mesh%gauss%n_w;
162 jglob = la%loc_to_glob(1,j)
164 DO ni = 1, mesh%gauss%n_w;
166 iglob = la%loc_to_glob(1,i)
171 DO l = 1, mesh%gauss%l_G
174 DO ni = 1, n_w; i = jj_loc(ni)
175 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
177 viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
182 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
184 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
185 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
186 + viscolm*mode**2*wwprod(ni,nj,l)/ray
191 CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
193 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
194 CALL matassemblyend(matrix,mat_final_assembly,ierr)
203 REAL(KIND=8),
INTENT(IN) :: mass, stab
204 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: heat_capa, visco
205 INTEGER,
INTENT(IN) :: mode
206 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
207 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
208 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
209 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
211 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
212 REAL(KIND=8) :: viscolm, xij
213 #include "petsc/finclude/petsc.h"
215 petscerrorcode :: ierr
216 CALL matzeroentries(matrix, ierr)
217 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
219 DO l = 1, mesh%gauss%l_G
220 DO ni = 1, mesh%gauss%n_w
221 DO nj = 1, mesh%gauss%n_w
222 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
228 DO m = 1, mesh%dom_me
229 jj_loc = mesh%jj(:,m)
232 DO nj = 1, mesh%gauss%n_w;
234 jglob = la%loc_to_glob(1,j)
236 DO ni = 1, mesh%gauss%n_w;
238 iglob = la%loc_to_glob(1,i)
243 DO l = 1, mesh%gauss%l_G
246 DO ni = 1, n_w; i = jj_loc(ni)
247 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
249 viscolm = (visco(m) + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
253 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
255 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
256 + heat_capa(m) * mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
257 + viscolm*mode**2*wwprod(ni,nj,l)/ray
262 CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
264 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
265 CALL matassemblyend(matrix,mat_final_assembly,ierr)
282 REAL(KIND=8),
INTENT(IN) :: mass
285 INTEGER ::l, m, ni, nj, i, j, ki, kj, n_w
286 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
287 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij
289 REAL(KIND=8),
DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
290 INTEGER,
DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
291 INTEGER :: ix, jx, iglob, jglob
292 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
293 #include "petsc/finclude/petsc.h"
295 petscerrorcode :: ierr
297 CALL matzeroentries(matrix, ierr)
298 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
303 DO l = 1, mesh%gauss%l_G
304 DO ni = 1, mesh%gauss%n_w
305 DO nj = 1, mesh%gauss%n_w
306 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
312 jj_loc = mesh%jj(:,m)
315 DO l = 1, mesh%gauss%l_G
317 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
318 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
320 DO nj = 1, mesh%gauss%n_w; j = jj_loc(nj)
321 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
322 aij(ni,nj) = aij(ni,nj) + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m)
331 iglob = la%loc_to_glob(ki,i)
337 jglob = la%loc_to_glob(kj,j)
341 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
347 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
350 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
351 CALL matassemblyend(matrix,mat_final_assembly,ierr)
368 INTEGER ,
INTENT(IN) :: type_op, mode
369 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab, stab_bdy_ns
372 INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
373 REAL(KIND=8) :: xij, viscolm
374 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
375 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
376 REAL(KIND=8) :: ray, eps1, eps2,
z, hm1, y, x, stab_normal
377 REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
378 REAL(KIND=8),
DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
379 INTEGER,
DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
380 REAL(KIND=8),
DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
381 INTEGER,
DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
382 INTEGER :: ix, jx, iglob, jglob
383 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
384 REAL(KIND=8) :: hm, viscomode, hh
386 #include "petsc/finclude/petsc.h"
388 petscerrorcode :: ierr
390 CALL matzeroentries(matrix, ierr)
391 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
395 n_ws = mesh%gauss%n_ws
397 IF (type_op == 1)
THEN
401 ELSEIF (type_op == 2)
THEN
405 ELSEIF (type_op == 3)
THEN
415 DO l = 1, mesh%gauss%l_G
416 DO ni = 1, mesh%gauss%n_w
417 DO nj = 1, mesh%gauss%n_w
418 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
424 jj_loc = mesh%jj(:,m)
429 DO l = 1, mesh%gauss%l_G
433 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
434 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
438 hm=0.5d0/inputs%m_max
439 viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m)
440 viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
442 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
444 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
448 DO k = 1, mesh%gauss%k_d
449 xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
453 z = ray * viscolm* xij &
454 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
455 + viscomode*mode**2*wwprod(ni,nj,l)/ray
456 cij(ni,nj) = cij(ni,nj) +
z
457 aij(ni,nj) = aij(ni,nj) +
z + viscomode*eps1*wwprod(ni,nj,l)/ray
459 bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
469 iglob = la%loc_to_glob(ki,i)
475 jglob = la%loc_to_glob(kj,j)
478 IF ((ki .LT. 3) .AND. (kj .LT. 3))
THEN
480 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
482 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
486 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
493 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
502 DO l = 1, mesh%gauss%l_G
506 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
507 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
510 viscolm = visco*mesh%gauss%rj(l,m)*ray
512 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
513 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
514 aij(ni,nj) = aij(ni,nj) &
515 + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2)
516 bij(ni,nj) = bij(ni,nj) &
517 + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2)
518 cij(ni,nj) = cij(ni,nj) &
519 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m))
520 dij(ni,nj) = dij(ni,nj) &
521 + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
522 -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray)
523 eij(ni,nj) = eij(ni,nj) &
524 + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray)
525 fij(ni,nj) = fij(ni,nj) &
526 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
537 iglob = la%loc_to_glob(ki,i)
550 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
553 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
556 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
562 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
565 mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
568 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
574 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
577 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
580 mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
583 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
587 IF (inputs%vv_nb_dirichlet_normal_velocity>0)
THEN
589 stab_normal = coeff_stab_normal*(1.d0+visco)
593 IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
598 DO ls = 1, mesh%gauss%l_Gs
600 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
601 IF (ray.LT.1.d-10) cycle
602 hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
603 x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
604 z = two*mesh%gauss%rjs(ls,ms)*ray*visco
606 DO ni = 1, mesh%gauss%n_ws
607 DO nj = 1, mesh%gauss%n_ws
608 y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
609 aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
610 -
z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
611 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
612 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
613 bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
614 -
z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
615 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
616 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
617 cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
618 -
z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
619 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
620 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
621 dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
622 -
z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
623 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
624 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
638 iglob = la%loc_to_glob(2*ki-1,i)
644 jglob = la%loc_to_glob(2*kj-1,j)
647 IF ((ki == 1) .AND. (kj == 1))
THEN
648 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
649 ELSE IF ((ki == 1) .AND. (kj==2))
THEN
650 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
651 ELSE IF ((ki == 2) .AND. (kj==1))
THEN
652 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
654 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
660 CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
664 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
665 CALL matassemblyend(matrix,mat_final_assembly,ierr)
683 INTEGER ,
INTENT(IN) :: type_op, mode
684 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab, stab_bdy_ns
692 INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
693 REAL(KIND=8) :: xij, viscolm, div_penal
694 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
695 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
696 REAL(KIND=8) :: ray, eps1, eps2,
z, hm1, y, x, stab_normal
697 REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
698 REAL(KIND=8),
DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
699 INTEGER,
DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
700 REAL(KIND=8),
DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
701 INTEGER,
DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
702 INTEGER :: ix, jx, iglob, jglob
703 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
704 REAL(KIND=8) :: viscomode, hm
706 #include "petsc/finclude/petsc.h"
708 petscerrorcode :: ierr
710 CALL matzeroentries(matrix, ierr)
711 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
715 n_ws = mesh%gauss%n_ws
717 IF (type_op == 1)
THEN
721 ELSEIF (type_op == 2)
THEN
725 ELSEIF (type_op == 3)
THEN
735 DO l = 1, mesh%gauss%l_G
736 DO ni = 1, mesh%gauss%n_w
737 DO nj = 1, mesh%gauss%n_w
738 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
744 jj_loc = mesh%jj(:,m)
749 DO l = 1, mesh%gauss%l_G
753 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
754 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
757 viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
758 hm=0.5d0/inputs%m_max
759 viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
761 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
763 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
767 DO k = 1, mesh%gauss%k_d
768 xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
772 z = ray * viscolm* xij &
773 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
774 + viscomode*mode**2*wwprod(ni,nj,l)/ray
775 cij(ni,nj) = cij(ni,nj) +
z
776 aij(ni,nj) = aij(ni,nj) +
z + viscomode*eps1*wwprod(ni,nj,l)/ray
778 bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
833 iglob = la%loc_to_glob(ki,i)
839 jglob = la%loc_to_glob(kj,j)
842 IF ((ki .LT. 3) .AND. (kj .LT. 3))
THEN
844 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
846 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
850 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
857 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
868 DO l = 1, mesh%gauss%l_G
872 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
873 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
876 viscolm = visco*mesh%gauss%rj(l,m)*ray
877 div_penal = inputs%div_stab_in_ns*mesh%gauss%rj(l,m)*ray
879 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
880 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
881 aij(ni,nj) = aij(ni,nj) &
882 + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2) &
883 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(mesh%gauss%dw(1,nj,l,m) &
884 + mesh%gauss%ww(nj,l)/ray)
885 bij(ni,nj) = bij(ni,nj) &
886 + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2) &
887 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(eps2*(mode/ray)*mesh%gauss%ww(nj,l))
888 cij(ni,nj) = cij(ni,nj) &
889 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m)) &
890 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*mesh%gauss%dw(2,nj,l,m)
891 dij(ni,nj) = dij(ni,nj) &
892 + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
893 -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray) &
894 + div_penal*wwprod(ni,nj,l)*(mode/ray)**2
895 eij(ni,nj) = eij(ni,nj) &
896 + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray) &
897 + div_penal*eps2*(mode/ray)*mesh%gauss%ww(ni,l)*mesh%gauss%dw(2,nj,l,m)
898 fij(ni,nj) = fij(ni,nj) &
899 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m)) &
900 + div_penal*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
917 iglob = la%loc_to_glob(ki,i)
930 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
933 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
936 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
942 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
945 mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
948 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
954 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
957 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
960 mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
963 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1039 IF (inputs%vv_nb_dirichlet_normal_velocity>0)
THEN
1041 stab_normal = coeff_stab_normal*(1.d0+visco)
1043 IF (minval(abs(mesh%sides(ms)- inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
1048 DO ls = 1, mesh%gauss%l_Gs
1050 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
1051 IF (ray.LT.1.d-10) cycle
1052 hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
1053 x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
1054 z = two*mesh%gauss%rjs(ls,ms)*ray*visco
1056 DO ni = 1, mesh%gauss%n_ws
1057 DO nj = 1, mesh%gauss%n_ws
1058 y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
1059 aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
1060 -
z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1061 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1062 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1063 bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
1064 -
z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1065 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1066 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1067 cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
1068 -
z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1069 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1070 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1071 dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
1072 -
z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1073 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1074 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1088 iglob = la%loc_to_glob(2*ki-1,i)
1090 idxn_s(ix) = iglob-1
1094 jglob = la%loc_to_glob(2*kj-1,j)
1096 jdxn_s(jx) = jglob-1
1097 IF ((ki == 1) .AND. (kj == 1))
THEN
1098 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
1099 ELSE IF ((ki == 1) .AND. (kj==2))
THEN
1100 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
1101 ELSE IF ((ki == 2) .AND. (kj==1))
THEN
1102 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
1104 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
1110 CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
1143 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1144 CALL matassemblyend(matrix,mat_final_assembly,ierr)
1156 REAL(KIND=8),
INTENT(IN) :: alpha
1158 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
1159 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1160 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1161 INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1162 REAL(KIND=8) :: al, x, ray
1163 #include "petsc/finclude/petsc.h"
1165 petscerrorcode :: ierr
1166 CALL matzeroentries(matrix, ierr)
1168 DO m = 1, mesh%dom_me
1169 jj_loc = mesh%jj(:,m)
1171 DO l = 1, mesh%gauss%l_G
1174 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1175 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1178 al = alpha *mesh%gauss%rj(l,m)*ray
1179 DO nj = 1, mesh%gauss%n_w;
1181 jglob = la%loc_to_glob(1,j)
1183 DO ni = 1, mesh%gauss%n_w;
1185 iglob = la%loc_to_glob(1,i)
1187 x = al*mesh%gauss%ww(nj,l)*mesh%gauss%ww(ni,l)
1188 mat_loc(ni,nj) = mat_loc(ni,nj) + x
1192 CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1194 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1195 CALL matassemblyend(matrix,mat_final_assembly,ierr)
1203 REAL(KIND=8),
INTENT(IN) :: visco
1205 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1206 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1207 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
1208 INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1209 REAL(KIND=8) :: al, x, ray
1210 #include "petsc/finclude/petsc.h"
1212 petscerrorcode :: ierr
1213 CALL matzeroentries(matrix, ierr)
1215 DO m = 1, mesh%dom_me
1216 jj_loc = mesh%jj(:,m)
1219 DO nj = 1, mesh%gauss%n_w;
1221 jglob = la%loc_to_glob(1,j)
1223 DO ni = 1, mesh%gauss%n_w;
1225 iglob = la%loc_to_glob(1,i)
1230 DO l = 1, mesh%gauss%l_G
1233 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1234 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1236 al = visco * mesh%gauss%rj(l,m)
1237 DO nj = 1, mesh%gauss%n_w;
1238 DO ni = 1, mesh%gauss%n_w;
1239 x =al*sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1240 mat_loc(ni,nj) = mat_loc(ni,nj) + x
1244 CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1246 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1247 CALL matassemblyend(matrix,mat_final_assembly,ierr)
1256 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab
1257 INTEGER,
INTENT(IN) :: mode
1258 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
1259 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1260 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1261 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
1263 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
1264 REAL(KIND=8) :: viscolm, xij, hm, viscomode, hh
1265 #include "petsc/finclude/petsc.h"
1267 petscerrorcode :: ierr
1268 CALL matzeroentries(matrix, ierr)
1269 CALL matsetoption(matrix, mat_row_oriented, petsc_false, ierr)
1271 DO l = 1, mesh%gauss%l_G
1272 DO ni = 1, mesh%gauss%n_w
1273 DO nj = 1, mesh%gauss%n_w
1274 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
1279 n_w = mesh%gauss%n_w
1280 DO m = 1, mesh%dom_me
1281 jj_loc = mesh%jj(:,m)
1284 DO nj = 1, mesh%gauss%n_w;
1286 jglob = la%loc_to_glob(1,j)
1288 DO ni = 1, mesh%gauss%n_w;
1290 iglob = la%loc_to_glob(1,i)
1295 DO l = 1, mesh%gauss%l_G
1298 DO ni = 1, n_w; i = jj_loc(ni)
1299 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1303 hm=0.5d0/inputs%m_max
1304 viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m)
1305 viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
1309 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1311 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
1312 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
1313 + viscomode*mode**2*wwprod(ni,nj,l)/ray
1318 CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
1320 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1321 CALL matassemblyend(matrix,mat_final_assembly,ierr)
subroutine qs_diff_mass_vect_m(type_op, LA, mesh, visco, mass, stab, mode, matrix)
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
subroutine qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, stab, mode, matrix)
subroutine qs_11_m(mesh, alpha, ia, ja, a0)
subroutine qs_mass_vect_3x3(LA, mesh, mass, matrix)
subroutine qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)
subroutine qs_diff_mass_scal_m_level(mesh, LA, visco, mass, stab, mode, matrix)
subroutine qs_00_m(mesh, alpha, ia, ja, a0)
subroutine error_petsc(string)
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 z
subroutine qs_diff_mass_vect_3x3_divpenal(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, mode, matrix)