7 #include "petsc/finclude/petsc.h" 19 chmp_vit, max_vel, my_par_cc, cc_list_dirichlet_sides, cc_per, nb_inter, &
38 REAL(KIND=8) :: time, dt
39 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
40 INTEGER,
INTENT(IN) :: nb_inter
45 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: cn_m1, cn
46 INTEGER,
DIMENSION(:),
INTENT(IN) :: cc_list_dirichlet_sides
47 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: chmp_vit
48 REAL(KIND=8),
INTENT(INOUT) :: max_vel
49 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: visc_entro_level
50 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: cc_global_D
51 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: cc_mode_global_js_D
52 LOGICAL,
SAVE :: once = .true., once_vel=.true.
53 LOGICAL,
SAVE :: re_init=.false.
54 INTEGER,
SAVE :: m_max_c
55 INTEGER,
DIMENSION(:),
POINTER,
SAVE :: cc_js_D
56 INTEGER,
SAVE :: my_petscworld_rank, nb_procs
57 REAL(KIND=8),
SAVE :: LES_coeff1_in_level
61 INTEGER,
POINTER,
DIMENSION(:) :: cc_1_ifrom
64 INTEGER :: bloc_size, m_max_pad
66 REAL(KIND=8),
DIMENSION(cc_mesh%np) :: ff
67 REAL(KIND=8),
DIMENSION(cc_mesh%np, 2) :: cn_p1
68 REAL(KIND=8),
DIMENSION(cc_mesh%np,2,SIZE(list_mode)) :: cext
69 REAL(KIND=8),
DIMENSION(cc_mesh%np,2,SIZE(list_mode)) :: cext_reg
70 REAL(KIND=8),
DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 2, SIZE(list_mode)) :: ff_conv
71 REAL(KIND=8),
DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 2, SIZE(list_mode)) :: ff_comp
72 REAL(KIND=8),
DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 6, SIZE(list_mode)) :: ff_entro
73 REAL(KIND=8),
DIMENSION(cc_mesh%gauss%l_G*cc_mesh%me, 2, SIZE(list_mode)) :: ff_phi_1mphi
74 REAL(KIND=8) :: int_mass_correct
75 REAL(KIND=8) :: tps, tps_tot, tps_cumul
76 REAL(KIND=8) :: one, zero, three
77 DATA zero, one, three/0.d0,1.d0,3.d0/
80 petscerrorcode :: ierr
81 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d
82 mat,
DIMENSION(:),
POINTER,
SAVE :: cc_mat
83 vec,
SAVE :: cb_1, cb_2, cx_1, cx_1_ghost
84 ksp,
DIMENSION(:),
POINTER,
SAVE :: cc_ksp
85 vec,
SAVE :: cb_reg_1, cb_reg_2
86 mat,
DIMENSION(:),
POINTER,
SAVE :: cc_reg_mat
87 ksp,
DIMENSION(:),
POINTER,
SAVE :: cc_reg_ksp
94 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
95 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
100 CALL veccreateghost(comm_one_d(1), n, &
101 petsc_determine,
SIZE(cc_1_ifrom), cc_1_ifrom, cx_1, ierr)
102 CALL vecghostgetlocalform(cx_1, cx_1_ghost, ierr)
103 CALL vecduplicate(cx_1, cb_1, ierr)
104 CALL vecduplicate(cx_1, cb_2, ierr)
105 CALL vecduplicate(cx_1, cb_reg_1, ierr)
106 CALL vecduplicate(cx_1, cb_reg_2, ierr)
110 m_max_c =
SIZE(list_mode)
118 ALLOCATE(cc_global_d(m_max_c))
120 ALLOCATE(cc_global_d(i)%DRL(
SIZE(cc_mode_global_js_d(i)%DIL)))
125 ALLOCATE(cc_mat(m_max_c),cc_ksp(m_max_c))
127 IF (
inputs%if_compression_mthd_JLG)
THEN 128 ALLOCATE(cc_reg_mat(m_max_c),cc_reg_ksp(m_max_c))
134 CALL qs_regul_m (cc_mesh, cc_1_la, 3.d0, i, mode, cc_reg_mat(i))
135 IF (cc_per%n_bord/=0)
THEN 139 CALL init_solver(my_par_cc,cc_reg_ksp(i),cc_reg_mat(i),comm_one_d(1),&
140 solver=my_par_cc%solver,precond=my_par_cc%precond)
145 IF (
inputs%if_level_set_P2)
THEN 146 les_coeff1_in_level=
inputs%LES_coeff1
148 les_coeff1_in_level=4.d0*
inputs%LES_coeff1
159 IF (my_petscworld_rank==0)
WRITE(*,*)
' Recomputing matrix for level set function' 160 IF (my_petscworld_rank==0)
WRITE(*,*)
' NEW MAX VEL test 1', time, max_vel
167 IF (
inputs%if_level_bdf2)
THEN 169 les_coeff1_in_level, i, mode, cc_mat(i))
172 les_coeff1_in_level, i, mode, cc_mat(i))
174 IF (cc_per%n_bord/=0)
THEN 178 CALL init_solver(my_par_cc,cc_ksp(i),cc_mat(i),comm_one_d(1),&
179 solver=my_par_cc%solver,precond=my_par_cc%precond, opt_re_init=re_init)
187 IF (
inputs%if_level_bdf2)
THEN 193 IF (
inputs%if_compression_mthd_JLG)
THEN 197 CALL qs_00 (cc_mesh,cc_1_la, cext(:,1,i), cb_reg_1)
198 CALL qs_00 (cc_mesh,cc_1_la, cext(:,2,i), cb_reg_2)
201 IF (cc_per%n_bord/=0)
THEN 202 CALL periodic_rhs_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cb_reg_1, cc_1_la)
203 CALL periodic_rhs_petsc(cc_per%n_bord, cc_per%list, cc_per%perlist, cb_reg_2, cc_1_la)
208 cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,1,cc_mesh%rr(:,cc_js_d), mode, time)
209 cc_global_d(i)%DRL(n+1:) = 0.d0
210 CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_reg_1)
211 cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,2,cc_mesh%rr(:,cc_js_d), mode, time)
212 cc_global_d(i)%DRL(n+1:) = 0.d0
213 CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_reg_2)
218 CALL solver(cc_reg_ksp(i),cb_reg_1,cx_1,reinit=.false.,
verbose=my_par_cc%verbose)
219 CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
220 CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
221 CALL extract(cx_1_ghost,1,1,cc_1_la,cext_reg(:,1,i))
224 CALL solver(cc_reg_ksp(i),cb_reg_2,cx_1,reinit=.false.,
verbose=my_par_cc%verbose)
225 CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
226 CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
227 CALL extract(cx_1_ghost,1,1,cc_1_la,cext_reg(:,2,i))
228 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
238 IF (
inputs%if_compression_mthd_JLG)
THEN 240 list_mode, cext, cext_reg, visc_entro_level,&
241 les_coeff1_in_level, nb_procs, ff_entro, ff_phi_1mphi)
244 chmp_vit, cext, nb_procs, ff_comp)
245 ff_conv=ff_conv-ff_comp
248 les_coeff1_in_level, nb_procs, ff_entro, ff_phi_1mphi)
251 IF (
inputs%if_mass_correction)
THEN 254 int_mass_correct = 0.d0
262 IF (
inputs%if_level_bdf2)
THEN 263 ff = (2/dt)*cn(:,1,i) - (1/(2*dt))*cn_m1(:,1,i) &
264 +source_in_level_set(nb_inter,1, cc_mesh%rr, mode, time)
266 mode, 1, cb_1, cext(:,1,i), &
267 ff_entro(:,:,i), -ff_phi_1mphi(:,1,i), int_mass_correct)
269 ff = (2/dt)*cn(:,2,i) - (1/(2*dt))*cn_m1(:,2,i) &
270 +source_in_level_set(nb_inter,2, cc_mesh%rr, mode, time)
272 mode, 2, cb_2, cext(:,2,i), &
273 ff_entro(:,:,i), -ff_phi_1mphi(:,2,i), int_mass_correct)
275 ff = (1/dt)*cn(:,1,i) + source_in_level_set(nb_inter,1, cc_mesh%rr, mode, time)
277 mode, 1, cb_1, cext(:,1,i), &
278 ff_entro(:,:,i), -ff_phi_1mphi(:,1,i), int_mass_correct)
280 ff = (1/dt)*cn(:,2,i) + source_in_level_set(nb_inter,2, cc_mesh%rr, mode, time)
282 mode, 2, cb_2, cext(:,2,i), &
283 ff_entro(:,:,i), -ff_phi_1mphi(:,2,i), int_mass_correct)
287 IF (cc_per%n_bord/=0)
THEN 294 cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,1,cc_mesh%rr(:,cc_js_d), mode, time)
295 cc_global_d(i)%DRL(n+1:) = 0.d0
296 CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_1)
297 cc_global_d(i)%DRL(1:n) = level_set_exact(nb_inter,2,cc_mesh%rr(:,cc_js_d), mode, time)
298 cc_global_d(i)%DRL(n+1:) = 0.d0
299 CALL dirichlet_rhs(cc_mode_global_js_d(i)%DIL-1,cc_global_d(i)%DRL,cb_2)
301 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
308 CALL solver(cc_ksp(i),cb_1,cx_1,reinit=.false.,
verbose=my_par_cc%verbose)
309 CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
310 CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
311 CALL extract(cx_1_ghost,1,1,cc_1_la,cn_p1(:,1))
314 CALL solver(cc_ksp(i),cb_2,cx_1,reinit=.false.,
verbose=my_par_cc%verbose)
315 CALL vecghostupdatebegin(cx_1,insert_values,scatter_forward,ierr)
316 CALL vecghostupdateend(cx_1,insert_values,scatter_forward,ierr)
317 CALL extract(cx_1_ghost,1,1,cc_1_la,cn_p1(:,2))
318 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
329 cn_m1(:,:,i) = cn(:,:,i)
332 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
338 bloc_size =
SIZE(cn,1)/nb_procs+1
339 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
364 INTEGER,
INTENT(IN) :: nb_procs
365 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
366 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: V_in, c_in
367 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: c_out
368 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: Gradc, W
369 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: Div, Cgauss
371 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
372 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
373 INTEGER :: m, l , i, mode, index, k
374 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: Vs
375 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: cs
376 INTEGER,
DIMENSION(:,:),
POINTER :: jj
377 INTEGER,
POINTER :: me
378 REAL(KIND=8) :: ray, tps
379 REAL(KIND=8),
DIMENSION(3) :: temps
380 INTEGER :: m_max_pad, bloc_size
382 mpi_comm :: communicator
389 DO i = 1,
SIZE(list_mode)
392 DO m = 1, mesh%dom_me
395 vs(:,k) = v_in(j_loc,k,i)
398 cs(:,k) = c_in(j_loc,k,i)
405 ray = sum(mesh%rr(1,j_loc)*
ww(:,l))
408 w(index,1,i) = sum(vs(:,1)*
ww(:,l))
409 w(index,3,i) = sum(vs(:,3)*
ww(:,l))
410 w(index,5,i) = sum(vs(:,5)*
ww(:,l))
412 w(index,2,i) = sum(vs(:,2)*
ww(:,l))
413 w(index,4,i) = sum(vs(:,4)*
ww(:,l))
414 w(index,6,i) = sum(vs(:,6)*
ww(:,l))
416 div(index,1,i) = sum(vs(:,1)*dw_loc(1,:)) + sum(vs(:,1)*
ww(:,l))/ray &
417 + (mode/ray)*sum(vs(:,4)*
ww(:,l)) + sum(vs(:,5)*dw_loc(2,:))
418 div(index,2,i) = sum(vs(:,2)*dw_loc(1,:)) + sum(vs(:,2)*
ww(:,l))/ray &
419 - (mode/ray)*sum(vs(:,3)*
ww(:,l)) + sum(vs(:,6)*dw_loc(2,:))
423 gradc(index,1,i) = sum(cs(:,1)*dw_loc(1,:))
424 gradc(index,2,i) = sum(cs(:,2)*dw_loc(1,:))
425 gradc(index,3,i) = mode/ray*sum(cs(:,2)*
ww(:,l))
426 gradc(index,4,i) = -mode/ray*sum(cs(:,1)*
ww(:,l))
427 gradc(index,5,i) = sum(cs(:,1)*dw_loc(2,:))
428 gradc(index,6,i) = sum(cs(:,2)*dw_loc(2,:))
430 cgauss(index,1,i) = sum(cs(:,1)*
ww(:,l))
431 cgauss(index,2,i) = sum(cs(:,2)*
ww(:,l))
442 bloc_size =
SIZE(gradc,1)/nb_procs+1
443 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
468 INTEGER,
INTENT(IN) :: nb_procs
469 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
470 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: V_in, c_in
471 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: c_out
472 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: Gradc, W
473 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: Cgauss
474 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
475 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
476 INTEGER :: m, l , i, mode, index, k
477 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: Vs
478 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: cs
479 INTEGER,
DIMENSION(:,:),
POINTER :: jj
480 INTEGER,
POINTER :: me
481 REAL(KIND=8) :: ray, tps, norm_vel_L2, Volume_3D
482 INTEGER :: m_max_pad, bloc_size
484 mpi_comm,
DIMENSION(2) :: communicator
491 DO i = 1,
SIZE(list_mode)
494 DO m = 1, mesh%dom_me
497 vs(:,k) = v_in(j_loc,k,i)
500 cs(:,k) = c_in(j_loc,k,i)
507 ray = sum(mesh%rr(1,j_loc)*
ww(:,l))
510 w(index,1,i) = sum(vs(:,1)*
ww(:,l))
511 w(index,3,i) = sum(vs(:,3)*
ww(:,l))
512 w(index,5,i) = sum(vs(:,5)*
ww(:,l))
514 w(index,2,i) = sum(vs(:,2)*
ww(:,l))
515 w(index,4,i) = sum(vs(:,4)*
ww(:,l))
516 w(index,6,i) = sum(vs(:,6)*
ww(:,l))
520 gradc(index,1,i) = sum(cs(:,1)*dw_loc(1,:))
521 gradc(index,2,i) = sum(cs(:,2)*dw_loc(1,:))
522 gradc(index,3,i) = mode/ray*sum(cs(:,2)*
ww(:,l))
523 gradc(index,4,i) = -mode/ray*sum(cs(:,1)*
ww(:,l))
524 gradc(index,5,i) = sum(cs(:,1)*dw_loc(2,:))
525 gradc(index,6,i) = sum(cs(:,2)*dw_loc(2,:))
527 cgauss(index,1,i) = sum(cs(:,1)*
ww(:,l))
528 cgauss(index,2,i) = sum(cs(:,2)*
ww(:,l))
533 bloc_size = mesh%gauss%l_G*mesh%me/nb_procs+1
534 bloc_size = mesh%gauss%l_G*(bloc_size/mesh%gauss%l_G)+mesh%gauss%l_G
535 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
539 volume_3d = volume_3d*2*acos(-1.d0)
540 norm_vel_l2 =
norm_sf(communicator,
'L2', mesh, list_mode, v_in)/volume_3d
544 mesh%hloc_gauss, mesh%gauss%l_G, nb_procs, bloc_size, m_max_pad)
559 REAL(KIND=8),
INTENT(OUT) :: RESLT
560 REAL(KIND=8) :: vol_loc, vol_out
561 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
562 INTEGER :: m, l , i , ni, code
565 mpi_comm :: communicator
568 DO m = 1, mesh%dom_me
570 DO l = 1, mesh%gauss%l_G
572 DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
573 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
575 vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
578 CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
582 SUBROUTINE qs_regul_m (mesh, LA, stab, i_mode, mode, matrix)
589 REAL(KIND=8),
INTENT(IN) :: stab
590 INTEGER,
INTENT(IN) :: mode, i_mode
591 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
592 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm, idxn
593 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w,mesh%gauss%l_G) :: wwprod
594 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,mesh%gauss%n_w) :: a_loc
596 INTEGER :: m, l, ni, nj, i, j, iglob, jglob, n_w
597 REAL(KIND=8) :: viscolm, xij, viscomode, hm, hh
600 petscerrorcode :: ierr
602 CALL matzeroentries (matrix, ierr)
603 CALL matsetoption (matrix, mat_row_oriented, petsc_false, ierr)
605 DO l = 1, mesh%gauss%l_G
606 DO ni = 1, mesh%gauss%n_w
607 DO nj = 1, mesh%gauss%n_w
608 wwprod(ni,nj,l) = mesh%gauss%ww(ni,l)*mesh%gauss%ww(nj,l)
614 DO m = 1, mesh%dom_me
615 jj_loc = mesh%jj(:,m)
618 DO nj = 1, mesh%gauss%n_w;
620 jglob = la%loc_to_glob(1,j)
622 DO ni = 1, mesh%gauss%n_w;
624 iglob = la%loc_to_glob(1,i)
629 DO l = 1, mesh%gauss%l_G
632 DO ni = 1, n_w; i = jj_loc(ni)
633 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
636 hm=min(mesh%hm(i_mode),hh)
640 viscolm = (stab*hh)**2*mesh%gauss%rj(l,m)
641 viscomode = (stab*hm)**2*mesh%gauss%rj(l,m)
645 xij = sum(mesh%gauss%dw(:,nj,l,m)*mesh%gauss%dw(:,ni,l,m))
647 a_loc(ni,nj) = a_loc(ni,nj) + ray * viscolm* xij &
648 + ray*wwprod(ni,nj,l)*mesh%gauss%rj(l,m) &
649 + viscomode*mode**2*wwprod(ni,nj,l)/ray
654 CALL matsetvalues(matrix,n_w,idxm,n_w,idxn,a_loc,add_values,ierr)
656 CALL matassemblybegin(matrix,mat_final_assembly,ierr)
657 CALL matassemblyend(matrix,mat_final_assembly,ierr)
662 coeff1_in_level, nb_procs, v_out, c_out)
671 INTEGER,
INTENT(IN) :: nb_procs
672 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
673 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: c_in
674 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: c_reg_in
675 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: visc_entro_real
676 REAL(KIND=8),
INTENT(IN) :: coeff1_in_level
677 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)),
INTENT(OUT) :: V_out
678 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)),
INTENT(OUT) :: c_out
679 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: Gradc_in, Gradc_reg_in
680 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)) :: c_in_gauss
681 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
682 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
683 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: c_in_loc
684 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: c_reg_in_loc
685 INTEGER,
DIMENSION(:,:),
POINTER :: jj
686 INTEGER,
POINTER :: me
687 REAL(KIND=8) :: ray, hh, hm
688 INTEGER :: m, l , i, mode, index, k
689 INTEGER :: m_max_pad, bloc_size
691 mpi_comm,
DIMENSION(2) :: communicator
697 DO i = 1,
SIZE(list_mode)
700 DO m = 1, mesh%dom_me
703 c_in_loc(:,k) = c_in(j_loc,k,i)
704 c_reg_in_loc(:,k) = c_reg_in(j_loc,k,i)
712 hh=mesh%hloc_gauss(index)
713 hm=min(mesh%hm(i),hh)
720 ray = sum(mesh%rr(1,j_loc)*
ww(:,l))
723 c_in_gauss(index,1,i) = sum(c_in_loc(:,1)*
ww(:,l))
724 c_in_gauss(index,2,i) = sum(c_in_loc(:,2)*
ww(:,l))
727 gradc_in(index,1,i) = sum(c_in_loc(:,1)*dw_loc(1,:))*hh
728 gradc_in(index,2,i) = sum(c_in_loc(:,2)*dw_loc(1,:))*hh
729 gradc_in(index,3,i) = mode/ray*sum(c_in_loc(:,2)*
ww(:,l))*hm
730 gradc_in(index,4,i) = -mode/ray*sum(c_in_loc(:,1)*
ww(:,l))*hm
731 gradc_in(index,5,i) = sum(c_in_loc(:,1)*dw_loc(2,:))*hh
732 gradc_in(index,6,i) = sum(c_in_loc(:,2)*dw_loc(2,:))*hh
735 gradc_reg_in(index,1,i) = sum(c_reg_in_loc(:,1)*dw_loc(1,:))*hh
736 gradc_reg_in(index,2,i) = sum(c_reg_in_loc(:,2)*dw_loc(1,:))*hh
737 gradc_reg_in(index,3,i) = mode/ray*sum(c_reg_in_loc(:,2)*
ww(:,l))*hm
738 gradc_reg_in(index,4,i) = -mode/ray*sum(c_reg_in_loc(:,1)*
ww(:,l))*hm
739 gradc_reg_in(index,5,i) = sum(c_reg_in_loc(:,1)*dw_loc(2,:))*hh
740 gradc_reg_in(index,6,i) = sum(c_reg_in_loc(:,2)*dw_loc(2,:))*hh
744 bloc_size =
SIZE(c_in_gauss,1)/nb_procs+1
746 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
751 mesh%hloc_gauss, coeff1_in_level, v_out, c_out, nb_procs, bloc_size, m_max_pad)
756 coeff1_in_level, nb_procs, v_out, c_out)
765 INTEGER,
INTENT(IN) :: nb_procs
766 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
767 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: c_in
768 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: visc_entro_real
769 REAL(KIND=8),
INTENT(IN) :: coeff1_in_level
770 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)),
INTENT(OUT) :: V_out
771 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)),
INTENT(OUT) :: c_out
772 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,6,SIZE(list_mode)) :: Gradc_in
773 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%dom_me,2,SIZE(list_mode)) :: c_in_gauss
774 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
775 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
776 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,2) :: c_in_loc
777 INTEGER,
DIMENSION(:,:),
POINTER :: jj
778 INTEGER,
POINTER :: me
779 REAL(KIND=8) :: ray, hh, hm
780 INTEGER :: m, l , i, mode, index, k
781 INTEGER :: m_max_pad, bloc_size
783 mpi_comm,
DIMENSION(2) :: communicator
789 DO i = 1,
SIZE(list_mode)
792 DO m = 1, mesh%dom_me
795 c_in_loc(:,k) = c_in(j_loc,k,i)
803 hh=mesh%hloc_gauss(index)
804 hm=min(mesh%hm(i),hh)
809 ray = sum(mesh%rr(1,j_loc)*
ww(:,l))
812 c_in_gauss(index,1,i) = sum(c_in_loc(:,1)*
ww(:,l))
813 c_in_gauss(index,2,i) = sum(c_in_loc(:,2)*
ww(:,l))
816 gradc_in(index,1,i) = sum(c_in_loc(:,1)*dw_loc(1,:))*hh
817 gradc_in(index,2,i) = sum(c_in_loc(:,2)*dw_loc(1,:))*hh
818 gradc_in(index,3,i) = mode/ray*sum(c_in_loc(:,2)*
ww(:,l))*hm
819 gradc_in(index,4,i) = -mode/ray*sum(c_in_loc(:,1)*
ww(:,l))*hm
820 gradc_in(index,5,i) = sum(c_in_loc(:,1)*dw_loc(2,:))*hh
821 gradc_in(index,6,i) = sum(c_in_loc(:,2)*dw_loc(2,:))*hh
826 bloc_size =
SIZE(c_in_gauss,1)/nb_procs+1
828 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
833 mesh%hloc_gauss, coeff1_in_level, v_out, c_out, nb_procs, bloc_size, m_max_pad)
842 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
843 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: c1_in
844 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: c2_in
845 REAL(KIND=8),
INTENT(OUT) :: int_out
846 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
848 REAL(KIND=8) :: int_c1_in, int_c1_in_loc, int_c1_in_F
849 REAL(KIND=8) :: int_c2_in, int_c2_in_loc, int_c2_in_F
850 INTEGER :: m, l , i, index, code
852 mpi_comm,
DIMENSION(2) :: communicator
857 DO i = 1,
SIZE(list_mode)
859 IF (list_mode(i)==0)
THEN 860 DO m = 1, mesh%dom_me
862 DO l = 1, mesh%gauss%l_G
865 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
868 int_c1_in_loc = int_c1_in_loc + c1_in(index,1,i)*ray*mesh%gauss%rj(l,m)
869 int_c2_in_loc = int_c2_in_loc + c2_in(index,1,i)*ray*mesh%gauss%rj(l,m)
875 int_c1_in_loc = int_c1_in_loc*2*acos(-1.d0)
876 CALL mpi_allreduce(int_c1_in_loc, int_c1_in_f, 1, mpi_double_precision, mpi_sum, &
877 communicator(2), code)
878 CALL mpi_allreduce(int_c1_in_f, int_c1_in, 1, mpi_double_precision, mpi_sum, &
879 communicator(1), code)
881 int_c2_in_loc = int_c2_in_loc*2*acos(-1.d0)
882 CALL mpi_allreduce(int_c2_in_loc, int_c2_in_f, 1, mpi_double_precision, mpi_sum, &
883 communicator(2), code)
884 CALL mpi_allreduce(int_c2_in_f, int_c2_in, 1, mpi_double_precision, mpi_sum, &
885 communicator(1), code)
887 IF (int_c2_in.LE.1.d-14)
THEN 890 int_out=-int_c1_in/int_c2_in
896 fcompr, ff_phi_1mphi, stab_mass)
903 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff, ff_gauss
904 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: level_set_ext
905 INTEGER ,
INTENT(IN) :: mode
906 INTEGER ,
INTENT(IN) ::
type 907 REAL(KIND=8),
DIMENSION(:,:),
INTENT(IN) :: fcompr
908 REAL(KIND=8),
DIMENSION(:),
INTENT(IN) :: ff_phi_1mphi
909 REAL(KIND=8),
INTENT(IN) :: stab_mass
910 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w) :: ff_loc, level_set_loc
911 INTEGER,
DIMENSION(mesh%gauss%n_w) :: jj_loc
912 REAL(KIND=8),
DIMENSION(3) :: fcomprl
913 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w) :: v_loc
914 INTEGER,
DIMENSION(mesh%gauss%n_w) :: idxm
915 INTEGER :: i, m, l, ni, iglob, index
916 REAL(KIND=8) :: fl, ray
919 petscerrorcode :: ierr
921 CALL vecset(vect, 0.d0, ierr)
924 DO m = 1, mesh%dom_me
925 jj_loc = mesh%jj(:,m)
927 level_set_loc = level_set_ext(jj_loc)
929 DO ni = 1, mesh%gauss%n_w
931 iglob = la%loc_to_glob(1,i)
936 DO l = 1, mesh%gauss%l_G
939 DO ni = 1, mesh%gauss%n_w
941 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
945 fl = (sum(ff_loc(:)*mesh%gauss%ww(:,l)) + ff_gauss(index)+stab_mass*ff_phi_1mphi(index))*mesh%gauss%rj(l,m)*ray
950 fcomprl(1) = fcompr(index,1)*mesh%gauss%rj(l,m)*ray
951 fcomprl(2) = -mode*fcompr(index,4)*mesh%gauss%rj(l,m)
952 fcomprl(3) = fcompr(index,5)*mesh%gauss%rj(l,m)*ray
954 ELSE IF (type==2)
THEN 955 fcomprl(1) = fcompr(index,2)*mesh%gauss%rj(l,m)*ray
956 fcomprl(2) = mode*fcompr(index,3)*mesh%gauss%rj(l,m)
957 fcomprl(3) = fcompr(index,6)*mesh%gauss%rj(l,m)*ray
959 CALL error_petsc(
'error in type while calling qs_00_level_set_gauss')
962 DO ni = 1, mesh%gauss%n_w
964 v_loc(ni) = v_loc(ni) + mesh%gauss%ww(ni,l) *fl
966 v_loc(ni) = v_loc(ni) + (fcomprl(1)*mesh%gauss%dw(1,ni,l,m) + fcomprl(2)*mesh%gauss%ww(ni,l) &
967 + fcomprl(3)*mesh%gauss%dw(2,ni,l,m))
971 CALL vecsetvalues(vect, mesh%gauss%n_w, idxm, v_loc, add_values, ierr)
973 CALL vecassemblybegin(vect,ierr)
974 CALL vecassemblyend(vect,ierr)
subroutine solver(my_ksp, b, x, reinit, verbose)
subroutine, public fft_par_compr_entro_visc_dcl(communicator, V1_in, V2_in, c1_in, c_in_real, hloc_gauss, coeff1_in_level, V_out, c_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine, public extract(xghost, ks, ke, LA, phi)
subroutine twod_volume(communicator, mesh, RESLT)
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
subroutine error_petsc(string)
real(kind=8) function user_time()
subroutine smb_ugradc_gauss_fft_par(communicator, mesh, list_mode, V_in, c_in, nb_procs, c_out)
subroutine scalar_with_bc_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_js_D, pp_mode_global_js_D)
subroutine, public three_level_level_set(comm_one_d, time, cc_1_LA, dt, list_mode, cc_mesh, cn_m1, cn, chmp_vit, max_vel, my_par_cc, cc_list_dirichlet_sides, cc_per, nb_inter, visc_entro_level)
subroutine smb_compression_gauss_fft_par(communicator, mesh, list_mode, V_in, c_in, nb_procs, c_out)
subroutine compute_int_mass_correct(communicator, mesh, list_mode, c1_in, c2_in, int_out)
subroutine smb_compr_visc_entro_gauss_fft_par(communicator, mesh, list_mode, c_in, c_reg_in, visc_entro_real, coeff1_in_level, nb_procs, V_out, c_out)
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine dirichlet_rhs(js_D, bs_D, b)
subroutine smb_visc_entro_gauss_fft_par(communicator, mesh, list_mode, c_in, visc_entro_real, coeff1_in_level, nb_procs, V_out, c_out)
real(kind=8) max_velocity_at_tn
subroutine dirichlet_nodes_parallel(mesh, list_dirichlet_sides, js_d)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine qs_00(mesh, LA, ff, vect)
subroutine, public fft_compression_level_set_dcl(communicator_F, communicator_S, V1_in, V2_in, c_in, c_out, hloc_gauss, l_G, nb_procs, bloc_size, m_max_pad, temps)
subroutine qs_00_level_set_gauss(mesh, LA, ff, ff_gauss, mode, type, vect, level_set_ext, fcompr, ff_phi_1mphi, stab_mass)
subroutine dirichlet_m_parallel(matrix, glob_js_D)
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
subroutine, public fft_par_entro_visc_dcl(communicator, V1_in, c1_in, c_in_real, hloc_gauss, coeff1_in_level, V_out, c_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine qs_diff_mass_scal_m_level(mesh, LA, visco, mass, stab, i_mode, mode, matrix)
subroutine, public fft_par_dot_prod_dcl(communicator, V1_in, V2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
real(kind=8), dimension(:,:,:,:), pointer dw
real(kind=8), dimension(:,:), pointer ww
subroutine qs_regul_m(mesh, LA, stab, i_mode, mode, matrix)
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)