8 #include "petsc/finclude/petsc.h" 18 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab
19 INTEGER,
INTENT(IN) :: type_op, mode
21 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
22 INTEGER,
DIMENSION(:),
ALLOCATABLE :: idxm, idxn
23 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
24 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc, b_loc
25 REAL(KIND=8),
DIMENSION(:,:),
ALLOCATABLE :: mat_loc
26 REAL(KIND=8) :: ray, eps1, eps2
27 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, k_max, n_w, ix, jx, ki, kj
28 REAL(KIND=8) :: viscolm, xij
31 petscerrorcode :: ierr
32 CALL matzeroentries (matrix, ierr)
33 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
36 IF (type_op == 1)
THEN 40 ELSEIF (type_op == 2)
THEN 44 ELSEIF (type_op == 3)
THEN 52 CALL error_petsc(
'BUG in qs_diff_mass_vect_M, type_op not correct')
55 IF (k_max/=
SIZE(la%loc_to_glob,1))
THEN 56 CALL error_petsc(
'BUG in qs_diff_mass_vect_petsc_M, k_max/=SIZE(LA%loc_to_glob,1)')
60 DO l = 1, mesh%gauss%l_G
63 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
68 ALLOCATE(mat_loc(k_max*n_w,k_max*n_w),idxm(k_max*n_w),idxn(k_max*n_w))
74 DO l = 1, mesh%gauss%l_G
77 DO ni = 1, n_w; i = jj_loc(ni)
78 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
80 viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
85 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
88 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
89 + viscolm*eps1*wwprod(ni,nj,l)/ray &
90 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
91 + viscolm*mode**2*wwprod(ni,nj,l)/ray
93 b_loc(ni,nj) = b_loc(ni,nj) + eps2*viscolm*2*mode*wwprod(ni,nj,l)/ray
102 iglob = la%loc_to_glob(ki,i)
108 jglob = la%loc_to_glob(kj,j)
112 mat_loc(ix,jx) = mat_loc(ix,jx) + a_loc(ni,nj)
114 mat_loc(ix,jx) = mat_loc(ix,jx) + b_loc(ni,nj)
120 CALL matsetvalues(matrix, k_max*n_w, idxm, k_max*n_w, idxn, mat_loc, add_values, ierr)
122 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
123 CALL matassemblyend(matrix,mat_final_assembly,ierr)
125 DEALLOCATE(mat_loc,idxm,idxn)
134 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab
135 INTEGER,
INTENT(IN) :: mode
136 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
137 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
138 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
139 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
141 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
142 REAL(KIND=8) :: viscolm, xij
145 petscerrorcode :: ierr
146 CALL matzeroentries (matrix, ierr)
147 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
149 DO l = 1, mesh%gauss%l_G
150 DO ni = 1, mesh%gauss%n_w
151 DO nj = 1, mesh%gauss%n_w
152 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
158 DO m = 1, mesh%dom_me
159 jj_loc = mesh%jj(:,m)
162 DO nj = 1, mesh%gauss%n_w;
164 jglob = la%loc_to_glob(1,j)
166 DO ni = 1, mesh%gauss%n_w;
168 iglob = la%loc_to_glob(1,i)
173 DO l = 1, mesh%gauss%l_G
176 DO ni = 1, n_w; i = jj_loc(ni)
177 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
179 viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
184 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
186 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
187 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
188 + viscolm*mode**2*wwprod(ni,nj,l)/ray
193 CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
195 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
196 CALL matassemblyend(matrix,mat_final_assembly,ierr)
201 convection_coeff, stab, mode, matrix)
207 REAL(KIND=8),
INTENT(IN) :: mass, stab
208 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: heat_capa, visco, convection_coeff
209 INTEGER,
DIMENSION(:),
INTENT(IN) :: temp_list_robin_sides
210 INTEGER,
INTENT(IN) :: mode
211 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
212 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
213 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
214 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
216 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w, n_ws, ms, ls, ib
217 REAL(KIND=8) :: viscolm, xij, x
218 REAL(KIND=8),
DIMENSION(mesh%gauss%n_ws,mesh%gauss%n_ws) :: h_phii_phij
219 INTEGER,
DIMENSION(1:1) :: coeff_index
223 petscerrorcode :: ierr
224 CALL matzeroentries (matrix, ierr)
225 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
227 DO l = 1, mesh%gauss%l_G
228 DO ni = 1, mesh%gauss%n_w
229 DO nj = 1, mesh%gauss%n_w
230 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
236 DO m = 1, mesh%dom_me
237 jj_loc = mesh%jj(:,m)
240 DO nj = 1, mesh%gauss%n_w;
242 jglob = la%loc_to_glob(1,j)
244 DO ni = 1, mesh%gauss%n_w;
246 iglob = la%loc_to_glob(1,i)
251 DO l = 1, mesh%gauss%l_G
254 DO ni = 1, n_w; i = jj_loc(ni)
255 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
257 viscolm = (visco(m) + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
261 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
263 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
264 + heat_capa(m) * mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
265 + viscolm*mode**2*wwprod(ni,nj,l)/ray
270 CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
275 IF (
SIZE(temp_list_robin_sides) > 0)
THEN 277 n_ws = mesh%gauss%n_ws
279 IF (minval(abs(temp_list_robin_sides - mesh%sides(ms))) > 0) cycle
281 coeff_index = minloc(abs(temp_list_robin_sides - mesh%sides(ms)))
282 DO ls = 1, mesh%gauss%l_Gs
284 ray = sum(mesh%rr(1,mesh%jjs(:,ms))* mesh%gauss%wws(:,ls))
285 x = convection_coeff(coeff_index(1)) * ray * mesh%gauss%rjs(ls,ms)
288 h_phii_phij(ni,nj) = h_phii_phij(ni,nj) + &
289 x * mesh%gauss%wws(ni,ls) * mesh%gauss%wws(nj,ls)
295 ib = la%loc_to_glob(1,i)
298 CALL matsetvalues(matrix, n_ws, idxn(1:n_ws), n_ws, idxn(1:n_ws), &
299 h_phii_phij, add_values, ierr)
306 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
307 CALL matassemblyend(matrix,mat_final_assembly,ierr)
324 REAL(KIND=8),
INTENT(IN) :: mass
327 INTEGER ::l, m, ni, nj, i, j, ki, kj, n_w
328 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
329 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij
331 REAL(KIND=8),
DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
332 INTEGER,
DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
333 INTEGER :: ix, jx, iglob, jglob
334 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
337 petscerrorcode :: ierr
339 CALL matzeroentries (matrix, ierr)
340 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
345 DO l = 1, mesh%gauss%l_G
346 DO ni = 1, mesh%gauss%n_w
347 DO nj = 1, mesh%gauss%n_w
348 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
354 jj_loc = mesh%jj(:,m)
357 DO l = 1, mesh%gauss%l_G
359 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
360 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
362 DO nj = 1, mesh%gauss%n_w; j = jj_loc(nj)
363 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
364 aij(ni,nj) = aij(ni,nj) + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m)
373 iglob = la%loc_to_glob(ki,i)
379 jglob = la%loc_to_glob(kj,j)
383 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
389 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
392 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
393 CALL matassemblyend(matrix,mat_final_assembly,ierr)
398 SUBROUTINE qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, i_mode, mode, matrix)
410 INTEGER ,
INTENT(IN) :: type_op, mode, i_mode
411 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab, stab_bdy_ns
414 INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
415 REAL(KIND=8) :: xij, viscolm
416 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
417 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
418 REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
419 REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
420 REAL(KIND=8),
DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
421 INTEGER,
DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
422 REAL(KIND=8),
DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
423 INTEGER,
DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
424 INTEGER :: ix, jx, iglob, jglob
425 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
426 REAL(KIND=8) :: hm, viscomode, hh
430 petscerrorcode :: ierr
432 CALL matzeroentries (matrix, ierr)
433 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
437 n_ws = mesh%gauss%n_ws
439 IF (type_op == 1)
THEN 443 ELSEIF (type_op == 2)
THEN 447 ELSEIF (type_op == 3)
THEN 460 DO l = 1, mesh%gauss%l_G
461 DO ni = 1, mesh%gauss%n_w
462 DO nj = 1, mesh%gauss%n_w
463 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
469 jj_loc = mesh%jj(:,m)
474 DO l = 1, mesh%gauss%l_G
478 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
479 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
482 hm=min(mesh%hm(i_mode),hh)
485 viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m)
486 viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
488 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
490 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
494 DO k = 1, mesh%gauss%k_d
495 xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
499 z = ray * viscolm* xij &
500 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
501 + viscomode*mode**2*wwprod(ni,nj,l)/ray
502 cij(ni,nj) = cij(ni,nj) + z
503 aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
505 bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
515 iglob = la%loc_to_glob(ki,i)
521 jglob = la%loc_to_glob(kj,j)
524 IF ((ki .LT. 3) .AND. (kj .LT. 3))
THEN 526 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
528 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
532 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
539 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
548 DO l = 1, mesh%gauss%l_G
552 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
553 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
556 viscolm = visco*mesh%gauss%rj(l,m)*ray
558 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
559 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
560 aij(ni,nj) = aij(ni,nj) &
561 + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2)
562 bij(ni,nj) = bij(ni,nj) &
563 + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2)
564 cij(ni,nj) = cij(ni,nj) &
565 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m))
566 dij(ni,nj) = dij(ni,nj) &
567 + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
568 -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray)
569 eij(ni,nj) = eij(ni,nj) &
570 + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray)
571 fij(ni,nj) = fij(ni,nj) &
572 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
583 iglob = la%loc_to_glob(ki,i)
596 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
599 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
602 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
608 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
611 mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
614 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
620 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
623 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
626 mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
629 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
633 IF (
inputs%vv_nb_dirichlet_normal_velocity>0)
THEN 635 stab_normal = coeff_stab_normal*(1.d0+visco)
639 IF (minval(abs(mesh%sides(ms)-
inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
644 DO ls = 1, mesh%gauss%l_Gs
646 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
647 IF (ray.LT.1.d-10) cycle
648 hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
649 x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
650 z = two*mesh%gauss%rjs(ls,ms)*ray*visco
652 DO ni = 1, mesh%gauss%n_ws
653 DO nj = 1, mesh%gauss%n_ws
654 y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
655 aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
656 - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
657 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
658 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
659 bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
660 - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
661 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
662 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
663 cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
664 - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
665 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
666 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
667 dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
668 - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
669 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
670 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
684 iglob = la%loc_to_glob(2*ki-1,i)
690 jglob = la%loc_to_glob(2*kj-1,j)
693 IF ((ki == 1) .AND. (kj == 1))
THEN 694 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
695 ELSE IF ((ki == 1) .AND. (kj==2))
THEN 696 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
697 ELSE IF ((ki == 2) .AND. (kj==1))
THEN 698 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
700 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
706 CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
710 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
711 CALL matassemblyend(matrix,mat_final_assembly,ierr)
729 INTEGER ,
INTENT(IN) :: type_op, mode, i_mode
730 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab, stab_bdy_ns
738 INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
739 REAL(KIND=8) :: xij, viscolm, div_penal
740 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
741 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
742 REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
743 REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
744 REAL(KIND=8),
DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
745 INTEGER,
DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
746 REAL(KIND=8),
DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
747 INTEGER,
DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
748 INTEGER :: ix, jx, iglob, jglob
749 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
750 REAL(KIND=8) :: viscomode, hm
754 petscerrorcode :: ierr
756 CALL matzeroentries (matrix, ierr)
757 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
761 n_ws = mesh%gauss%n_ws
763 IF (type_op == 1)
THEN 767 ELSEIF (type_op == 2)
THEN 771 ELSEIF (type_op == 3)
THEN 784 DO l = 1, mesh%gauss%l_G
785 DO ni = 1, mesh%gauss%n_w
786 DO nj = 1, mesh%gauss%n_w
787 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
793 jj_loc = mesh%jj(:,m)
798 DO l = 1, mesh%gauss%l_G
802 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
803 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
806 viscolm = (visco + stab*mesh%hloc(m))*mesh%gauss%rj(l,m)
809 hm=min(mesh%hm(i_mode),mesh%hloc(m))
810 viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
812 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
814 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
818 DO k = 1, mesh%gauss%k_d
819 xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
823 z = ray * viscolm* xij &
824 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
825 + viscomode*mode**2*wwprod(ni,nj,l)/ray
826 cij(ni,nj) = cij(ni,nj) + z
827 aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
829 bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
884 iglob = la%loc_to_glob(ki,i)
890 jglob = la%loc_to_glob(kj,j)
893 IF ((ki .LT. 3) .AND. (kj .LT. 3))
THEN 895 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
897 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
901 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
908 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
919 DO l = 1, mesh%gauss%l_G
923 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
924 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
927 viscolm = visco*mesh%gauss%rj(l,m)*ray
930 div_penal =
inputs%div_stab_in_ns/
inputs%Re*mesh%gauss%rj(l,m)*ray
932 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
933 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
934 aij(ni,nj) = aij(ni,nj) &
935 + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2) &
936 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(mesh%gauss%dw(1,nj,l,m) &
937 + mesh%gauss%ww(nj,l)/ray)
938 bij(ni,nj) = bij(ni,nj) &
939 + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2) &
940 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(eps2*(mode/ray)*mesh%gauss%ww(nj,l))
941 cij(ni,nj) = cij(ni,nj) &
942 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m)) &
943 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*mesh%gauss%dw(2,nj,l,m)
944 dij(ni,nj) = dij(ni,nj) &
945 + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
946 -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray) &
947 + div_penal*wwprod(ni,nj,l)*(mode/ray)**2
948 eij(ni,nj) = eij(ni,nj) &
949 + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray) &
950 + div_penal*eps2*(mode/ray)*mesh%gauss%ww(ni,l)*mesh%gauss%dw(2,nj,l,m)
951 fij(ni,nj) = fij(ni,nj) &
952 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m)) &
953 + div_penal*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
970 iglob = la%loc_to_glob(ki,i)
983 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
986 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
989 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
995 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
998 mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
1001 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
1007 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
1010 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
1013 mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
1016 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1092 IF (
inputs%vv_nb_dirichlet_normal_velocity>0)
THEN 1094 stab_normal = coeff_stab_normal*(1.d0+visco)
1096 IF (minval(abs(mesh%sides(ms)-
inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
1101 DO ls = 1, mesh%gauss%l_Gs
1103 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
1104 IF (ray.LT.1.d-10) cycle
1105 hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
1106 x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
1107 z = two*mesh%gauss%rjs(ls,ms)*ray*visco
1109 DO ni = 1, mesh%gauss%n_ws
1110 DO nj = 1, mesh%gauss%n_ws
1111 y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
1112 aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
1113 - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1114 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1115 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1116 bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
1117 - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1118 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1119 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1120 cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
1121 - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1122 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1123 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1124 dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
1125 - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1126 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1127 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1141 iglob = la%loc_to_glob(2*ki-1,i)
1143 idxn_s(ix) = iglob-1
1147 jglob = la%loc_to_glob(2*kj-1,j)
1149 jdxn_s(jx) = jglob-1
1150 IF ((ki == 1) .AND. (kj == 1))
THEN 1151 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
1152 ELSE IF ((ki == 1) .AND. (kj==2))
THEN 1153 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
1154 ELSE IF ((ki == 2) .AND. (kj==1))
THEN 1155 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
1157 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
1163 CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
1196 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1197 CALL matassemblyend(matrix,mat_final_assembly,ierr)
1204 stab_art_comp, i_mode, mode, matrix)
1218 INTEGER ,
INTENT(IN) :: type_op, mode, i_mode
1219 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab, stab_bdy_ns, stab_art_comp
1222 INTEGER :: k, l, m, ni, nj, i, j, np, ki, kj, k_max, ls, ms, n_w, n_ws
1223 REAL(KIND=8) :: xij, viscolm, div_penal
1224 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1225 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: aij, bij, cij, dij, eij, fij
1226 REAL(KIND=8) :: ray, eps1, eps2, z, hm1, y, x, stab_normal
1227 REAL(KIND=8) :: two = 2.d0, coeff_stab_normal=10.d0
1228 REAL(KIND=8),
DIMENSION(3*mesh%gauss%n_w,3*mesh%gauss%n_w) :: mat_loc
1229 INTEGER,
DIMENSION(3*mesh%gauss%n_w) :: idxn, jdxn
1230 REAL(KIND=8),
DIMENSION(2*mesh%gauss%n_ws,2*mesh%gauss%n_ws) :: mat_loc_s
1231 INTEGER,
DIMENSION(2*mesh%gauss%n_ws) :: idxn_s, jdxn_s
1232 INTEGER :: ix, jx, iglob, jglob
1233 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
1234 REAL(KIND=8) :: viscomode, hh, hm
1238 petscerrorcode :: ierr
1240 CALL matzeroentries (matrix, ierr)
1241 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
1243 np =
SIZE(mesh%rr,2)
1244 n_w = mesh%gauss%n_w
1245 n_ws = mesh%gauss%n_ws
1247 IF (type_op == 1)
THEN 1251 ELSEIF (type_op == 2)
THEN 1255 ELSEIF (type_op == 3)
THEN 1268 DO l = 1, mesh%gauss%l_G
1269 DO ni = 1, mesh%gauss%n_w
1270 DO nj = 1, mesh%gauss%n_w
1271 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
1277 jj_loc = mesh%jj(:,m)
1282 DO l = 1, mesh%gauss%l_G
1286 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
1287 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1291 hm=min(mesh%hm(i_mode),hh)
1294 viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m)
1295 viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
1297 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
1299 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
1303 DO k = 1, mesh%gauss%k_d
1304 xij = xij + mesh%gauss%dw(k,nj,l,m) * mesh%gauss%dw(k,ni,l,m)
1308 z = ray * viscolm* xij &
1309 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
1310 + viscomode*mode**2*wwprod(ni,nj,l)/ray
1311 cij(ni,nj) = cij(ni,nj) + z
1312 aij(ni,nj) = aij(ni,nj) + z + viscomode*eps1*wwprod(ni,nj,l)/ray
1314 bij(ni,nj) = bij(ni,nj) + eps2*viscomode*2*mode*wwprod(ni,nj,l)/ray
1324 iglob = la%loc_to_glob(ki,i)
1330 jglob = la%loc_to_glob(kj,j)
1333 IF ((ki .LT. 3) .AND. (kj .LT. 3))
THEN 1335 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
1337 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
1341 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
1348 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1358 DO l = 1, mesh%gauss%l_G
1362 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni,m)
1363 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1366 viscolm = visco*mesh%gauss%rj(l,m)*ray
1367 div_penal = stab_art_comp*mesh%gauss%rj(l,m)*ray
1368 DO nj = 1, mesh%gauss%n_w; j = mesh%jj(nj, m)
1369 DO ni = 1, mesh%gauss%n_w; i = mesh%jj(ni, m)
1370 aij(ni,nj) = aij(ni,nj) &
1371 + viscolm*(mesh%gauss%dw(1,ni,l,m)*mesh%gauss%dw(1,nj,l,m) + wwprod(ni,nj,l)/ray**2) &
1372 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(mesh%gauss%dw(1,nj,l,m) &
1373 + mesh%gauss%ww(nj,l)/ray)
1374 bij(ni,nj) = bij(ni,nj) &
1375 + viscolm*(-eps2*mode*mesh%gauss%ww(ni,l)*mesh%gauss%dw(1,nj,l,m)/ray+eps2*mode*wwprod(ni,nj,l)/ray**2) &
1376 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*(eps2*(mode/ray)*mesh%gauss%ww(nj,l))
1377 cij(ni,nj) = cij(ni,nj) &
1378 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(1,nj,l,m)) &
1379 + div_penal*(mesh%gauss%dw(1,ni,l,m) + mesh%gauss%ww(ni,l)/ray)*mesh%gauss%dw(2,nj,l,m)
1380 dij(ni,nj) = dij(ni,nj) &
1381 + viscolm*(-mesh%gauss%dw(1,ni,l,m)*mesh%gauss%ww(nj,l)/ray+(mode/ray)**2*wwprod(ni,nj,l) &
1382 -mesh%gauss%dw(1,nj,l,m)*mesh%gauss%ww(ni,l)/ray) &
1383 + div_penal*wwprod(ni,nj,l)*(mode/ray)**2
1384 eij(ni,nj) = eij(ni,nj) &
1385 + viscolm*(-eps2*mode*mesh%gauss%dw(2,ni,l,m)*mesh%gauss%ww(nj,l)/ray) &
1386 + div_penal*eps2*(mode/ray)*mesh%gauss%ww(ni,l)*mesh%gauss%dw(2,nj,l,m)
1387 fij(ni,nj) = fij(ni,nj) &
1388 + viscolm*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m)) &
1389 + div_penal*(mesh%gauss%dw(2,ni,l,m)*mesh%gauss%dw(2,nj,l,m))
1400 iglob = la%loc_to_glob(ki,i)
1401 ix = (ki-1)*n_w + ni
1413 mat_loc(ix,jx) = mat_loc(ix,jx) + aij(ni,nj)
1416 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(ni,nj)
1419 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(ni,nj)
1425 mat_loc(ix,jx) = mat_loc(ix,jx) + bij(nj,ni)
1428 mat_loc(ix,jx) = mat_loc(ix,jx) + dij(ni,nj)
1431 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(ni,nj)
1437 mat_loc(ix,jx) = mat_loc(ix,jx) + cij(nj,ni)
1440 mat_loc(ix,jx) = mat_loc(ix,jx) + eij(nj,ni)
1443 mat_loc(ix,jx) = mat_loc(ix,jx) + fij(ni,nj)
1446 CALL matsetvalues(matrix, 3*n_w, idxn, 3*n_w, jdxn, mat_loc, add_values, ierr)
1451 IF (
inputs%vv_nb_dirichlet_normal_velocity>0)
THEN 1453 stab_normal = coeff_stab_normal*(1.d0+visco)
1455 IF (minval(abs(mesh%sides(ms)-
inputs%vv_list_dirichlet_normal_velocity_sides)).NE.0) cycle
1460 DO ls = 1, mesh%gauss%l_Gs
1462 ray = sum(mesh%rr(1,mesh%jjs(:,ms))*mesh%gauss%wws(:,ls))
1463 IF (ray.LT.1.d-10) cycle
1464 hm1 = stab_bdy_ns/sum(mesh%gauss%rjs(:,ms))
1465 x = two*stab_normal*hm1*mesh%gauss%rjs(ls,ms)*ray
1466 z = two*mesh%gauss%rjs(ls,ms)*ray*visco
1468 DO ni = 1, mesh%gauss%n_ws
1469 DO nj = 1, mesh%gauss%n_ws
1470 y = x * mesh%gauss%wws(ni,ls)*mesh%gauss%wws(nj,ls)
1471 aij(ni,nj) = aij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)**2 &
1472 - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1473 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1474 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1475 bij(ni,nj) = bij(ni,nj) + y *mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms) &
1476 - z*mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%wws(ni,ls) &
1477 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1478 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1479 cij(ni,nj) = cij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%rnorms(1,ls,ms) &
1480 - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1481 *(mesh%gauss%rnorms(1,ls,ms)**2*mesh%gauss%dw_s(1,nj,ls,ms) &
1482 + mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(2,nj,ls,ms))
1483 dij(ni,nj) = dij(ni,nj) + y *mesh%gauss%rnorms(2,ls,ms)**2 &
1484 - z*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%wws(ni,ls) &
1485 *(mesh%gauss%rnorms(1,ls,ms)*mesh%gauss%rnorms(2,ls,ms)*mesh%gauss%dw_s(1,nj,ls,ms) &
1486 + mesh%gauss%rnorms(2,ls,ms)**2*mesh%gauss%dw_s(2,nj,ls,ms))
1500 iglob = la%loc_to_glob(2*ki-1,i)
1502 idxn_s(ix) = iglob-1
1506 jglob = la%loc_to_glob(2*kj-1,j)
1508 jdxn_s(jx) = jglob-1
1509 IF ((ki == 1) .AND. (kj == 1))
THEN 1510 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + aij(ni,nj) + aij(nj,ni)
1511 ELSE IF ((ki == 1) .AND. (kj==2))
THEN 1512 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + bij(ni,nj) + cij(nj,ni)
1513 ELSE IF ((ki == 2) .AND. (kj==1))
THEN 1514 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + cij(ni,nj) + bij(nj,ni)
1516 mat_loc_s(ix,jx) = mat_loc_s(ix,jx) + dij(ni,nj) + dij(nj,ni)
1522 CALL matsetvalues(matrix, 2*n_ws, idxn_s, 2*n_ws, jdxn_s, mat_loc_s, add_values, ierr)
1526 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1527 CALL matassemblyend(matrix,mat_final_assembly,ierr)
1532 SUBROUTINE qs_00_m (mesh, alpha, LA, matrix)
1536 REAL(KIND=8),
INTENT(IN) :: alpha
1538 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
1539 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1540 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1541 INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1542 REAL(KIND=8) :: al, x, ray
1545 petscerrorcode :: ierr
1546 CALL matzeroentries (matrix, ierr)
1548 DO m = 1, mesh%dom_me
1549 jj_loc = mesh%jj(:,m)
1551 DO l = 1, mesh%gauss%l_G
1554 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1555 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1558 al = alpha *mesh%gauss%rj(l,m)*ray
1559 DO nj = 1, mesh%gauss%n_w;
1561 jglob = la%loc_to_glob(1,j)
1563 DO ni = 1, mesh%gauss%n_w;
1565 iglob = la%loc_to_glob(1,i)
1567 x = al*mesh%gauss%ww(nj,l)*mesh%gauss%ww(ni,l)
1568 mat_loc(ni,nj) = mat_loc(ni,nj) + x
1572 CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1574 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1575 CALL matassemblyend(matrix,mat_final_assembly,ierr)
1578 SUBROUTINE qs_11_m (mesh, visco, LA, matrix)
1583 REAL(KIND=8),
INTENT(IN) :: visco
1585 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: mat_loc
1586 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1587 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
1588 INTEGER :: m, l, ni, nj, i, j, iglob, jglob
1589 REAL(KIND=8) :: al, x, ray
1592 petscerrorcode :: ierr
1593 CALL matzeroentries (matrix, ierr)
1595 DO m = 1, mesh%dom_me
1596 jj_loc = mesh%jj(:,m)
1599 DO nj = 1, mesh%gauss%n_w;
1601 jglob = la%loc_to_glob(1,j)
1603 DO ni = 1, mesh%gauss%n_w;
1605 iglob = la%loc_to_glob(1,i)
1610 DO l = 1, mesh%gauss%l_G
1613 DO ni = 1, mesh%gauss%n_w; i = jj_loc(ni)
1614 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1616 al = visco * mesh%gauss%rj(l,m)
1617 DO nj = 1, mesh%gauss%n_w;
1618 DO ni = 1, mesh%gauss%n_w;
1619 x =al*sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1620 mat_loc(ni,nj) = mat_loc(ni,nj) + x
1624 CALL matsetvalues(matrix, mesh%gauss%n_w, idxm, mesh%gauss%n_w, idxn, mat_loc, add_values, ierr)
1626 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1627 CALL matassemblyend(matrix,mat_final_assembly,ierr)
1636 REAL(KIND=8),
INTENT(IN) :: visco, mass, stab
1637 INTEGER,
INTENT(IN) :: mode, i_mode
1638 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
1639 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
1640 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
1641 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
1643 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
1644 REAL(KIND=8) :: viscolm, xij, hm, viscomode, hh
1647 petscerrorcode :: ierr
1648 CALL matzeroentries (matrix, ierr)
1649 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
1651 DO l = 1, mesh%gauss%l_G
1652 DO ni = 1, mesh%gauss%n_w
1653 DO nj = 1, mesh%gauss%n_w
1654 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
1659 n_w = mesh%gauss%n_w
1660 DO m = 1, mesh%dom_me
1661 jj_loc = mesh%jj(:,m)
1664 DO nj = 1, mesh%gauss%n_w;
1666 jglob = la%loc_to_glob(1,j)
1668 DO ni = 1, mesh%gauss%n_w;
1670 iglob = la%loc_to_glob(1,i)
1675 DO l = 1, mesh%gauss%l_G
1678 DO ni = 1, n_w; i = jj_loc(ni)
1679 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
1682 hm=min(mesh%hm(i_mode),hh)
1685 viscolm = (visco + stab*hh)*mesh%gauss%rj(l,m)
1686 viscomode = (visco + stab*hm)*mesh%gauss%rj(l,m)
1690 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
1692 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
1693 + mass*ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
1694 + viscomode*mode**2*wwprod(ni,nj,l)/ray
1699 CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
1701 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
1702 CALL matassemblyend(matrix,mat_final_assembly,ierr)
subroutine qs_11_m(mesh, visco, LA, matrix)
subroutine qs_diff_mass_scal_m_variant(mesh, LA, heat_capa, visco, mass, temp_list_robin_sides, convection_coeff, stab, mode, matrix)
subroutine error_petsc(string)
subroutine qs_diff_mass_scal_m(mesh, LA, visco, mass, stab, mode, matrix)
subroutine qs_diff_mass_vect_m(type_op, LA, mesh, visco, mass, stab, mode, matrix)
subroutine qs_diff_mass_vect_3x3_divpenal(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, i_mode, mode, matrix)
subroutine qs_mass_vect_3x3(LA, mesh, mass, matrix)
subroutine qs_00_m(mesh, alpha, LA, matrix)
subroutine qs_diff_mass_vect_3x3(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, i_mode, mode, matrix)
subroutine qs_diff_mass_vect_3x3_divpenal_art_comp(type_op, LA, mesh, visco, mass, stab, stab_bdy_ns, stab_art_comp, i_mode, mode, matrix)
subroutine qs_diff_mass_scal_m_level(mesh, LA, visco, mass, stab, i_mode, mode, matrix)