7 #include "petsc/finclude/petsc.h" 14 SUBROUTINE bdf2_ns_stress_bc_with_u(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, dt, Re, list_mode, pp_mesh, &
15 vv_mesh, incpn_m1, incpn, pn_m1, pn, un_m1, un, &
16 chmp_mag, bn_p2, opt_tempn)
39 REAL(KIND=8) :: time, dt, Re
40 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
41 TYPE(
mesh_type),
INTENT(IN) :: pp_mesh, vv_mesh
44 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: incpn_m1, incpn
45 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: pn_m1, pn
46 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: un_m1, un
47 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN),
OPTIONAL:: opt_tempn
48 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: chmp_mag, Bn_p2
49 REAL(KIND=8),
DIMENSION(vv_mesh%np,6,SIZE(list_mode)):: chmp_mag_aux
51 INTEGER,
SAVE :: m_max_c
54 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: pp_global_D
55 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: pp_mode_global_js_D
58 TYPE(
dyn_int_line),
DIMENSION(:),
POINTER,
SAVE :: vv_mode_global_js_D
59 TYPE(
dyn_real_line),
DIMENSION(:),
ALLOCATABLE,
SAVE :: vel_global_D
60 LOGICAL,
SAVE :: once = .true.
61 INTEGER,
SAVE :: my_petscworld_rank
62 REAL(KIND=8),
DIMENSION(:,:,:,:),
ALLOCATABLE,
SAVE :: visco_entro_sym_grad_u
66 INTEGER,
POINTER,
DIMENSION(:) :: pp_1_ifrom, vv_3_ifrom
68 INTEGER :: i, k, m, n, n1, n2, n3, n123
69 INTEGER :: nb_procs, code, nu_mat, mode
70 REAL(KIND=8) :: moyenne
72 REAL(KIND=8),
DIMENSION(pp_mesh%np, 2, SIZE(list_mode)) :: div
73 REAL(KIND=8),
DIMENSION(pp_mesh%np, 2) :: pn_p1, phi
74 REAL(KIND=8),
DIMENSION(vv_mesh%np, 6) :: un_p1
75 REAL(KIND=8),
DIMENSION(vv_mesh%np, 6, SIZE(list_mode)) :: un_m2
76 REAL(KIND=8),
DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rotv_v, rotb_b, rotb_b_aux
77 REAL(KIND=8),
DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: mag_force
78 REAL(KIND=8),
DIMENSION(vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: rhs_gauss
79 REAL(KIND=8),
DIMENSION(vv_mesh%np,6,SIZE(list_mode)) :: uext
80 REAL(KIND=8),
DIMENSION(vv_mesh%np,2,SIZE(list_mode)) :: nu_tilde
81 REAL(KIND=8),
DIMENSION(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,SIZE(list_mode)) :: two_nu_tilde_grads_u_ext
82 REAL(KIND=8),
DIMENSION(vv_mesh%np) :: vel_loc, vel_tot
83 REAL(KIND=8) :: tps, tps_tot, tps_cumul, coeff, vloc, cfl, cfl_max, norm
85 REAL(KIND=8) :: one, zero, three
86 DATA zero, one, three/0.d0,1.d0,3.d0/
89 petscerrorcode :: ierr
90 mpi_comm,
DIMENSION(:),
POINTER :: comm_one_d
91 mat,
DIMENSION(:),
POINTER,
SAVE :: vel_mat
92 mat,
DIMENSION(:),
POINTER,
SAVE :: press_mat
93 mat,
SAVE :: mass_mat, mass_mat0
94 vec,
SAVE :: vb_3_145, vb_3_236, vx_3, vx_3_ghost
95 vec,
SAVE :: pb_1, pb_2, px_1, px_1_ghost
96 ksp,
DIMENSION(:),
POINTER,
SAVE :: vel_ksp, press_ksp
97 ksp,
SAVE :: mass_ksp, mass_ksp0
102 CALL mpi_comm_rank(petsc_comm_world,my_petscworld_rank,code)
107 CALL veccreateghost(comm_one_d(1), n, &
108 petsc_determine,
SIZE(vv_3_ifrom), vv_3_ifrom, vx_3, ierr)
109 CALL vecghostgetlocalform(vx_3, vx_3_ghost, ierr)
110 CALL vecduplicate(vx_3, vb_3_145, ierr)
111 CALL vecduplicate(vx_3, vb_3_236, ierr)
115 CALL veccreateghost(comm_one_d(1), n, &
116 petsc_determine,
SIZE(pp_1_ifrom), pp_1_ifrom, px_1, ierr)
117 CALL vecghostgetlocalform(px_1, px_1_ghost, ierr)
118 CALL vecduplicate(px_1, pb_1, ierr)
119 CALL vecduplicate(px_1, pb_2, ierr)
123 m_max_c =
SIZE(list_mode)
131 ALLOCATE(pp_global_d(m_max_c))
133 ALLOCATE(pp_global_d(i)%DRL(
SIZE(pp_mode_global_js_d(i)%DIL)))
139 vv_js_d, vv_mode_global_js_d)
141 ALLOCATE(vel_global_d(m_max_c))
143 ALLOCATE(vel_global_d(i)%DRL(
SIZE(vv_mode_global_js_d(i)%DIL)))
150 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 154 IF (list_mode(i)==0) cycle
160 IF (minval(list_mode)==0)
THEN 163 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 167 IF (list_mode(i).NE.0) cycle
176 ALLOCATE(visco_entro_sym_grad_u(3,vv_mesh%gauss%l_G*vv_mesh%dom_me,6,
SIZE(list_mode)))
177 visco_entro_sym_grad_u = 0.d0
178 ALLOCATE(vel_mat(2*m_max_c),vel_ksp(2*m_max_c))
179 ALLOCATE(press_mat(m_max_c),press_ksp(m_max_c))
186 inputs%LES_coeff1,
inputs%stab_bdy_ns, i, mode, vel_mat(nu_mat))
187 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 189 vel_mat(nu_mat), vv_3_la)
192 CALL init_solver(
inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
197 inputs%LES_coeff1,
inputs%stab_bdy_ns, i, mode, vel_mat(nu_mat))
198 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 200 vel_mat(nu_mat), vv_3_la)
203 CALL init_solver(
inputs%my_par_vv,vel_ksp(nu_mat),vel_mat(nu_mat),comm_one_d(1),&
212 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 230 IF (
inputs%type_pb==
'mhd')
THEN 232 IF (
inputs%if_quasi_static_approx)
THEN 235 chmp_mag_aux(:,:,i) = h_b_quasi_static(
'H', vv_mesh%rr, mode)
240 chmp_mag_aux(:,:,i) = h_b_quasi_static(
'B', vv_mesh%rr, mode)
243 rotb_b = rotb_b + rotb_b_aux
247 rotv_v = rotv_v - rotb_b
250 IF (
inputs%type_pb==
'fhd')
THEN 252 IF (
inputs%if_helmholtz_force)
THEN 257 rotv_v = rotv_v - mag_force
260 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
264 IF (
inputs%verbose_CFL)
THEN 267 IF (list_mode(i)==0)
THEN 272 vel_loc = vel_loc + coeff*(un(:,1,i)**2+un(:,2,i)**2+un(:,3,i)**2 &
273 +un(:,4,i)**2+un(:,5,i)**2+un(:,6,i)**2)
275 CALL mpi_comm_size(comm_one_d(2),nb_procs,code)
276 CALL mpi_allreduce(vel_loc,vel_tot,vv_mesh%np,mpi_double_precision, mpi_sum, comm_one_d(2), code)
277 vel_tot = sqrt(vel_tot)
279 DO m = 1, vv_mesh%dom_me
280 vloc = maxval(vel_tot(vv_mesh%jj(:,m)))
282 cfl = max(vloc*dt/min(vv_mesh%hloc(m),maxval(vv_mesh%hm)),cfl)
284 CALL mpi_allreduce(cfl,cfl_max,1,mpi_double_precision, mpi_max, comm_one_d(1), code)
291 IF (
PRESENT(opt_tempn))
THEN 293 (4*un-un_m1)/(2*
inputs%dt), pn, (4.d0*incpn-incpn_m1)/3.d0, &
294 rotv_v, rhs_gauss, opt_tempn=opt_tempn)
297 (4*un-un_m1)/(2*
inputs%dt), pn, (4.d0*incpn-incpn_m1)/3.d0, &
303 IF (
inputs%if_variable_visco)
THEN 306 IF (.NOT. (
inputs%verbose_CFL))
THEN 307 CALL mpi_comm_size(comm_one_d(2), nb_procs, code)
310 nu_tilde, 2 * un - un_m1, two_nu_tilde_grads_u_ext)
319 pn_p1(:,:) = pn(:,:,i)
323 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236, &
324 opt_tensor=-visco_entro_sym_grad_u(:,:,:,i))
325 ELSE IF (
inputs%if_variable_visco)
THEN 326 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236, &
327 opt_tensor=-two_nu_tilde_grads_u_ext(:,:,:,i))
329 CALL rhs_3x3(vv_mesh, vv_3_la, mode, rhs_gauss(:,:,i), vb_3_145, vb_3_236)
332 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 333 CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_236, vv_3_la)
334 CALL periodic_rhs_petsc(vvz_per%n_bord, vvz_per%list, vvz_per%perlist, vb_3_145, vv_3_la)
337 n1 =
SIZE(vv_js_d(1)%DIL)
338 n2 =
SIZE(vv_js_d(2)%DIL)
339 n3 =
SIZE(vv_js_d(3)%DIL)
341 vel_global_d(i)%DRL(1:n1) = vv_exact(1,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
342 vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(4,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
343 vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(5,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
344 vel_global_d(i)%DRL(n123+1:) = 0.d0
345 CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_145)
346 vel_global_d(i)%DRL(1:n1) = vv_exact(2,vv_mesh%rr(:,vv_js_d(1)%DIL), mode, time)
347 vel_global_d(i)%DRL(n1+1:n1+n2) = vv_exact(3,vv_mesh%rr(:,vv_js_d(2)%DIL), mode, time)
348 vel_global_d(i)%DRL(n1+n2+1:n123)= vv_exact(6,vv_mesh%rr(:,vv_js_d(3)%DIL), mode, time)
349 CALL dirichlet_rhs(vv_mode_global_js_d(i)%DIL-1,vel_global_d(i)%DRL,vb_3_236)
350 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
358 CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
359 CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
360 CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,1))
361 CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,4))
362 CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,5))
367 CALL vecghostupdatebegin(vx_3,insert_values,scatter_forward,ierr)
368 CALL vecghostupdateend(vx_3,insert_values,scatter_forward,ierr)
369 CALL extract(vx_3_ghost,1,1,vv_3_la,un_p1(:,2))
370 CALL extract(vx_3_ghost,2,2,vv_3_la,un_p1(:,3))
371 CALL extract(vx_3_ghost,3,3,vv_3_la,un_p1(:,6))
373 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
382 IF (
inputs%my_periodic%nb_periodic_pairs/=0)
THEN 389 pp_global_d(i)%DRL = 0.d0
390 CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_1)
391 CALL dirichlet_rhs(pp_mode_global_js_d(i)%DIL-1,pp_global_d(i)%DRL,pb_2)
396 CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
397 CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
398 CALL extract(px_1_ghost,1,1,pp_1_la,phi(:,1))
401 CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
402 CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
403 CALL extract(px_1_ghost,1,1,pp_1_la,phi(:,2))
405 phi = -phi*(1.5d0/dt)
406 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
416 CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
417 CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
418 CALL extract(px_1_ghost,1,1,pp_1_la,div(:,1,i))
424 CALL vecghostupdatebegin(px_1,insert_values,scatter_forward,ierr)
425 CALL vecghostupdateend(px_1,insert_values,scatter_forward,ierr)
426 CALL extract(px_1_ghost,1,1,pp_1_la,div(:,2,i))
434 pn_p1(:,k) = pn_p1(:,k) + phi(:,k) - 2*div(:,k,i)/re -
inputs%div_stab_in_ns/re*div(:,k,i)
436 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
443 CALL moy(comm_one_d(1),pp_mesh, pn_p1(:,1),moyenne)
444 pn_p1(:,1) = pn_p1(:,1)-moyenne
458 un_m2(:,:,i) = un_m1(:,:,i)
461 un_m1(:,:,i) = un(:,:,i)
463 pn_m1(:,:,i) = pn(:,:,i)
465 incpn_m1(:,:,i) = incpn(:,:,i)
467 tps =
user_time() - tps; tps_cumul=tps_cumul+tps
472 IF (
inputs%verbose_divergence)
THEN 473 norm =
norm_sf(comm_one_d,
'H1', vv_mesh, list_mode, un)
480 IF (
PRESENT(opt_tempn))
THEN 482 un, un_m1, un_m2, pn_m1, rotv_v, visco_entro_sym_grad_u, opt_tempn)
485 un, un_m1, un_m2, pn_m1, rotv_v, visco_entro_sym_grad_u)
505 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
506 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: V_in
507 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(OUT) :: V_out
508 LOGICAL,
INTENT(IN) :: precession_in
509 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: RotV, W, Om
510 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
511 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
512 INTEGER :: m, l , i, mode, index, k
513 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: Vs, Omega_s
515 REAL(KIND=8) :: tps, PI
516 REAL(KIND=8),
DIMENSION(3) :: temps
517 REAL(KIND=8),
DIMENSION(:,:,:),
ALLOCATABLE,
SAVE :: Omega
518 LOGICAL,
SAVE :: once=.true.
520 INTEGER :: nb_procs, bloc_size, m_max_pad, code
521 mpi_comm :: communicator
525 ALLOCATE(omega(mesh%np,6,
SIZE(list_mode)))
528 DO i=1,
SIZE(list_mode)
529 IF (list_mode(i) == 1)
THEN 531 omega(:,1,i) =
inputs%taux_precession * sin(
inputs%angle_s_pi*pi)
532 omega(:,4,i) = -
inputs%taux_precession * sin(
inputs%angle_s_pi*pi)
533 ELSE IF (list_mode(i) == 0)
THEN 535 omega(:,5,i) =
inputs%taux_precession * cos(
inputs%angle_s_pi*pi)
541 DO i = 1,
SIZE(list_mode)
547 vs(:,k) = v_in(j_loc,k,i)
548 omega_s(:,k) = omega(mesh%jj(:,m),k,i)
551 DO l = 1, mesh%gauss%l_G
553 dw_loc = mesh%gauss%dw(:,:,l,m)
556 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
559 w(index,1,i) = sum(vs(:,1)*mesh%gauss%ww(:,l))
560 w(index,3,i) = sum(vs(:,3)*mesh%gauss%ww(:,l))
561 w(index,5,i) = sum(vs(:,5)*mesh%gauss%ww(:,l))
563 w(index,2,i) = sum(vs(:,2)*mesh%gauss%ww(:,l))
564 w(index,4,i) = sum(vs(:,4)*mesh%gauss%ww(:,l))
565 w(index,6,i) = sum(vs(:,6)*mesh%gauss%ww(:,l))
567 om(index,1,i) = sum(omega_s(:,1)*mesh%gauss%ww(:,l))
568 om(index,3,i) = sum(omega_s(:,3)*mesh%gauss%ww(:,l))
569 om(index,5,i) = sum(omega_s(:,5)*mesh%gauss%ww(:,l))
571 om(index,2,i) = sum(omega_s(:,2)*mesh%gauss%ww(:,l))
572 om(index,4,i) = sum(omega_s(:,4)*mesh%gauss%ww(:,l))
573 om(index,6,i) = sum(omega_s(:,6)*mesh%gauss%ww(:,l))
576 rotv(index,1,i) = mode/ray*w(index,6,i) &
577 -sum(vs(:,3)*dw_loc(2,:))
578 rotv(index,4,i) = sum(vs(:,2)*dw_loc(2,:)) &
579 -sum(vs(:,6)*dw_loc(1,:))
580 rotv(index,5,i) = 1/ray*w(index,3,i) &
581 +sum(vs(:,3)*dw_loc(1,:)) &
582 -mode/ray*w(index,2,i)
584 rotv(index,2,i) =-mode/ray*w(index,5,i) &
585 -sum(vs(:,4)*dw_loc(2,:))
586 rotv(index,3,i) = sum(vs(:,1)*dw_loc(2,:)) &
587 -sum(vs(:,5)*dw_loc(1,:))
588 rotv(index,6,i) = 1/ray*w(index,4,i) &
589 +sum(vs(:,4)*dw_loc(1,:))&
590 +mode/ray*w(index,1,i)
596 IF (.NOT.precession_in)
THEN 602 rotv = rotv + 2.d0 * om
606 CALL mpi_comm_size(communicator, nb_procs, code)
607 bloc_size =
SIZE(rotv,1)/nb_procs+1
608 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
622 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
623 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: V_in, W_in
624 REAL(KIND=8),
DIMENSION(:,:,:) :: V_out
625 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: RotV, V, W
626 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
627 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
628 INTEGER :: m, l , i, mode, index, k
629 REAL(KIND=8),
DIMENSION(mesh%gauss%n_w,6) :: Vs, Ws
632 REAL(KIND=8),
DIMENSION(3) :: temps
634 INTEGER :: nb_procs, bloc_size, m_max_pad, code
635 mpi_comm :: communicator
638 DO i = 1,
SIZE(list_mode)
644 vs(:,k) = v_in(j_loc,k,i)
645 ws(:,k) = w_in(j_loc,k,i)
648 DO l = 1, mesh%gauss%l_G
650 dw_loc = mesh%gauss%dw(:,:,l,m)
653 ray = sum(mesh%rr(1,j_loc)*mesh%gauss%ww(:,l))
656 w(index,1,i) = sum(ws(:,1)*mesh%gauss%ww(:,l))
657 w(index,3,i) = sum(ws(:,3)*mesh%gauss%ww(:,l))
658 w(index,5,i) = sum(ws(:,5)*mesh%gauss%ww(:,l))
660 w(index,2,i) = sum(ws(:,2)*mesh%gauss%ww(:,l))
661 w(index,4,i) = sum(ws(:,4)*mesh%gauss%ww(:,l))
662 w(index,6,i) = sum(ws(:,6)*mesh%gauss%ww(:,l))
663 v(index,1,i) = sum(vs(:,1)*mesh%gauss%ww(:,l))
664 v(index,3,i) = sum(vs(:,3)*mesh%gauss%ww(:,l))
665 v(index,5,i) = sum(vs(:,5)*mesh%gauss%ww(:,l))
667 v(index,2,i) = sum(vs(:,2)*mesh%gauss%ww(:,l))
668 v(index,4,i) = sum(vs(:,4)*mesh%gauss%ww(:,l))
669 v(index,6,i) = sum(vs(:,6)*mesh%gauss%ww(:,l))
672 rotv(index,1,i) = mode/ray*v(index,6,i) &
673 -sum(vs(:,3)*dw_loc(2,:))
674 rotv(index,4,i) = sum(vs(:,2)*dw_loc(2,:)) &
675 -sum(vs(:,6)*dw_loc(1,:))
676 rotv(index,5,i) = 1/ray*v(index,3,i) &
677 +sum(vs(:,3)*dw_loc(1,:)) &
678 -mode/ray*v(index,2,i)
680 rotv(index,2,i) =-mode/ray*v(index,5,i) &
681 -sum(vs(:,4)*dw_loc(2,:))
682 rotv(index,3,i) = sum(vs(:,1)*dw_loc(2,:)) &
683 -sum(vs(:,5)*dw_loc(1,:))
684 rotv(index,6,i) = 1/ray*v(index,4,i) &
685 +sum(vs(:,4)*dw_loc(1,:))&
686 +mode/ray*v(index,1,i)
697 CALL mpi_comm_size(communicator, nb_procs, code)
698 bloc_size =
SIZE(rotv,1)/nb_procs+1
699 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
718 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
719 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: scal_in, vect_in
720 REAL(KIND=8),
DIMENSION(:,:,:) :: vect_out
721 REAL(KIND=8),
DIMENSION(mesh%np,2,SIZE(list_mode)) :: chi_coeff, vect_in_square
722 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: chi_coeff_gauss
723 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: grad_vect_in_square
724 INTEGER :: i, mode, index, m, k, l
725 INTEGER,
DIMENSION(:,:),
POINTER :: jj
726 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
729 INTEGER :: nb_procs, bloc_size, m_max_pad, code
730 mpi_comm :: communicator
735 CALL mpi_comm_size(communicator, nb_procs, code)
736 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
740 bloc_size =
SIZE(chi_coeff)/nb_procs+1
741 CALL fft_par_scal_funct(communicator, chi_coeff, chi_coeff_law, nb_procs, bloc_size, m_max_pad)
746 DO i = 1,
SIZE(list_mode)
748 DO m = 1, mesh%dom_me
752 chi_coeff_gauss(index,k,i) = sum(chi_coeff(jj(:,m),k,i)*
ww(:,l))
759 bloc_size =
SIZE(vect_in,1)/nb_procs+1
760 CALL fft_par_dot_prod_dcl(communicator, vect_in, vect_in, vect_in_square, nb_procs, bloc_size, m_max_pad)
763 DO i = 1,
SIZE(list_mode)
766 DO m = 1, mesh%dom_me
770 rad = sum(mesh%rr(1,jj(:,m))*
ww(:,l))
771 grad_vect_in_square(index,1,i) = sum(vect_in_square(jj(:,m),1,i)*dw_loc(1,:))
772 grad_vect_in_square(index,2,i) = sum(vect_in_square(jj(:,m),2,i)*dw_loc(1,:))
773 grad_vect_in_square(index,3,i) = mode/rad * sum(vect_in_square(jj(:,m),2,i)*
ww(:,l))
774 grad_vect_in_square(index,4,i) = - mode/rad * sum(vect_in_square(jj(:,m),1,i)*
ww(:,l))
775 grad_vect_in_square(index,5,i) = sum(vect_in_square(jj(:,m),1,i)*dw_loc(2,:))
776 grad_vect_in_square(index,6,i) = sum(vect_in_square(jj(:,m),2,i)*dw_loc(2,:))
782 bloc_size =
SIZE(chi_coeff_gauss,1)/nb_procs+1
783 CALL fft_scalar_vect_dcl(communicator, 0.5*grad_vect_in_square, chi_coeff_gauss, vect_out, 1, nb_procs, bloc_size, m_max_pad)
794 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
795 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(IN) :: scal_in, vect_in
796 REAL(KIND=8),
DIMENSION(:,:,:) :: vect_out
797 REAL(KIND=8),
DIMENSION(mesh%np,2,SIZE(list_mode)) :: vect_in_square, chi_coeff
798 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,2,SIZE(list_mode)) :: vect_in_square_gauss
799 REAL(KIND=8),
DIMENSION(mesh%gauss%l_G*mesh%me,6,SIZE(list_mode)) :: grad_chi_coeff
800 INTEGER :: i, mode, index, m, k, l
801 INTEGER,
DIMENSION(:,:),
POINTER :: jj
802 REAL(KIND=8),
DIMENSION(mesh%gauss%k_d,mesh%gauss%n_w) :: dw_loc
805 INTEGER :: nb_procs, bloc_size, m_max_pad, code
806 mpi_comm :: communicator
811 CALL mpi_comm_size(communicator, nb_procs, code)
812 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
815 bloc_size =
SIZE(vect_in,1)/nb_procs+1
816 CALL fft_par_dot_prod_dcl(communicator, vect_in, vect_in, vect_in_square, nb_procs, bloc_size, m_max_pad)
821 DO i = 1,
SIZE(list_mode)
823 DO m = 1, mesh%dom_me
827 vect_in_square_gauss(index,k,i) = sum(vect_in_square(jj(:,m),k,i)*
ww(:,l))
835 bloc_size =
SIZE(chi_coeff)/nb_procs+1
836 CALL fft_par_scal_funct(communicator, chi_coeff, chi_coeff_law, nb_procs, bloc_size, m_max_pad)
839 DO i = 1,
SIZE(list_mode)
842 DO m = 1, mesh%dom_me
846 rad = sum(mesh%rr(1,jj(:,m))*
ww(:,l))
847 grad_chi_coeff(index,1,i) = sum(chi_coeff(jj(:,m),1,i)*dw_loc(1,:))
848 grad_chi_coeff(index,2,i) = sum(chi_coeff(jj(:,m),2,i)*dw_loc(1,:))
849 grad_chi_coeff(index,3,i) = mode/rad * sum(chi_coeff(jj(:,m),2,i)*
ww(:,l))
850 grad_chi_coeff(index,4,i) = - mode/rad * sum(chi_coeff(jj(:,m),1,i)*
ww(:,l))
851 grad_chi_coeff(index,5,i) = sum(chi_coeff(jj(:,m),1,i)*dw_loc(2,:))
852 grad_chi_coeff(index,6,i) = sum(chi_coeff(jj(:,m),2,i)*dw_loc(2,:))
858 bloc_size =
SIZE(grad_chi_coeff,1)/nb_procs+1
859 CALL fft_scalar_vect_dcl(communicator, grad_chi_coeff, -0.5*vect_in_square_gauss, vect_out, 1, nb_procs, bloc_size, m_max_pad)
863 SUBROUTINE smb_nu_tilde_fft_par(communicator,list_mode,scal_inout) ! MODIFICATION: computation of nu tilde function of temperature
869 INTEGER,
DIMENSION(:),
INTENT(IN) :: list_mode
870 REAL(KIND=8),
DIMENSION(:,:,:),
INTENT(INOUT) :: scal_inout
872 INTEGER :: nb_procs, bloc_size, m_max_pad, code
873 mpi_comm :: communicator
876 CALL mpi_comm_size(communicator, nb_procs, code)
877 m_max_pad = 3*
SIZE(list_mode)*nb_procs/2
879 bloc_size =
SIZE(scal_inout)/nb_procs+1
880 CALL fft_par_scal_funct(communicator, scal_inout, nu_tilde_law, nb_procs, bloc_size, m_max_pad)
884 SUBROUTINE moy(communicator,mesh,p,RESLT)
890 REAL(KIND=8),
DIMENSION(:) ,
INTENT(IN) :: p
891 REAL(KIND=8) ,
INTENT(OUT) :: RESLT
892 REAL(KIND=8) :: vol_loc, vol_out, r_loc, r_out
893 INTEGER :: m, l , i , ni, code
894 INTEGER,
DIMENSION(mesh%gauss%n_w) :: j_loc
896 mpi_comm :: communicator
900 DO m = 1, mesh%dom_me
902 DO l = 1, mesh%gauss%l_G
904 DO ni = 1, mesh%gauss%n_w; i = j_loc(ni)
905 ray = ray + mesh%rr(1,i)*mesh%gauss%ww(ni,l)
908 r_loc = r_loc + sum(p(j_loc(:))*mesh%gauss%ww(:,l))*ray*mesh%gauss%rj(l,m)
909 vol_loc = vol_loc + ray*mesh%gauss%rj(l,m)
914 CALL mpi_allreduce(r_loc,r_out,1,mpi_double_precision, mpi_sum, communicator, code)
915 CALL mpi_allreduce(vol_loc,vol_out,1,mpi_double_precision, mpi_sum, communicator, code)
916 reslt = r_out / vol_out
subroutine solver(my_ksp, b, x, reinit, verbose)
subroutine smb_cross_prod_gauss_sft_par(communicator, mesh, list_mode, V_in, V_out, precession_in)
subroutine smb_helmholtz_force_gauss_fft_par(communicator, mesh, list_mode, scal_in, vect_in, vect_out)
subroutine, public extract(xghost, ks, ke, LA, phi)
subroutine, public create_my_ghost(mesh, LA, ifrom)
subroutine smb_kelvin_force_gauss_fft_par(communicator, mesh, list_mode, scal_in, vect_in, vect_out)
subroutine, public fft_par_scal_funct(communicator, c1_inout, funct, nb_procs, bloc_size, m_max_pad, temps)
subroutine create_local_petsc_matrix(communicator, LA, matrix, clean)
real(kind=8) function user_time()
subroutine, public fft_par_cross_prod_dcl(communicator, V1_in, V2_in, V_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine, public smb_explicit_diffu_sym(communicator, mesh, list_mode, nb_procs, visc_dyn, vel, V_out)
subroutine qs_diff_mass_scal_m(mesh, LA, 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, public rhs_ns_gauss_3x3(vv_mesh, pp_mesh, communicator, list_mode, time, V1m, pn, pn_inc, rotv_v, rhs_gauss, opt_tempn)
subroutine qs_01_div_hybrid_generic(type_fe_velocity, vv_mesh, pp_mesh, pp_LA, mode, gg, pb_1, pb_2)
subroutine init_solver(my_par, my_ksp, matrix, communicator, solver, precond, opt_re_init)
subroutine dirichlet_rhs(js_D, bs_D, b)
subroutine vector_glob_js_d(vv_mesh, list_mode, vv_3_LA, vv_list_dirichlet_sides, vv_js_D, vv_mode_global_js_D)
subroutine, public periodic_matrix_petsc(n_bord, list, perlist, matrix, LA)
subroutine, public bdf2_ns_stress_bc_with_u(comm_one_d, time, vv_3_LA, pp_1_LA, vvz_per, pp_per, dt, Re, list_mode, pp_mesh, vv_mesh, incpn_m1, incpn, pn_m1, pn, un_m1, un, chmp_mag, Bn_p2, opt_tempn)
subroutine dirichlet_m_parallel(matrix, glob_js_D)
type(my_verbose), public talk_to_me
real(kind=8) function norm_sf(communicator, norm_type, mesh, list_mode, v)
subroutine, public fft_scalar_vect_dcl(communicator, V1_in, V2_in, V_out, pb, nb_procs, bloc_size, m_max_pad, temps)
subroutine rhs_3x3(mesh, vv_3_LA, mode, rhs_in, vb_145, vb_236, opt_tensor, opt_tensor_scaln_bdy)
subroutine, public compute_entropy_viscosity(comm_one_d, vv_3_LA, vv_mesh, pp_mesh, time, list_mode, vvz_per, un, un_m1, un_m2, pn_m1, rotv_v_m1, visco_entro_grad_u, opt_tempn)
subroutine moy(communicator, mesh, p, RESLT)
subroutine, public fft_par_dot_prod_dcl(communicator, V1_in, V2_in, c_out, nb_procs, bloc_size, m_max_pad, temps)
subroutine smb_nu_tilde_fft_par(communicator, list_mode, scal_inout)
real(kind=8), dimension(:,:,:,:), pointer dw
real(kind=8), dimension(:,:), pointer ww
subroutine, public periodic_rhs_petsc(n_bord, list, perlist, v_rhs, LA)
subroutine smb_curlh_cross_b_gauss_sft_par(communicator, mesh, list_mode, V_in, W_in, V_out)
subroutine scalar_without_glob_js_d(pp_mesh, list_mode, pp_1_LA, pp_mode_global_js_D)