From 68b2ad19ae99f3dd94110b99a2046b36f0c6e35f Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Tue, 5 Nov 2024 11:36:48 -0500 Subject: [PATCH 1/7] add viscous and surface_tension flags --- src/common/m_mpi_common.fpp | 4 +- src/common/m_variables_conversion.fpp | 8 +- src/post_process/m_checker.fpp | 2 +- src/post_process/m_global_parameters.fpp | 6 +- src/post_process/m_start_up.f90 | 2 + src/pre_process/m_assign_variables.fpp | 2 +- src/pre_process/m_global_parameters.fpp | 6 +- src/pre_process/m_mpi_proxy.fpp | 2 +- src/pre_process/m_start_up.fpp | 2 + src/simulation/m_data_output.fpp | 30 +++--- src/simulation/m_global_parameters.fpp | 26 ++++-- src/simulation/m_mpi_proxy.fpp | 7 +- src/simulation/m_rhs.fpp | 54 +++++------ src/simulation/m_riemann_solvers.fpp | 112 +++++++++++------------ src/simulation/m_sim_helpers.f90 | 12 +-- src/simulation/m_start_up.fpp | 12 ++- src/simulation/m_time_steppers.fpp | 2 +- src/simulation/m_viscous.fpp | 20 ++-- 18 files changed, 166 insertions(+), 143 deletions(-) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index b02180f00..42748e57b 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -245,7 +245,7 @@ contains MPI_DOUBLE_PRECISION, MPI_MAX, 0, & MPI_COMM_WORLD, ierr) - if (any(Re_size > 0)) then + if (viscous) then call MPI_REDUCE(vcfl_max_loc, vcfl_max_glb, 1, & MPI_DOUBLE_PRECISION, MPI_MAX, 0, & MPI_COMM_WORLD, ierr) @@ -258,7 +258,7 @@ contains icfl_max_glb = icfl_max_loc - if (any(Re_size > 0)) then + if (viscous) then vcfl_max_glb = vcfl_max_loc Rc_min_glb = Rc_min_loc end if diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 7dfc4c07b..23538b9dd 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -340,7 +340,7 @@ contains #ifdef MFC_SIMULATION ! Computing the shear and bulk Reynolds numbers from species analogs - if (any(Re_size > 0)) then + if (viscous) then if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 do i = 1, 2 @@ -530,7 +530,7 @@ contains G_K = max(0d0, G_K) end if - if (any(Re_size > 0)) then + if (viscous) then do i = 1, 2 Re_K(i) = dflt_real @@ -596,7 +596,7 @@ contains qv_K = qvs(1) end if - if (any(Re_size > 0)) then + if (viscous) then if (num_fluids == 1) then ! need to consider case with num_fluids >= 2 do i = 1, 2 @@ -661,7 +661,7 @@ contains #ifdef MFC_SIMULATION - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) do i = 1, 2 do j = 1, Re_size(i) diff --git a/src/post_process/m_checker.fpp b/src/post_process/m_checker.fpp index 3ebdf3f87..2b2bf88e6 100644 --- a/src/post_process/m_checker.fpp +++ b/src/post_process/m_checker.fpp @@ -115,7 +115,7 @@ contains !> Checks constraints on surface tension parameters (cf_wrt and sigma) subroutine s_check_inputs_surface_tension - @:PROHIBIT(cf_wrt .and. f_is_default(sigma), & + @:PROHIBIT(cf_wrt .and. .not. surface_tension, & "cf_wrt can only be enabled if the surface coefficient is set") end subroutine s_check_inputs_surface_tension diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index e74f179ce..5814376e0 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -275,6 +275,7 @@ module m_global_parameters !> @name surface tension coefficient !> @{ real(kind(0d0)) :: sigma + logical :: surface_tension !> #} !> @name Index variables used for m_variables_conversion @@ -397,6 +398,7 @@ contains poly_sigma = dflt_real sigR = dflt_real sigma = dflt_real + surface_tension = .false. adv_n = .false. end subroutine s_assign_default_values_to_user_inputs @@ -534,7 +536,7 @@ contains sys_size = stress_idx%end end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -559,7 +561,7 @@ contains sys_size = internalEnergies_idx%end alf_idx = 1 ! dummy, cannot actually have a void fraction - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 946d6cfc5..17e69bc72 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -114,6 +114,8 @@ subroutine s_read_input_file if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true. + if (.not. f_is_default(sigma)) surface_tension = .true. + else call s_mpi_abort('File post_process.inp is missing. Exiting ...') end if diff --git a/src/pre_process/m_assign_variables.fpp b/src/pre_process/m_assign_variables.fpp index f9280d86d..a2b2598bb 100644 --- a/src/pre_process/m_assign_variables.fpp +++ b/src/pre_process/m_assign_variables.fpp @@ -663,7 +663,7 @@ contains end do end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then q_prim_vf(c_idx)%sf(j, k, l) = eta*patch_icpp(patch_id)%cf_val + & (1d0 - eta)*patch_icpp(smooth_patch_id)%cf_val end if diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 4ce08a0c4..b399c6779 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -228,6 +228,7 @@ module m_global_parameters !> @name Surface Tension Modeling !> @{ real(kind(0d0)) :: sigma + logical :: surface_tension !> @} !> @name Index variables used for m_variables_conversion @@ -400,6 +401,7 @@ contains Re_inv = dflt_real Web = dflt_real poly_sigma = dflt_real + surface_tension = .false. adv_n = .false. @@ -615,7 +617,7 @@ contains sys_size = stress_idx%end end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -639,7 +641,7 @@ contains internalEnergies_idx%end = adv_idx%end + num_fluids sys_size = internalEnergies_idx%end - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 09c20034e..26b2d188e 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -56,7 +56,7 @@ contains & 'perturb_flow', 'perturb_sph', 'mixlayer_vel_profile', & & 'mixlayer_perturb', 'bubbles', 'polytropic', 'polydisperse', & & 'qbmm', 'file_per_process', 'adv_n', 'ib' , 'cfl_adap_dt', & - & 'cfl_const_dt', 'cfl_dt'] + & 'cfl_const_dt', 'cfl_dt', 'surface_tension'] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor call MPI_BCAST(fluid_rho(1), num_fluids_max, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 7712db0a0..362e10901 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -170,6 +170,8 @@ contains if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true. + if (.not. f_is_default(sigma)) surface_tension = .true. + else call s_mpi_abort('File pre_process.inp is missing. Exiting ...') end if diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index bb66c3722..69171c12d 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -148,7 +148,7 @@ contains ! Generating table header for the stability criteria to be outputted if (cfl_dt) then - if (any(Re_size > 0)) then + if (viscous) then write (1, '(A)') '==== Time-steps ====== dt ===== Time ======= ICFL '// & 'Max ==== VCFL Max ====== Rc Min =======' else @@ -156,7 +156,7 @@ contains '============== ICFL Max =============' end if else - if (any(Re_size > 0)) then + if (viscous) then write (1, '(A)') '==== Time-steps ====== Time ======= ICFL '// & 'Max ==== VCFL Max ====== Rc Min =======' else @@ -266,7 +266,7 @@ contains ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0d0, c) - if (any(Re_size > 0)) then + if (viscous) then call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf, vcfl_sf, Rc_sf) else call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl_sf) @@ -282,13 +282,13 @@ contains #ifdef CRAY_ACC_WAR !$acc update host(icfl_sf) - if (any(Re_size > 0)) then + if (viscous) then !$acc update host(vcfl_sf, Rc_sf) end if icfl_max_loc = maxval(icfl_sf) - if (any(Re_size > 0)) then + if (viscous) then vcfl_max_loc = maxval(vcfl_sf) Rc_min_loc = minval(Rc_sf) end if @@ -297,7 +297,7 @@ contains icfl_max_loc = maxval(icfl_sf) !$acc end kernels - if (any(Re_size > 0)) then + if (viscous) then !$acc kernels vcfl_max_loc = maxval(vcfl_sf) Rc_min_loc = minval(Rc_sf) @@ -317,21 +317,21 @@ contains Rc_min_glb) else icfl_max_glb = icfl_max_loc - if (any(Re_size > 0)) vcfl_max_glb = vcfl_max_loc - if (any(Re_size > 0)) Rc_min_glb = Rc_min_loc + if (viscous) vcfl_max_glb = vcfl_max_loc + if (viscous) Rc_min_glb = Rc_min_loc end if ! Determining the stability criteria extrema over all the time-steps if (icfl_max_glb > icfl_max) icfl_max = icfl_max_glb - if (any(Re_size > 0)) then + if (viscous) then if (vcfl_max_glb > vcfl_max) vcfl_max = vcfl_max_glb if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb end if ! Outputting global stability criteria extrema at current time-step if (proc_rank == 0) then - if (any(Re_size > 0)) then + if (viscous) then write (1, '(6X,I8,F10.6,6X,6X,F10.6,6X,F9.6,6X,F9.6,6X,F10.6)') & t_step, dt, t_step*dt, icfl_max_glb, & vcfl_max_glb, & @@ -348,7 +348,7 @@ contains call s_mpi_abort('ICFL is greater than 1.0. Exiting ...') end if - if (any(Re_size > 0)) then + if (viscous) then if (vcfl_max_glb /= vcfl_max_glb) then call s_mpi_abort('VCFL is NaN. Exiting ...') elseif (vcfl_max_glb > 1d0) then @@ -1578,8 +1578,8 @@ contains write (3, '(A)') '' write (3, '(A,F9.6)') 'ICFL Max: ', icfl_max - if (any(Re_size > 0)) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max - if (any(Re_size > 0)) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min + if (viscous) write (3, '(A,F9.6)') 'VCFL Max: ', vcfl_max + if (viscous) write (3, '(A,F10.6)') 'Rc Min: ', Rc_min call cpu_time(run_time) @@ -1615,7 +1615,7 @@ contains @:ALLOCATE_GLOBAL(icfl_sf(0:m, 0:n, 0:p)) icfl_max = 0d0 - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(vcfl_sf(0:m, 0:n, 0:p)) @:ALLOCATE_GLOBAL(Rc_sf (0:m, 0:n, 0:p)) @@ -1652,7 +1652,7 @@ contains ! Deallocating the ICFL, VCFL, CCFL, and Rc stability criteria @:DEALLOCATE_GLOBAL(icfl_sf) - if (any(Re_size > 0)) then + if (viscous) then @:DEALLOCATE_GLOBAL(vcfl_sf, Rc_sf) end if diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index dd5808ebb..fd2499c84 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -163,6 +163,9 @@ module m_global_parameters logical :: hypoelasticity !< hypoelasticity modeling logical, parameter :: chemistry = .${chemistry}$. !< Chemistry modeling logical :: cu_tensor + logical :: viscous !< Viscous effects + logical :: shear_stress !< Shear stresses + logical :: bulk_stress !< Bulk stresses !$acc declare create(chemistry) @@ -459,7 +462,8 @@ module m_global_parameters !> @name Surface tension parameters !> @{ real(kind(0d0)) :: sigma - !$acc declare create(sigma) + logical :: surface_tension + !$acc declare create(sigma, surface_tension) !> @} integer :: momxb, momxe @@ -562,6 +566,9 @@ contains weno_flat = .true. riemann_flat = .true. rdma_mpi = .false. + viscous = .false. + shear_stress = .false. + bulk_stress = .false. #:if not MFC_CASE_OPTIMIZATION mapped_weno = .false. @@ -650,6 +657,7 @@ contains ! Surface tension sigma = dflt_real + surface_tension = .false. ! Cuda aware MPI cu_tensor = .false. @@ -888,7 +896,7 @@ contains sys_size = stress_idx%end end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -906,7 +914,7 @@ contains internalEnergies_idx%end = adv_idx%end + num_fluids sys_size = internalEnergies_idx%end - if (.not. f_is_default(sigma)) then + if (surface_tension) then c_idx = sys_size + 1 sys_size = c_idx end if @@ -973,11 +981,15 @@ contains if (fluid_pp(i)%Re(2) > 0) Re_size(2) = Re_size(2) + 1 end do - !$acc update device(Re_size) + if (any(Re_size > 0d0)) viscous = .true. + if (Re_size(1) > 0d0) shear_stress = .true. + if (Re_size(2) > 0d0) bulk_stress = .true. + + !$acc update device(Re_size, viscous, shear_stress, bulk_stress) ! Bookkeeping the indexes of any viscous fluids and any pairs of ! fluids whose interface will support effects of surface tension - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_idx(1:2, 1:maxval(Re_size))) @@ -1049,7 +1061,7 @@ contains ! sufficient boundary conditions data as to iterate the solution in ! the physical computational domain from one time-step iteration to ! the next one - if (any(Re_size > 0)) then + if (viscous) then buff_size = 2*weno_polyn + 2 ! else if (hypoelasticity) then !TODO: check if necessary ! buff_size = 2*weno_polyn + 2 @@ -1193,7 +1205,7 @@ contains ! Deallocating the variables bookkeeping the indexes of any viscous ! fluids and any pairs of fluids whose interfaces supported effects ! of surface tension - if (any(Re_size > 0)) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_idx) end if diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 7da0b04d6..bc0552b1d 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -143,7 +143,7 @@ contains v_size = sys_size end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then nVars = num_dims + 1 if (n > 0) then if (p > 0) then @@ -203,7 +203,8 @@ contains & 'polydisperse', 'qbmm', 'acoustic_source', 'probe_wrt', 'integral_wrt', & & 'prim_vars_wrt', 'weno_avg', 'file_per_process', 'relax', & & 'adv_n', 'adap_dt', 'ib', 'bodyForces', 'bf_x', 'bf_y', 'bf_z', & - & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt' ] + & 'cfl_adap_dt', 'cfl_const_dt', 'cfl_dt', 'surface_tension', & + & 'viscous', 'shear_stress', 'bulk_stress' ] call MPI_BCAST(${VAR}$, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) #:endfor @@ -2343,7 +2344,7 @@ contains @:DEALLOCATE_GLOBAL(ib_buff_send, ib_buff_recv) end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then @:DEALLOCATE_GLOBAL(c_divs_buff_send, c_divs_buff_recv) end if diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index ae541ec64..0945728eb 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -246,7 +246,7 @@ contains @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end)) end do - if (.not. f_is_default(sigma)) then + if (surface_tension) then ! This assumes that the color function advection equation is ! the last equation. If this changes then this logic will ! need updated @@ -274,14 +274,14 @@ contains !$acc enter data attach(q_prim_qp%vf(l)%sf) end do - if (.not. f_is_default(sigma)) then + if (surface_tension) then q_prim_qp%vf(c_idx)%sf => & q_cons_qp%vf(c_idx)%sf !$acc enter data copyin(q_prim_qp%vf(c_idx)%sf) !$acc enter data attach(q_prim_qp%vf(c_idx)%sf) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(tau_Re_vf(1:sys_size)) do i = 1, num_dims @:ALLOCATE(tau_Re_vf(cont_idx%end + i)%sf(idwbuff(1)%beg:idwbuff(1)%end, & @@ -377,11 +377,11 @@ contains @:ALLOCATE_GLOBAL(dq_prim_dy_qp(1:1)) @:ALLOCATE_GLOBAL(dq_prim_dz_qp(1:1)) - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE(dq_prim_dx_qp(1)%vf(1:sys_size)) @:ALLOCATE(dq_prim_dy_qp(1)%vf(1:sys_size)) @:ALLOCATE(dq_prim_dz_qp(1)%vf(1:sys_size)) - if (any(Re_size > 0)) then + if (viscous) then do l = mom_idx%beg, mom_idx%end @:ALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf( & @@ -446,7 +446,7 @@ contains @:ALLOCATE_GLOBAL(dqR_prim_dy_n(1:num_dims)) @:ALLOCATE_GLOBAL(dqR_prim_dz_n(1:num_dims)) - if (any(Re_size > 0)) then + if (viscous) then do i = 1, num_dims @:ALLOCATE(dqL_prim_dx_n(i)%vf(1:sys_size)) @:ALLOCATE(dqL_prim_dy_n(i)%vf(1:sys_size)) @@ -455,7 +455,7 @@ contains @:ALLOCATE(dqR_prim_dy_n(i)%vf(1:sys_size)) @:ALLOCATE(dqR_prim_dz_n(i)%vf(1:sys_size)) - if (any(Re_size > 0)) then + if (viscous) then do l = mom_idx%beg, mom_idx%end @:ALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf( & @@ -502,7 +502,7 @@ contains end if ! END: Allocation/Association of d K_prim_ds_n ================== - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then @:ALLOCATE_GLOBAL(dqL_rsx_vf(idwbuff(1)%beg:idwbuff(1)%end, & idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end)) @@ -568,7 +568,7 @@ contains & idwbuff(3)%beg:idwbuff(3)%end)) end do - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then do l = mom_idx%beg, E_idx @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & @@ -643,11 +643,11 @@ contains end do !$acc update device(gamma_min, pres_inf) - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) end if - if (any(Re_size > 0)) then + if (viscous) then do i = 1, 2 do j = 1, Re_size(i) Res(i, j) = fluid_pp(Re_idx(i, j))%Re(i) @@ -793,7 +793,7 @@ contains if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb, rhs_pb, mv, rhs_mv, idwbuff(1), idwbuff(2), idwbuff(3), nbub) call nvtxStartRange("Viscous") - if (any(Re_size > 0)) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + if (viscous) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n, & qL_prim, & qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & @@ -805,7 +805,7 @@ contains call nvtxEndRange call nvtxStartRange("Surface_Tension") - if (.not. f_is_default(sigma)) call s_get_capilary(q_prim_qp%vf) + if (surface_tension) call s_get_capilary(q_prim_qp%vf) call nvtxEndRange ! Dimensional Splitting Loop ======================================= @@ -815,7 +815,7 @@ contains call nvtxStartRange("RHS-WENO") - if (f_is_default(sigma)) then + if (.not. surface_tension) then ! Reconstruct densitiess iv%beg = 1; iv%end = sys_size call s_reconstruct_cell_boundary_values( & @@ -927,7 +927,7 @@ contains ! RHS additions for viscosity call nvtxStartRange("RHS_add_phys") - if (any(Re_size > 0d0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then call s_compute_additional_physics_rhs(id, & q_prim_qp%vf, & rhs_vf, & @@ -1604,7 +1604,7 @@ contains if (idir == 1) then ! x-direction - if (.not. f_is_default(sigma)) then + if (surface_tension) then !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p do k = 0, n @@ -1636,7 +1636,7 @@ contains elseif (idir == 2) then ! y-direction - if (.not. f_is_default(sigma)) then + if (surface_tension) then !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p do k = 0, n @@ -1652,7 +1652,7 @@ contains end if if (cyl_coord .and. ((bc_y%beg == -2) .or. (bc_y%beg == -14))) then - if (any(Re_size > 0)) then + if (viscous) then if (p > 0) then call s_compute_viscous_stress_tensor(q_prim_vf, & dq_prim_dx_vf(mom_idx%beg:mom_idx%end), & @@ -1736,7 +1736,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(2) gang vector default(present) do l = 0, p do j = 0, m @@ -1771,7 +1771,7 @@ contains elseif (idir == 3) then ! z-direction - if (.not. f_is_default(sigma)) then + if (surface_tension) then !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p do k = 0, n @@ -2035,7 +2035,7 @@ contains pi_inf = pi_inf + alpha(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re(i) = dflt_real @@ -2248,7 +2248,7 @@ contains @:DEALLOCATE_GLOBAL(qL_rsz_vf, qR_rsz_vf) end if - if (any(Re_size > 0) .and. weno_Re_flux) then + if (viscous .and. weno_Re_flux) then @:DEALLOCATE_GLOBAL(dqL_rsx_vf, dqR_rsx_vf) if (n > 0) then @@ -2265,7 +2265,7 @@ contains deallocate (alf_sum%sf) end if - if (any(Re_size > 0)) then + if (viscous) then do l = mom_idx%beg, mom_idx%end @:DEALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf) end do @@ -2289,9 +2289,9 @@ contains @:DEALLOCATE(dq_prim_dz_qp(1)%vf) end if - if (any(Re_size > 0)) then + if (viscous) then do i = num_dims, 1, -1 - if (any(Re_size > 0)) then + if (viscous) then do l = mom_idx%beg, mom_idx%end @:DEALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf) @@ -2339,7 +2339,7 @@ contains @:DEALLOCATE(flux_gsrc_n(i)%vf(l)%sf) end do - if (any(Re_size > 0)) then + if (viscous) then do l = mom_idx%beg, E_idx @:DEALLOCATE(flux_src_n(i)%vf(l)%sf) end do @@ -2363,7 +2363,7 @@ contains @:DEALLOCATE_GLOBAL(flux_n, flux_src_n, flux_gsrc_n) - if (any(Re_size > 0) .and. cyl_coord) then + if (viscous .and. cyl_coord) then do i = 1, num_dims @:DEALLOCATE(tau_re_vf(cont_idx%end + i)%sf) end do diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 6e5123746..7c0b03c9d 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -450,7 +450,7 @@ contains qv_R = qv_R + alpha_rho_R(i)*qvs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real @@ -584,7 +584,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, c_sum_Yi_Phi, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -821,7 +821,7 @@ contains #:endfor - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then call s_compute_viscous_source_flux( & @@ -1077,7 +1077,7 @@ contains alpha_R(i) = qR_prim_rs${XYZ}$_vf(j + 1, k, l, advxb + i - 1) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real @@ -1131,7 +1131,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, 0d0, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -1197,7 +1197,7 @@ contains end do flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_L + pres_L)*vel_L(dir_idx(1)) - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qL_prim_rs${XYZ}$_vf(j, k, l, c_idx)*s_S end if @@ -1231,7 +1231,7 @@ contains end do flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_R + pres_R)*vel_R(dir_idx(1)) - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qR_prim_rs${XYZ}$_vf(j + 1, k, l, c_idx)*s_S end if @@ -1272,7 +1272,7 @@ contains end do flux_rs${XYZ}$_vf(j, k, l, E_idx) = (E_Star + p_Star)*s_S - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qL_prim_rs${XYZ}$_vf(j, k, l, c_idx)*s_S end if @@ -1315,7 +1315,7 @@ contains ! Compute the star velocities for the non-conservative terms end do - if (.not. f_is_default(sigma)) then + if (surface_tension) then flux_rs${XYZ}$_vf(j, k, l, c_idx) = & qR_prim_rs${XYZ}$_vf(j + 1, k, l, c_idx)*s_S end if @@ -1683,7 +1683,7 @@ contains qv_R = qvs(1) end if - if (any(Re_size > 0)) then + if (viscous) then if (num_fluids == 1) then ! Need to consider case with num_fluids >= 2 !$acc loop seq do i = 1, 2 @@ -1846,7 +1846,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, 0d0, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -2149,7 +2149,7 @@ contains qv_R = qv_R + qR_prim_rs${XYZ}$_vf(j + 1, k, l, i)*qvs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_L(i) = dflt_real @@ -2255,7 +2255,7 @@ contains call s_compute_speed_of_sound(pres_R, rho_avg, gamma_avg, pi_inf_R, H_avg, alpha_R, & vel_avg_rms, c_sum_Yi_Phi, c_avg) - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_avg_rs${XYZ}$_vf(j, k, l, i) = 2d0/(1d0/Re_L(i) + 1d0/Re_R(i)) @@ -2461,7 +2461,7 @@ contains #:endfor ! Computing HLLC flux and source flux for Euler system of equations - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then call s_compute_viscous_source_flux( & qL_prim_vf(momxb:momxe), & @@ -2487,7 +2487,7 @@ contains end if end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then call s_compute_capilary_source_flux( & q_prim_vf, & vel_src_rsx_vf, & @@ -2520,11 +2520,11 @@ contains end do !$acc update device(Gs) - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) end if - if (any(Re_size > 0)) then + if (viscous) then do i = 1, 2 do j = 1, Re_size(i) Res(i, j) = fluid_pp(Re_idx(i, j))%Re(i) @@ -2570,7 +2570,7 @@ contains @:ALLOCATE_GLOBAL(mom_sp_rsx_vf(is1%beg:is1%end + 1, is2%beg:is2%end, is3%beg:is3%end, 1:4)) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_avg_rsx_vf(is1%beg:is1%end, & is2%beg:is2%end, & is3%beg:is3%end, 1:2)) @@ -2598,7 +2598,7 @@ contains @:ALLOCATE_GLOBAL(mom_sp_rsy_vf(is1%beg:is1%end + 1, is2%beg:is2%end, is3%beg:is3%end, 1:4)) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_avg_rsy_vf(is1%beg:is1%end, & is2%beg:is2%end, & is3%beg:is3%end, 1:2)) @@ -2626,7 +2626,7 @@ contains @:ALLOCATE_GLOBAL(mom_sp_rsz_vf(is1%beg:is1%end + 1, is2%beg:is2%end, is3%beg:is3%end, 1:4)) end if - if (any(Re_size > 0)) then + if (viscous) then @:ALLOCATE_GLOBAL(Re_avg_rsz_vf(is1%beg:is1%end, & is2%beg:is2%end, & is3%beg:is3%end, 1:2)) @@ -2725,7 +2725,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe do l = isz%beg, isz%end @@ -2780,7 +2780,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe @@ -2839,7 +2839,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe @@ -2889,7 +2889,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe @@ -2942,7 +2942,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe do k = isy%beg, isy%end @@ -2986,7 +2986,7 @@ contains end do end do - if (any(Re_size > 0)) then + if (viscous) then !$acc parallel loop collapse(3) gang vector default(present) do i = momxb, momxe do k = isy%beg, isy%end @@ -3061,7 +3061,7 @@ contains if (norm_dir == 1) then - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then !$acc parallel loop collapse(4) gang vector default(present) do i = momxb, E_idx @@ -3094,7 +3094,7 @@ contains ! Reshaping Inputted Data in y-direction =========================== elseif (norm_dir == 2) then - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then !$acc parallel loop collapse(4) gang vector default(present) do i = momxb, E_idx do l = is3%beg, is3%end @@ -3125,7 +3125,7 @@ contains ! Reshaping Inputted Data in z-direction =========================== else - if (any(Re_size > 0) .or. (.not. f_is_default(sigma))) then + if (viscous .or. (surface_tension)) then !$acc parallel loop collapse(4) gang vector default(present) do i = momxb, E_idx do j = is1%beg, is1%end @@ -3216,7 +3216,7 @@ contains ! Viscous Stresses in z-direction ================================== if (norm_dir == 1) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3242,7 +3242,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3270,7 +3270,7 @@ contains if (n == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3312,7 +3312,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3344,7 +3344,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3385,7 +3385,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( avg_vel, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3415,7 +3415,7 @@ contains ! Viscous Stresses in r-direction ================================== elseif (norm_dir == 2) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end @@ -3465,7 +3465,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3500,7 +3500,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3545,7 +3545,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3575,7 +3575,7 @@ contains ! Viscous Stresses in theta-direction ================================== else - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3640,7 +3640,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(avg_vel, dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3740,7 +3740,7 @@ contains ! Viscous Stresses in x-direction ================================== if (norm_dir == 1) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3766,7 +3766,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3794,7 +3794,7 @@ contains if (n == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3835,7 +3835,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3863,7 +3863,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3903,7 +3903,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3933,7 +3933,7 @@ contains ! Viscous Stresses in y-direction ================================== elseif (norm_dir == 2) then - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -3978,7 +3978,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4009,7 +4009,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4050,7 +4050,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4080,7 +4080,7 @@ contains ! Viscous Stresses in z-direction ================================== else - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4137,7 +4137,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private( dvel_avg_dx, dvel_avg_dy, dvel_avg_dz, tau_Re) do l = isz%beg, isz%end do k = isy%beg, isy%end @@ -4363,7 +4363,7 @@ contains ! to convert mixture or species variables to the mixture variables ! s_convert_to_mixture_variables => null() - if (Re_size(1) > 0) then + if (shear_stress) then @:DEALLOCATE_GLOBAL(Re_avg_rsx_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsx_vf) @@ -4376,7 +4376,7 @@ contains if (n == 0) return - if (Re_size(1) > 0) then + if (shear_stress) then @:DEALLOCATE_GLOBAL(Re_avg_rsy_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsy_vf) @@ -4389,7 +4389,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then + if (shear_stress) then @:DEALLOCATE_GLOBAL(Re_avg_rsz_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsz_vf) diff --git a/src/simulation/m_sim_helpers.f90 b/src/simulation/m_sim_helpers.f90 index 9b1f1ce21..a4af91281 100644 --- a/src/simulation/m_sim_helpers.f90 +++ b/src/simulation/m_sim_helpers.f90 @@ -118,7 +118,7 @@ subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl dz(l)/(abs(vel(3)) + c)) end if - if (any(Re_size > 0)) then + if (viscous) then if (grid_geometry == 3) then vcfl_sf(j, k, l) = maxval(dt/Re_l/rho) & @@ -145,7 +145,7 @@ subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl icfl_sf(j, k, l) = dt/min(dx(j)/(abs(vel(1)) + c), & dy(k)/(abs(vel(2)) + c)) - if (any(Re_size > 0)) then + if (viscous) then vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/min(dx(j), dy(k))**2d0 @@ -159,7 +159,7 @@ subroutine s_compute_stability_from_dt(vel, c, rho, Re_l, j, k, l, icfl_sf, vcfl !1D icfl_sf(j, k, l) = (dt/dx(j))*(abs(vel(1)) + c) - if (any(Re_size > 0)) then + if (viscous) then vcfl_sf(j, k, l) = maxval(dt/Re_l/rho)/dx(j)**2d0 @@ -213,7 +213,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) dz(l)/(abs(vel(3)) + c)) end if - if (any(Re_size > 0)) then + if (viscous) then if (grid_geometry == 3) then vcfl_dt = cfl_target*(min(dx(j), dy(k), fltr_dtheta)**2d0) & /minval(1/(rho*Re_l)) @@ -228,7 +228,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) icfl_dt = cfl_target*min(dx(j)/(abs(vel(1)) + c), & dy(k)/(abs(vel(2)) + c)) - if (any(Re_size > 0)) then + if (viscous) then vcfl_dt = cfl_target*(min(dx(j), dy(k))**2d0)/maxval((1/Re_l)/rho) end if @@ -236,7 +236,7 @@ subroutine s_compute_dt_from_cfl(vel, c, max_dt, rho, Re_l, j, k, l) !1D icfl_dt = cfl_target*(dx(j)/(abs(vel(1)) + c)) - if (any(Re_size > 0)) then + if (viscous) then vcfl_dt = cfl_target*(dx(j)**2d0)/minval(1/(rho*Re_l)) end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index f40429b53..b0ac1eb01 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -198,6 +198,8 @@ contains if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true. + if (.not. f_is_default(sigma)) surface_tension = .true. + else call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') end if @@ -1325,13 +1327,13 @@ contains call s_initialize_acoustic_src() end if - if (any(Re_size > 0)) then + if (viscous) then call s_initialize_viscous_module() end if call s_initialize_rhs_module() - if (.not. f_is_default(sigma)) call s_initialize_surface_tension_module() + if (surface_tension) call s_initialize_surface_tension_module() #if defined(MFC_OpenACC) && defined(MFC_MEMORY_DUMP) call acc_present_dump() @@ -1470,7 +1472,7 @@ contains !$acc update device(R_n, R_v, phi_vn, phi_nv, Pe_c, Tw, pv, M_n, M_v, k_n, k_v, pb0, mass_n0, mass_v0, Pe_T, Re_trans_T, Re_trans_c, Im_trans_T, Im_trans_c, omegaN , mul0, ss, gamma_v, mu_v, gamma_m, gamma_n, mu_n, gam) !$acc update device(acoustic_source, num_source) - !$acc update device(sigma) + !$acc update device(sigma, surface_tension) !$acc update device(dx, dy, dz, x_cb, x_cc, y_cb, y_cc, z_cb, z_cc) @@ -1507,11 +1509,11 @@ contains call s_finalize_mpi_proxy_module() call s_finalize_global_parameters_module() if (relax) call s_finalize_relaxation_solver_module() - if (any(Re_size > 0)) then + if (viscous) then call s_finalize_viscous_module() end if - if (.not. f_is_default(sigma)) call s_finalize_surface_tension_module() + if (surface_tension) call s_finalize_surface_tension_module() if (bodyForces) call s_finalize_body_forces_module() ! Terminating MPI execution environment diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 97cc9fdac..3c282f47b 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -189,7 +189,7 @@ contains end do end if - if (.not. f_is_default(sigma)) then + if (surface_tension) then @:ALLOCATE(q_prim_vf(c_idx)%sf(idwbuff(1)%beg:idwbuff(1)%end, & idwbuff(2)%beg:idwbuff(2)%end, & idwbuff(3)%beg:idwbuff(3)%end)) diff --git a/src/simulation/m_viscous.fpp b/src/simulation/m_viscous.fpp index 139fccc7c..3e2c3bf70 100644 --- a/src/simulation/m_viscous.fpp +++ b/src/simulation/m_viscous.fpp @@ -95,7 +95,7 @@ contains end do end do end do - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -161,7 +161,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -202,7 +202,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -268,7 +268,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -306,7 +306,7 @@ contains if (p == 0) return - if (Re_size(1) > 0) then ! Shear stresses + if (shear_stress) then ! Shear stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -372,7 +372,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -414,7 +414,7 @@ contains end do end if - if (Re_size(2) > 0) then ! Bulk stresses + if (bulk_stress) then ! Bulk stresses !$acc parallel loop collapse(3) gang vector default(present) private(alpha_visc, alpha_rho_visc, Re_visc, tau_Re ) do l = is3_viscous%beg, is3_viscous%end do k = -1, 1 @@ -480,7 +480,7 @@ contains pi_inf_visc = pi_inf_visc + alpha_visc(i)*pi_infs(i) end do - if (any(Re_size > 0)) then + if (viscous) then !$acc loop seq do i = 1, 2 Re_visc(i) = dflt_real @@ -1023,7 +1023,7 @@ contains is1_viscous, is2_viscous, is3_viscous) end if - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then if (norm_dir == 2) then !$acc parallel loop collapse(4) gang vector default(present) @@ -1124,7 +1124,7 @@ contains is1_viscous, is2_viscous, is3_viscous) end if - if (any(Re_size > 0)) then + if (viscous) then if (weno_Re_flux) then if (norm_dir == 2) then !$acc parallel loop collapse(4) gang vector default(present) From 9c10f6c8e6b9425fa7a9e0c2681ec94000484f67 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Tue, 5 Nov 2024 12:23:24 -0500 Subject: [PATCH 2/7] format --- src/simulation/m_rhs.fpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 0945728eb..bea7aa49d 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -794,14 +794,14 @@ contains call nvtxStartRange("Viscous") if (viscous) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n, & - qL_prim, & - qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n, & - qR_prim, & - q_prim_qp, & - dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, & - idwbuff(1), idwbuff(2), idwbuff(3)) + dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n, & + qL_prim, & + qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n, & + qR_prim, & + q_prim_qp, & + dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, & + idwbuff(1), idwbuff(2), idwbuff(3)) call nvtxEndRange call nvtxStartRange("Surface_Tension") From a477c426f3e95709a0b24fd731564723ab60e509 Mon Sep 17 00:00:00 2001 From: wilfonba Date: Tue, 5 Nov 2024 12:54:04 -0500 Subject: [PATCH 3/7] fix frontier --- src/simulation/m_global_parameters.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index fd2499c84..e31eed308 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -186,7 +186,7 @@ module m_global_parameters !$acc declare create(num_dims, weno_polyn, weno_order, weno_num_stencils, num_fluids, wenojs, mapped_weno, wenoz, teno, wenoz_q) #:endif - !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach) + !$acc declare create(mpp_lim, model_eqns, mixture_err, alt_soundspeed, avg_state, mp_weno, weno_eps, teno_CT, hypoelasticity, low_Mach, viscous, shear_stress, bulk_stress) logical :: relax !< activate phase change integer :: relax_model !< Relaxation model From c8ad8c0b0ed430e25641332035d9e6711f1c91e2 Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Tue, 5 Nov 2024 15:19:16 -0500 Subject: [PATCH 4/7] requested changes --- docs/documentation/case.md | 6 ++++++ src/post_process/m_start_up.f90 | 4 +--- src/pre_process/m_start_up.fpp | 5 ++--- src/simulation/m_global_parameters.fpp | 1 - src/simulation/m_riemann_solvers.fpp | 6 +++--- src/simulation/m_start_up.fpp | 5 ++--- toolchain/mfc/run/case_dicts.py | 6 +++++- toolchain/mfc/test/cases.py | 2 +- 8 files changed, 20 insertions(+), 15 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index aaf463322..ab9635763 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -371,6 +371,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r | `n_start` | Integer | Save file from which to start simulation | | `t_save` | Real | Time duration between data output | | `t_stop` | Real | Simulation stop time | +| `surface_tension` | Logical | Activate surface tension. Requires `surface_tension = T` | +| `viscous` | Logical | Activate viscosity Requires `viscous = T` | - \* Options that work only with `model_eqns = 2`. - † Options that work only with ``cyl_coord = 'F'``. @@ -445,6 +447,10 @@ If this option is false, velocity gradient is computed using finite difference s - `weno_avg` it activates the arithmetic average of the left and right, WENO-reconstructed, cell-boundary values. This option requires `weno_Re_flux` to be true because cell boundary values are only utilized when employing the scalar divergence method in the computation of velocity gradients. +- `surface_tension` activates surface tension when set to `'T'`. Requires `sigma` to be set. + +- `viscous` activates viscosity when set to `'T'`. Requires `Re(1)` and `Re(2)` to be set. + #### Constant Time-Stepping - `dt` specifies the constant time step size that is used in simulation. diff --git a/src/post_process/m_start_up.f90 b/src/post_process/m_start_up.f90 index 17e69bc72..47832450a 100644 --- a/src/post_process/m_start_up.f90 +++ b/src/post_process/m_start_up.f90 @@ -83,7 +83,7 @@ subroutine s_read_input_file polydisperse, poly_sigma, file_per_process, relax, & relax_model, cf_wrt, sigma, adv_n, ib, & cfl_adap_dt, cfl_const_dt, t_save, t_stop, n_start, & - cfl_target + cfl_target, surface_tension ! Inquiring the status of the post_process.inp file file_loc = 'post_process.inp' @@ -114,8 +114,6 @@ subroutine s_read_input_file if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true. - if (.not. f_is_default(sigma)) surface_tension = .true. - else call s_mpi_abort('File post_process.inp is missing. Exiting ...') end if diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 362e10901..641752644 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -141,7 +141,8 @@ contains sigR, sigV, dist_type, rhoRV, R0_type, & file_per_process, relax, relax_model, & palpha_eps, ptgalpha_eps, ib, num_ibs, patch_ib, & - sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, n_start_old + sigma, adv_n, cfl_adap_dt, cfl_const_dt, n_start, & + n_start_old, surface_tension ! Inquiring the status of the pre_process.inp file file_loc = 'pre_process.inp' @@ -170,8 +171,6 @@ contains if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true. - if (.not. f_is_default(sigma)) surface_tension = .true. - else call s_mpi_abort('File pre_process.inp is missing. Exiting ...') end if diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index e31eed308..63fc21c73 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -981,7 +981,6 @@ contains if (fluid_pp(i)%Re(2) > 0) Re_size(2) = Re_size(2) + 1 end do - if (any(Re_size > 0d0)) viscous = .true. if (Re_size(1) > 0d0) shear_stress = .true. if (Re_size(2) > 0d0) bulk_stress = .true. diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 7c0b03c9d..682859b70 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -4363,7 +4363,7 @@ contains ! to convert mixture or species variables to the mixture variables ! s_convert_to_mixture_variables => null() - if (shear_stress) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_avg_rsx_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsx_vf) @@ -4376,7 +4376,7 @@ contains if (n == 0) return - if (shear_stress) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_avg_rsy_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsy_vf) @@ -4389,7 +4389,7 @@ contains if (p == 0) return - if (shear_stress) then + if (viscous) then @:DEALLOCATE_GLOBAL(Re_avg_rsz_vf) end if @:DEALLOCATE_GLOBAL(vel_src_rsz_vf) diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index b0ac1eb01..633f0e8ce 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -164,7 +164,8 @@ contains pi_fac, adv_n, adap_dt, bf_x, bf_y, bf_z, & k_x, k_y, k_z, w_x, w_y, w_z, p_x, p_y, p_z, & g_x, g_y, g_z, n_start, t_save, t_stop, & - cfl_adap_dt, cfl_const_dt, cfl_target + cfl_adap_dt, cfl_const_dt, cfl_target, & + viscous, surface_tension ! Checking that an input file has been provided by the user. If it ! has, then the input file is read in, otherwise, simulation exits. @@ -198,8 +199,6 @@ contains if (cfl_adap_dt .or. cfl_const_dt) cfl_dt = .true. - if (.not. f_is_default(sigma)) surface_tension = .true. - else call s_mpi_abort(trim(file_path)//' is missing. Exiting ...') end if diff --git a/toolchain/mfc/run/case_dicts.py b/toolchain/mfc/run/case_dicts.py index 7738ca65d..49be5732b 100644 --- a/toolchain/mfc/run/case_dicts.py +++ b/toolchain/mfc/run/case_dicts.py @@ -87,7 +87,8 @@ def analytic(self): 'num_ibs': ParamType.INT, 'cfl_dt': ParamType.LOG, 'n_start': ParamType.INT, - 'n_start_old': ParamType.INT + 'n_start_old': ParamType.INT, + 'surface_tension': ParamType.LOG, }) for ib_id in range(1, 10+1): @@ -225,6 +226,8 @@ def analytic(self): 't_save': ParamType.REAL, 'cfl_target': ParamType.REAL, 'low_Mach': ParamType.INT, + 'surface_tension': ParamType.LOG, + 'viscous': ParamType.LOG, }) for var in [ 'diffusion', 'reactions' ]: @@ -336,6 +339,7 @@ def analytic(self): 't_save': ParamType.REAL, 't_stop': ParamType.REAL, 'n_start': ParamType.INT, + 'surface_tension': ParamType.LOG, }) for cmp_id in range(1,3+1): diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index cf7208a01..99a9371ac 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -648,7 +648,7 @@ def alter_viscosity(dimInfo): # Viscosity & bubbles checks if len(dimInfo[0]) > 0: stack.push("Viscosity -> Bubbles", - {"fluid_pp(1)%Re(1)": 50, "bubbles": 'T'}) + {"fluid_pp(1)%Re(1)": 50, "bubbles": 'T', "viscous": 'T'}) stack.push('', { 'nb' : 1, 'fluid_pp(1)%gamma' : 0.16, 'fluid_pp(1)%pi_inf': 3515.0, From 35ec4c327b9fb822db9ad8c2cd388936c7dddc9e Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 5 Nov 2024 16:03:09 -0500 Subject: [PATCH 5/7] Update case.md --- docs/documentation/case.md | 80 +++++++++++++++++++------------------- 1 file changed, 40 insertions(+), 40 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index ab9635763..86f299c12 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -371,8 +371,8 @@ Details of implementation of viscosity in MFC can be found in [Coralic (2015)](r | `n_start` | Integer | Save file from which to start simulation | | `t_save` | Real | Time duration between data output | | `t_stop` | Real | Simulation stop time | -| `surface_tension` | Logical | Activate surface tension. Requires `surface_tension = T` | -| `viscous` | Logical | Activate viscosity Requires `viscous = T` | +| `surface_tension` | Logical | Activate surface tension | +| `viscous` | Logical | Activate viscosity | - \* Options that work only with `model_eqns = 2`. - † Options that work only with ``cyl_coord = 'F'``. @@ -447,37 +447,37 @@ If this option is false, velocity gradient is computed using finite difference s - `weno_avg` it activates the arithmetic average of the left and right, WENO-reconstructed, cell-boundary values. This option requires `weno_Re_flux` to be true because cell boundary values are only utilized when employing the scalar divergence method in the computation of velocity gradients. -- `surface_tension` activates surface tension when set to `'T'`. Requires `sigma` to be set. +- `surface_tension` activates surface tension when set to ``'T'``. Requires `sigma` to be set. -- `viscous` activates viscosity when set to `'T'`. Requires `Re(1)` and `Re(2)` to be set. +- `viscous` activates viscosity when set to ``'T'``. Requires `Re(1)` and `Re(2)` to be set. #### Constant Time-Stepping -- `dt` specifies the constant time step size that is used in simulation. -The value of `dt` needs to be sufficiently small such that the Courant-Friedrichs-Lewy (CFL) condition is satisfied. +- `dt` specifies the constant time step size used in the simulation. +The value of `dt` needs to be sufficiently small to satisfy the Courant-Friedrichs-Lewy (CFL) condition. -- `t_step_start` and `t_step_end` define the time steps at which simulation starts and ends, respectively. +- `t_step_start` and `t_step_end` define the time steps at which the simulation starts and ends. `t_step_save` is the time step interval for data output during simulation. To newly start the simulation, set `t_step_start = 0`. -To restart simulation from $k$-th time step, set `t_step_start = k`, see [Restarting Cases](running.md#restarting-cases). +To restart the simulation from $k$-th time step, set `t_step_start = k`; see [Restarting Cases](running.md#restarting-cases). ##### Adaptive Time-Stepping - `cfl_adap_dt` enables adaptive time stepping with a constant CFL when true -- `cfl_const_dt` enables constant dt time-stepping where dt results in a specified CFL for the initial condition +- `cfl_const_dt` enables constant `dt` time-stepping where `dt` results in a specified CFL for the initial condition - `cfl_target` specifies the target CFL value - `n_start` specifies the save file to start at -- `t_save` specifies the time interval between data output during simulation +- `t_save` specifies the time interval between data output during the simulation - `t_stop` specifies at what time the simulation should stop To newly start the simulation, set `n_start = 0`. -To restart simulation from $k$-th time step, see [Restarting Cases](running.md#restarting-cases). +To restart the simulation from $k$-th time step, see [Restarting Cases](running.md#restarting-cases). ### 7. Formatted Output @@ -518,23 +518,23 @@ The table lists formatted database output parameters. The parameters define vari - `parallel_io` activates parallel input/output (I/O) of data files. It is highly recommended to activate this option in a parallel environment. With parallel I/O, MFC inputs and outputs a single file throughout pre-process, simulation, and post-process, regardless of the number of processors used. -Parallel I/O enables the use of different number of processors in each of the processes (i.e., simulation data generated using 1000 processors can be post-processed using a single processor). +Parallel I/O enables the use of different numbers of processors in each process (e.g., simulation data generated using 1000 processors can be post-processed using a single processor). - `file_per_process` deactivates shared file MPI-IO and activates file per process MPI-IO. The default behavior is to use a shared file. -File per process is useful when running on 10's of thousands of ranks. +File per process is useful when running on >10K ranks. If `file_per_process` is true, then pre_process, simulation, and post_process must be run with the same number of ranks. -- `cons_vars_wrt` and `prim_vars_wrt` activate output of conservative and primitive state variables into the database, respectively. +- `cons_vars_wrt` and `prim_vars_wrt` activate the output of conservative and primitive state variables into the database. -- `[variable's name]_wrt` activates output of the each specified variable into the database. +- `[variable's name]_wrt` activates the output of each specified variable into the database. - `schlieren_alpha(i)` specifies the intensity of the numerical Schlieren of $i$-th component. -- `fd_order` specifies the order of the finite difference scheme that is used to compute the vorticity from the velocity field and the numerical schlieren from the density field by an integer of 1, 2, and 4. -`fd_order = 1`, `2`, and `4` correspond to the first, second, and fourth-order finite difference schemes, respectively. +- `fd_order` specifies the order of the finite difference scheme used to compute the vorticity from the velocity field and the numerical schlieren from the density field using an integer of 1, 2, and 4. +`fd_order = 1`, `2`, and `4` correspond to the first, second, and fourth-order finite difference schemes. -- `probe_wrt` activates output of state variables at coordinates specified by `probe(i)%[x;y,z]`. +- `probe_wrt` activates the output of state variables at coordinates specified by `probe(i)%[x;y,z]`. ### 8. Acoustic Source {#acoustic-source} @@ -597,15 +597,15 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon - `%%foc_length` specifies the focal length of the transducer for transducer waves. It is the distance from the transducer to the focal point. -- `%%aperture` specifies the aperture of the transducer. It is the diameter of the projection of the transducer arc onto the y-axis (2D) or spherical cap onto the y-z plane (3D). To simulate a transducer enclosing half of the circle/sphere, set the aperture to double the focal length. For transducer array, it is the total aperture of the array. +- `%%aperture` specifies the aperture of the transducer. It is the diameter of the projection of the transducer arc onto the y-axis (2D) or spherical cap onto the y-z plane (3D). Set the aperture to double the focal length to simulate a transducer enclosing half of the circle/sphere. For the transducer array, it is the total aperture of the array. - `%%num_elements` specifies the number of transducer elements in a transducer array. -- `%%element_on` specifies the element number of the transducer array that is on. The element number starts from 1. If all elements are on, set `%%element_on` to 0. +- `%%element_on` specifies the element number of the transducer array that is on. The element number starts from 1, if all elements are on, set `%%element_on` to 0. -- `%%element_spacing_angle` specifies the spacing angle between adjacent transducer in radians. The total aperture (`%%aperture`) is set, so each transducer element is smaller if `%%element_spacing_angle` is larger. +- `%%element_spacing_angle` specifies the spacing angle between adjacent transducers in radians. The total aperture (`%%aperture`) is set, so each transducer element is smaller if `%%element_spacing_angle` is larger. -- `%%element_polygon_ratio` specifies the ratio of the polygon side length to the aperture diameter of each transducer element in a circular 3D transducer array. The polygon side length is calculated by using the total aperture (`%%aperture`) as the circumcicle diameter, and `%%num_elements` as the number of sides of the polygon. The ratio is used specify the aperture size of each transducer element in the array, as a ratio of the total aperture. +- `%%element_polygon_ratio` specifies the ratio of the polygon side length to the aperture diameter of each transducer element in a circular 3D transducer array. The polygon side length is calculated by using the total aperture (`%%aperture`) as the circumcircle diameter and `%%num_elements` as the number of sides of the polygon. The ratio is used to specify the aperture size of each transducer element in the array as a ratio of the total aperture. - `%%rotate_angle` specifies the rotation angle of the 3D circular transducer array along the x-axis (principal axis). It is optional and defaults to 0. @@ -633,16 +633,16 @@ Details of the transducer acoustic source model can be found in [Maeda and Colon | `mu_v` † | Real | Viscosity | | `k_v` † | Real | Thermal conductivity | | `qbmm` | Logical | Quadrature by method of moments| -| `dist_type` | Integer | Joint probability density function for bubble radius and velocity (only when qbmm is true)| -| `sigR` | Real | Standard deviation for probability density function of bubble radius (only when qbmm is true) | -| `sigV` | Real | Standard deviation for probability density function of bubble velocity (only when qbmm is true) | -| `rhoRV` | Real | Correlation coefficient for joint probability density function of bubble radius and velocity (only when qbmm is true) | +| `dist_type` | Integer | Joint probability density function for bubble radius and velocity (only for ``qbmm = 'T'``)| +| `sigR` | Real | Standard deviation for the probability density function of bubble radius (only for ``qbmm = 'T'``) | +| `sigV` | Real | Standard deviation for the probability density function of bubble velocity (only for ``qbmm = 'T'``) | +| `rhoRV` | Real | Correlation coefficient for the joint probability density function of bubble radius and velocity (only for ``qbmm = 'T'``) | -These options work only for gas-liquid two component flows. +These options work only for gas-liquid two-component flows. Component indexes are required to be 1 for liquid and 2 for gas. -- \* These parameters should be pretended with patch index $1$ that is filled with liquid: `fluid_pp(1)%`. -- † These parameters should be pretended with patch indexes that are respectively filled with liquid and gas: `fluid_pp(1)%` and `fluid_pp(2)%`. +- \* These parameters should be prepended with patch index $1$ that is filled with liquid: `fluid_pp(1)%`. +- † These parameters should be prepended with patch indexes filled with liquid and gas: `fluid_pp(1)%` and `fluid_pp(2)%`. This table lists the ensemble-averaged bubble model parameters. @@ -652,9 +652,9 @@ This table lists the ensemble-averaged bubble model parameters. `bubble_model = 1`, `2`, and `3` correspond to the Gilmore, Keller-Miksis, and Rayleigh-Plesset models. - `polytropic` activates polytropic gas compression in the bubble. -When `polytropic` is set `False`, the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on ([Preston et al., 2007](references.md#Preston07)). +When ``polytropic = 'F'``, the gas compression is modeled as non-polytropic due to heat and mass transfer across the bubble wall with constant heat and mass transfer coefficients based on ([Preston et al., 2007](references.md#Preston07)). -- `polydisperse` activates polydispersity in the bubble model by means of a probability density function (PDF) of the equiliibrium bubble radius. +- `polydisperse` activates polydispersity in the bubble model through a probability density function (PDF) of the equilibrium bubble radius. - `thermal` specifies a model for heat transfer across the bubble interface by an integer from 1 through 3. `thermal = 1`, `2`, and `3` correspond to no heat transfer (adiabatic gas compression), isothermal heat transfer, and heat transfer with a constant heat transfer coefficient based on [Preston et al., 2007](references.md#Preston07), respectively. @@ -679,17 +679,17 @@ Implementation of the parameters into the model follow [Ando (2010)](references. - `dist_type` specifies the initial joint PDF of initial bubble radius and bubble velocity required in qbmm. `dist_type = 1` and `2` correspond to binormal and lognormal-normal distributions respectively. -- `sigR` specifies the standard deviation of the PDF of bubble radius required in qbmm. +- `sigR` specifies the standard deviation of the PDF of bubble radius required in the QBMM feature. -- `sigV` specifies the standard deviation of the PDF of bubble velocity required in qbmm. +- `sigV` specifies the standard deviation of the PDF of bubble velocity required in the QBMM feature. -- `rhoRV` specifies the correlation coefficient of the joint PDF of bubble radius and bubble velocity required in qbmm. +- `rhoRV` specifies the correlation coefficient of the joint PDF of bubble radius and bubble velocity required in the QBMM feature. ### 10. Velocity Field Setup | Parameter | Type | Description | | ---: | :----: | :--- | -| `perturb_flow` | Logical | Perturb the initlial velocity field by random noise | +| `perturb_flow` | Logical | Perturb the initial velocity field by random noise | | `perturb_flow_fluid` | Integer | Fluid density whose flow is to be perturbed | | `perturb_flow_mag` | Real | Set the magnitude of flow perturbations | | `perturb_sph` | Logical | Perturb the initial partial density by random noise | @@ -704,7 +704,7 @@ The parameters are optionally used to define initial velocity profiles and pertu - `perturb_flow` activates the perturbation of initial velocity by random noise. -- `perturb_flow_fluid` specifies the fluid component whose flow is to be perturbed. +- `perturb_flow_fluid` specifies the fluid component whose flow will be perturbed. - `perturb_flow` activates the perturbation of initial velocity by random noise. @@ -712,7 +712,7 @@ The parameters are optionally used to define initial velocity profiles and pertu - `perturb_sph_fluid` specifies the fluid component whose partial density is to be perturbed. -- `mixlayer_vel_profile` activates setting of the mean streamwise velocity to hyperbolic tangent profile. This option works only for 2D and 3D cases. +- `mixlayer_vel_profile` activates setting the mean streamwise velocity to a hyperbolic tangent profile. This option works only for 2D and 3D cases. - `mixlayer_vel_coef` is a parameter for the hyperbolic tangent profile of a mixing layer when `mixlayer_vel_profile = 'T'`. The mean streamwise velocity profile is given as: @@ -734,7 +734,7 @@ $$ u = patch\_icpp(1)\%vel(1) * tanh(y\_cc * mixlayer\_vel\_profile) $$ - `relax_model` Specifies the phase change model to be used: [5] enables pT-equilibrium, while [6] activates pTg-equilibrium (if criteria are met). -- `palpha_eps` Specifies the tolerance used for the Newton Solvers used in the pT-equilibrium model. +- `palpha_eps` Specifies the tolerance for the Newton Solvers used in the pT-equilibrium model. - `ptgalpha_eps` Specifies the tolerance used for the Newton Solvers used in the pTg-equilibrium model. @@ -846,7 +846,7 @@ Each patch requires a different set of parameters, which are also listed in this | 10 | Annular Transducer Array | 2D-Axisym | #9 requirements | | 11 | Circular Transducer Array | 3D | #7 requirements, `%%element_polygon_ratio`, and `%%rotate_angle` | -Details of the required parameters for each acoustic support type are listed in [Acoustic Source](#acoustic-source). +The required parameters for each acoustic support type are listed in [Acoustic Source](#acoustic-source). The acoustic support number (`#`) corresponds to the acoustic support type `Acoustic(i)%%support`, where $i$ is the acoustic source index. For each `%%parameter`, prepend the parameter with `acoustic(i)%`. @@ -855,7 +855,7 @@ Additional requirements for all acoustic support types: - `num_source` must be set to the total number of acoustic sources. -- `%%support` must be set to the acoustic support number listed in the table. +- `%%support` must be set to the acoustic support number in the table. - `%%dipole` is only supported for planar sources. From c62f460c6475c1ce383d47579c07dc497cf6bb80 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 5 Nov 2024 17:35:54 -0500 Subject: [PATCH 6/7] Update m_rhs.fpp --- src/simulation/m_rhs.fpp | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index bea7aa49d..adb5ccb48 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -248,8 +248,8 @@ contains if (surface_tension) then ! This assumes that the color function advection equation is - ! the last equation. If this changes then this logic will - ! need updated + ! the last equation. If this changes, then this logic will + ! need updating do l = adv_idx%end + 1, sys_size - 1 @:ALLOCATE(q_prim_qp%vf(l)%sf(idwbuff(1)%beg:idwbuff(1)%end, idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end)) end do @@ -455,8 +455,6 @@ contains @:ALLOCATE(dqR_prim_dy_n(i)%vf(1:sys_size)) @:ALLOCATE(dqR_prim_dz_n(i)%vf(1:sys_size)) - if (viscous) then - do l = mom_idx%beg, mom_idx%end @:ALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & @@ -494,15 +492,11 @@ contains end do end if - end if - @:ACC_SETUP_VFs(dqL_prim_dx_n(i), dqL_prim_dy_n(i), dqL_prim_dz_n(i)) @:ACC_SETUP_VFs(dqR_prim_dx_n(i), dqR_prim_dy_n(i), dqR_prim_dz_n(i)) end do - end if - ! END: Allocation/Association of d K_prim_ds_n ================== + ! END: Allocation/Association of d K_prim_ds_n ================== - if (viscous) then if (weno_Re_flux) then @:ALLOCATE_GLOBAL(dqL_rsx_vf(idwbuff(1)%beg:idwbuff(1)%end, & idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end)) @@ -535,11 +529,8 @@ contains idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end)) end if - end if end if - ! ================================================================== - ! Allocation of gm_alphaK_n ===================================== @:ALLOCATE_GLOBAL(gm_alphaL_n(1:num_dims)) @:ALLOCATE_GLOBAL(gm_alphaR_n(1:num_dims)) @@ -568,7 +559,7 @@ contains & idwbuff(3)%beg:idwbuff(3)%end)) end do - if (viscous .or. (surface_tension)) then + if (viscous .or. surface_tension) then do l = mom_idx%beg, E_idx @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & @@ -927,7 +918,7 @@ contains ! RHS additions for viscosity call nvtxStartRange("RHS_add_phys") - if (viscous .or. (surface_tension)) then + if (viscous .or. surface_tension) then call s_compute_additional_physics_rhs(id, & q_prim_qp%vf, & rhs_vf, & From 63e005eb75ccfb41d1e13823b535e0ea374bde3c Mon Sep 17 00:00:00 2001 From: Ben Wilfong <48168887+wilfonba@users.noreply.github.com> Date: Tue, 5 Nov 2024 17:45:41 -0500 Subject: [PATCH 7/7] format --- src/simulation/m_rhs.fpp | 2942 +++++++++++++++++++------------------- 1 file changed, 1471 insertions(+), 1471 deletions(-) diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index adb5ccb48..cef03b63d 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -455,42 +455,42 @@ contains @:ALLOCATE(dqR_prim_dy_n(i)%vf(1:sys_size)) @:ALLOCATE(dqR_prim_dz_n(i)%vf(1:sys_size)) + do l = mom_idx%beg, mom_idx%end + @:ALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + @:ALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do + + if (n > 0) then do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf( & + @:ALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & & idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) - @:ALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf( & + @:ALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & & idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) end do + end if - if (n > 0) then - do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - @:ALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - end do - end if - - if (p > 0) then - do l = mom_idx%beg, mom_idx%end - @:ALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - @:ALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - end do - end if + if (p > 0) then + do l = mom_idx%beg, mom_idx%end + @:ALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + @:ALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do + end if @:ACC_SETUP_VFs(dqL_prim_dx_n(i), dqL_prim_dy_n(i), dqL_prim_dz_n(i)) @:ACC_SETUP_VFs(dqR_prim_dx_n(i), dqR_prim_dy_n(i), dqR_prim_dz_n(i)) @@ -529,1218 +529,1253 @@ contains idwbuff(2)%beg:idwbuff(2)%end, idwbuff(3)%beg:idwbuff(3)%end, mom_idx%beg:mom_idx%end)) end if - end if - - ! Allocation of gm_alphaK_n ===================================== - @:ALLOCATE_GLOBAL(gm_alphaL_n(1:num_dims)) - @:ALLOCATE_GLOBAL(gm_alphaR_n(1:num_dims)) - ! ================================================================== + end if - ! Allocation/Association of flux_n, flux_src_n, and flux_gsrc_n === - @:ALLOCATE_GLOBAL(flux_n(1:num_dims)) - @:ALLOCATE_GLOBAL(flux_src_n(1:num_dims)) - @:ALLOCATE_GLOBAL(flux_gsrc_n(1:num_dims)) + ! Allocation of gm_alphaK_n ===================================== + @:ALLOCATE_GLOBAL(gm_alphaL_n(1:num_dims)) + @:ALLOCATE_GLOBAL(gm_alphaR_n(1:num_dims)) + ! ================================================================== - do i = 1, num_dims + ! Allocation/Association of flux_n, flux_src_n, and flux_gsrc_n === + @:ALLOCATE_GLOBAL(flux_n(1:num_dims)) + @:ALLOCATE_GLOBAL(flux_src_n(1:num_dims)) + @:ALLOCATE_GLOBAL(flux_gsrc_n(1:num_dims)) - @:ALLOCATE(flux_n(i)%vf(1:sys_size)) - @:ALLOCATE(flux_src_n(i)%vf(1:sys_size)) - @:ALLOCATE(flux_gsrc_n(i)%vf(1:sys_size)) + do i = 1, num_dims - if (i == 1) then - do l = 1, sys_size - @:ALLOCATE(flux_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - @:ALLOCATE(flux_gsrc_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - end do + @:ALLOCATE(flux_n(i)%vf(1:sys_size)) + @:ALLOCATE(flux_src_n(i)%vf(1:sys_size)) + @:ALLOCATE(flux_gsrc_n(i)%vf(1:sys_size)) - if (viscous .or. surface_tension) then - do l = mom_idx%beg, E_idx - @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & + if (i == 1) then + do l = 1, sys_size + @:ALLOCATE(flux_n(i)%vf(l)%sf( & & idwbuff(1)%beg:idwbuff(1)%end, & & idwbuff(2)%beg:idwbuff(2)%end, & & idwbuff(3)%beg:idwbuff(3)%end)) + @:ALLOCATE(flux_gsrc_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) end do - end if - @:ALLOCATE(flux_src_n(i)%vf(adv_idx%beg)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) + if (viscous .or. surface_tension) then + do l = mom_idx%beg, E_idx + @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do + end if - if (riemann_solver == 1) then - do l = adv_idx%beg + 1, adv_idx%end - @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) - end do - end if + @:ALLOCATE(flux_src_n(i)%vf(adv_idx%beg)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) - if (chemistry) then - do l = chemxb, chemxe - @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & - & idwbuff(1)%beg:idwbuff(1)%end, & - & idwbuff(2)%beg:idwbuff(2)%end, & - & idwbuff(3)%beg:idwbuff(3)%end)) + if (riemann_solver == 1) then + do l = adv_idx%beg + 1, adv_idx%end + @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do + end if + + if (chemistry) then + do l = chemxb, chemxe + @:ALLOCATE(flux_src_n(i)%vf(l)%sf( & + & idwbuff(1)%beg:idwbuff(1)%end, & + & idwbuff(2)%beg:idwbuff(2)%end, & + & idwbuff(3)%beg:idwbuff(3)%end)) + end do + end if + + else + do l = 1, sys_size + @:ALLOCATE(flux_gsrc_n(i)%vf(l)%sf( & + idwbuff(1)%beg:idwbuff(1)%end, & + idwbuff(2)%beg:idwbuff(2)%end, & + idwbuff(3)%beg:idwbuff(3)%end)) end do end if - else - do l = 1, sys_size - @:ALLOCATE(flux_gsrc_n(i)%vf(l)%sf( & - idwbuff(1)%beg:idwbuff(1)%end, & - idwbuff(2)%beg:idwbuff(2)%end, & - idwbuff(3)%beg:idwbuff(3)%end)) - end do - end if + @:ACC_SETUP_VFs(flux_n(i), flux_src_n(i), flux_gsrc_n(i)) - @:ACC_SETUP_VFs(flux_n(i), flux_src_n(i), flux_gsrc_n(i)) + if (i == 1) then + if (riemann_solver /= 1) then + do l = adv_idx%beg + 1, adv_idx%end + flux_src_n(i)%vf(l)%sf => flux_src_n(i)%vf(adv_idx%beg)%sf - if (i == 1) then - if (riemann_solver /= 1) then - do l = adv_idx%beg + 1, adv_idx%end - flux_src_n(i)%vf(l)%sf => flux_src_n(i)%vf(adv_idx%beg)%sf + !$acc enter data attach(flux_src_n(i)%vf(l)%sf) + end do + end if + else + do l = 1, sys_size + flux_n(i)%vf(l)%sf => flux_n(1)%vf(l)%sf + flux_src_n(i)%vf(l)%sf => flux_src_n(1)%vf(l)%sf - !$acc enter data attach(flux_src_n(i)%vf(l)%sf) + !$acc enter data attach(flux_n(i)%vf(l)%sf,flux_src_n(i)%vf(l)%sf) end do end if - else - do l = 1, sys_size - flux_n(i)%vf(l)%sf => flux_n(1)%vf(l)%sf - flux_src_n(i)%vf(l)%sf => flux_src_n(1)%vf(l)%sf - - !$acc enter data attach(flux_n(i)%vf(l)%sf,flux_src_n(i)%vf(l)%sf) - end do - end if - end do + end do - ! END: Allocation/Association of flux_n, flux_src_n, and flux_gsrc_n === + ! END: Allocation/Association of flux_n, flux_src_n, and flux_gsrc_n === - if (alt_soundspeed) then - @:ALLOCATE_GLOBAL(blkmod1(0:m, 0:n, 0:p), blkmod2(0:m, 0:n, 0:p), alpha1(0:m, 0:n, 0:p), alpha2(0:m, 0:n, 0:p), Kterm(0:m, 0:n, 0:p)) - end if + if (alt_soundspeed) then + @:ALLOCATE_GLOBAL(blkmod1(0:m, 0:n, 0:p), blkmod2(0:m, 0:n, 0:p), alpha1(0:m, 0:n, 0:p), alpha2(0:m, 0:n, 0:p), Kterm(0:m, 0:n, 0:p)) + end if - @:ALLOCATE_GLOBAL(gamma_min(1:num_fluids), pres_inf(1:num_fluids)) + @:ALLOCATE_GLOBAL(gamma_min(1:num_fluids), pres_inf(1:num_fluids)) - do i = 1, num_fluids - gamma_min(i) = 1d0/fluid_pp(i)%gamma + 1d0 - pres_inf(i) = fluid_pp(i)%pi_inf/(1d0 + fluid_pp(i)%gamma) - end do - !$acc update device(gamma_min, pres_inf) + do i = 1, num_fluids + gamma_min(i) = 1d0/fluid_pp(i)%gamma + 1d0 + pres_inf(i) = fluid_pp(i)%pi_inf/(1d0 + fluid_pp(i)%gamma) + end do + !$acc update device(gamma_min, pres_inf) - if (viscous) then - @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) - end if + if (viscous) then + @:ALLOCATE_GLOBAL(Res(1:2, 1:maxval(Re_size))) + end if - if (viscous) then - do i = 1, 2 - do j = 1, Re_size(i) - Res(i, j) = fluid_pp(Re_idx(i, j))%Re(i) + if (viscous) then + do i = 1, 2 + do j = 1, Re_size(i) + Res(i, j) = fluid_pp(Re_idx(i, j))%Re(i) + end do end do - end do - !$acc update device(Res, Re_idx, Re_size) - end if + !$acc update device(Res, Re_idx, Re_size) + end if - ! Associating procedural pointer to the subroutine that will be - ! utilized to calculate the solution of a given Riemann problem - if (riemann_solver == 1) then - s_riemann_solver => s_hll_riemann_solver - elseif (riemann_solver == 2) then - s_riemann_solver => s_hllc_riemann_solver - end if + ! Associating procedural pointer to the subroutine that will be + ! utilized to calculate the solution of a given Riemann problem + if (riemann_solver == 1) then + s_riemann_solver => s_hll_riemann_solver + elseif (riemann_solver == 2) then + s_riemann_solver => s_hllc_riemann_solver + end if - ! Associating the procedural pointer to the appropriate subroutine - ! that will be utilized in the conversion to the mixture variables - if (model_eqns == 1) then ! Gamma/pi_inf model - s_convert_to_mixture_variables => & - s_convert_mixture_to_mixture_variables - else if (bubbles) then ! Volume fraction for bubbles - s_convert_to_mixture_variables => & - s_convert_species_to_mixture_variables_bubbles - else ! Volume fraction model - s_convert_to_mixture_variables => & - s_convert_species_to_mixture_variables - end if + ! Associating the procedural pointer to the appropriate subroutine + ! that will be utilized in the conversion to the mixture variables + if (model_eqns == 1) then ! Gamma/pi_inf model + s_convert_to_mixture_variables => & + s_convert_mixture_to_mixture_variables + else if (bubbles) then ! Volume fraction for bubbles + s_convert_to_mixture_variables => & + s_convert_species_to_mixture_variables_bubbles + else ! Volume fraction model + s_convert_to_mixture_variables => & + s_convert_species_to_mixture_variables + end if - !$acc parallel loop collapse(4) gang vector default(present) - do id = 1, num_dims - do i = 1, sys_size - do l = startz, p - startz - do k = starty, n - starty - do j = startx, m - startx - flux_gsrc_n(id)%vf(i)%sf(j, k, l) = 0d0 + !$acc parallel loop collapse(4) gang vector default(present) + do id = 1, num_dims + do i = 1, sys_size + do l = startz, p - startz + do k = starty, n - starty + do j = startx, m - startx + flux_gsrc_n(id)%vf(i)%sf(j, k, l) = 0d0 + end do end do end do end do end do - end do - if (bubbles) then - @:ALLOCATE_GLOBAL(nbub(0:m, 0:n, 0:p)) - end if + if (bubbles) then + @:ALLOCATE_GLOBAL(nbub(0:m, 0:n, 0:p)) + end if - end subroutine s_initialize_rhs_module + end subroutine s_initialize_rhs_module - subroutine s_compute_rhs(q_cons_vf, q_prim_vf, rhs_vf, pb, rhs_pb, mv, rhs_mv, t_step, time_avg) + subroutine s_compute_rhs(q_cons_vf, q_prim_vf, rhs_vf, pb, rhs_pb, mv, rhs_mv, t_step, time_avg) - type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf - type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf - type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb - real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv - integer, intent(in) :: t_step - real(kind(0d0)), intent(inout) :: time_avg + type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf + type(scalar_field), dimension(sys_size), intent(inout) :: q_prim_vf + type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf + real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: pb, rhs_pb + real(kind(0d0)), dimension(startx:, starty:, startz:, 1:, 1:), intent(inout) :: mv, rhs_mv + integer, intent(in) :: t_step + real(kind(0d0)), intent(inout) :: time_avg - real(kind(0d0)) :: t_start, t_finish - real(kind(0d0)) :: gp_sum + real(kind(0d0)) :: t_start, t_finish + real(kind(0d0)) :: gp_sum - real(kind(0d0)) :: top, bottom !< Numerator and denominator when evaluating flux limiter function - real(kind(0d0)), dimension(num_fluids) :: myalpha_rho, myalpha + real(kind(0d0)) :: top, bottom !< Numerator and denominator when evaluating flux limiter function + real(kind(0d0)), dimension(num_fluids) :: myalpha_rho, myalpha - real(kind(0d0)) :: tmp1, tmp2, tmp3, tmp4, & - c_gas, c_liquid, & - Cpbw, Cpinf, Cpinf_dot, & - myH, myHdot, rddot, alf_gas + real(kind(0d0)) :: tmp1, tmp2, tmp3, tmp4, & + c_gas, c_liquid, & + Cpbw, Cpinf, Cpinf_dot, & + myH, myHdot, rddot, alf_gas - real(kind(0d0)) :: n_tait, B_tait, angle, angle_z + real(kind(0d0)) :: n_tait, B_tait, angle, angle_z - real(kind(0d0)), dimension(nb) :: Rtmp, Vtmp - real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav - real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: nbub - integer :: ndirs + real(kind(0d0)), dimension(nb) :: Rtmp, Vtmp + real(kind(0d0)) :: myR, myV, alf, myP, myRho, R2Vav + real(kind(0d0)), dimension(0:m, 0:n, 0:p) :: nbub + integer :: ndirs - real(kind(0d0)) :: sound - real(kind(0d0)) :: start, finish - real(kind(0d0)) :: s2, const_sos, s1 + real(kind(0d0)) :: sound + real(kind(0d0)) :: start, finish + real(kind(0d0)) :: s2, const_sos, s1 - integer :: i, c, j, k, l, q, ii, id !< Generic loop iterators - integer :: term_index + integer :: i, c, j, k, l, q, ii, id !< Generic loop iterators + integer :: term_index - call nvtxStartRange("Compute_RHS") + call nvtxStartRange("Compute_RHS") - call cpu_time(t_start) - ! Association/Population of Working Variables ====================== - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, sys_size - do l = idwbuff(3)%beg, idwbuff(3)%end - do k = idwbuff(2)%beg, idwbuff(2)%end - do j = idwbuff(1)%beg, idwbuff(1)%end - q_cons_qp%vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) + call cpu_time(t_start) + ! Association/Population of Working Variables ====================== + !$acc parallel loop collapse(4) gang vector default(present) + do i = 1, sys_size + do l = idwbuff(3)%beg, idwbuff(3)%end + do k = idwbuff(2)%beg, idwbuff(2)%end + do j = idwbuff(1)%beg, idwbuff(1)%end + q_cons_qp%vf(i)%sf(j, k, l) = q_cons_vf(i)%sf(j, k, l) + end do + end do end do end do - end do - end do - ! ================================================================== + ! ================================================================== - ! Converting Conservative to Primitive Variables ================== + ! Converting Conservative to Primitive Variables ================== - if (mpp_lim .and. bubbles) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = idwbuff(3)%beg, idwbuff(3)%end - do k = idwbuff(2)%beg, idwbuff(2)%end - do j = idwbuff(1)%beg, idwbuff(1)%end - alf_sum%sf(j, k, l) = 0d0 - !$acc loop seq - do i = advxb, advxe - 1 - alf_sum%sf(j, k, l) = alf_sum%sf(j, k, l) + q_cons_qp%vf(i)%sf(j, k, l) - end do - !$acc loop seq - do i = advxb, advxe - 1 - q_cons_qp%vf(i)%sf(j, k, l) = q_cons_qp%vf(i)%sf(j, k, l)*(1.d0 - q_cons_qp%vf(alf_idx)%sf(j, k, l)) & - /alf_sum%sf(j, k, l) + if (mpp_lim .and. bubbles) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = idwbuff(3)%beg, idwbuff(3)%end + do k = idwbuff(2)%beg, idwbuff(2)%end + do j = idwbuff(1)%beg, idwbuff(1)%end + alf_sum%sf(j, k, l) = 0d0 + !$acc loop seq + do i = advxb, advxe - 1 + alf_sum%sf(j, k, l) = alf_sum%sf(j, k, l) + q_cons_qp%vf(i)%sf(j, k, l) + end do + !$acc loop seq + do i = advxb, advxe - 1 + q_cons_qp%vf(i)%sf(j, k, l) = q_cons_qp%vf(i)%sf(j, k, l)*(1.d0 - q_cons_qp%vf(alf_idx)%sf(j, k, l)) & + /alf_sum%sf(j, k, l) + end do + end do end do end do - end do - end do - end if + end if - call nvtxStartRange("RHS-CONVERT") - call s_convert_conservative_to_primitive_variables( & - q_cons_qp%vf, & - q_prim_qp%vf, & - idwint, & - gm_alpha_qp%vf) - call nvtxEndRange + call nvtxStartRange("RHS-CONVERT") + call s_convert_conservative_to_primitive_variables( & + q_cons_qp%vf, & + q_prim_qp%vf, & + idwint, & + gm_alpha_qp%vf) + call nvtxEndRange - call nvtxStartRange("RHS-MPI") - call s_populate_variables_buffers(q_prim_qp%vf, pb, mv) - call nvtxEndRange + call nvtxStartRange("RHS-MPI") + call s_populate_variables_buffers(q_prim_qp%vf, pb, mv) + call nvtxEndRange - if (cfl_dt) then - if (mytime >= t_stop) return - else - if (t_step == t_step_stop) return - end if + if (cfl_dt) then + if (mytime >= t_stop) return + else + if (t_step == t_step_stop) return + end if - ! ================================================================== + ! ================================================================== - if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb, rhs_pb, mv, rhs_mv, idwbuff(1), idwbuff(2), idwbuff(3), nbub) + if (qbmm) call s_mom_inv(q_cons_qp%vf, q_prim_qp%vf, mom_sp, mom_3d, pb, rhs_pb, mv, rhs_mv, idwbuff(1), idwbuff(2), idwbuff(3), nbub) - call nvtxStartRange("Viscous") - if (viscous) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n, & - qL_prim, & - qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n, & - qR_prim, & - q_prim_qp, & - dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, & - idwbuff(1), idwbuff(2), idwbuff(3)) - call nvtxEndRange + call nvtxStartRange("Viscous") + if (viscous) call s_get_viscous(qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n, & + qL_prim, & + qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n, & + qR_prim, & + q_prim_qp, & + dq_prim_dx_qp, dq_prim_dy_qp, dq_prim_dz_qp, & + idwbuff(1), idwbuff(2), idwbuff(3)) + call nvtxEndRange - call nvtxStartRange("Surface_Tension") - if (surface_tension) call s_get_capilary(q_prim_qp%vf) - call nvtxEndRange - ! Dimensional Splitting Loop ======================================= + call nvtxStartRange("Surface_Tension") + if (surface_tension) call s_get_capilary(q_prim_qp%vf) + call nvtxEndRange + ! Dimensional Splitting Loop ======================================= - do id = 1, num_dims + do id = 1, num_dims - ! Reconstructing Primitive/Conservative Variables =============== + ! Reconstructing Primitive/Conservative Variables =============== - call nvtxStartRange("RHS-WENO") + call nvtxStartRange("RHS-WENO") - if (.not. surface_tension) then - ! Reconstruct densitiess - iv%beg = 1; iv%end = sys_size - call s_reconstruct_cell_boundary_values( & - q_prim_qp%vf(1:sys_size), & - qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - id) - else - iv%beg = 1; iv%end = E_idx - 1 - call s_reconstruct_cell_boundary_values( & - q_prim_qp%vf(iv%beg:iv%end), & - qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - id) - - iv%beg = E_idx; iv%end = E_idx - call s_reconstruct_cell_boundary_values_first_order( & - q_prim_qp%vf(E_idx), & - qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - id) - - iv%beg = E_idx + 1; iv%end = sys_size - call s_reconstruct_cell_boundary_values( & - q_prim_qp%vf(iv%beg:iv%end), & - qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - id) - end if + if (.not. surface_tension) then + ! Reconstruct densitiess + iv%beg = 1; iv%end = sys_size + call s_reconstruct_cell_boundary_values( & + q_prim_qp%vf(1:sys_size), & + qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + id) + else + iv%beg = 1; iv%end = E_idx - 1 + call s_reconstruct_cell_boundary_values( & + q_prim_qp%vf(iv%beg:iv%end), & + qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + id) + + iv%beg = E_idx; iv%end = E_idx + call s_reconstruct_cell_boundary_values_first_order( & + q_prim_qp%vf(E_idx), & + qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + id) + + iv%beg = E_idx + 1; iv%end = sys_size + call s_reconstruct_cell_boundary_values( & + q_prim_qp%vf(iv%beg:iv%end), & + qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + id) + end if - ! Reconstruct viscous derivatives for viscosity - if (weno_Re_flux) then - iv%beg = momxb; iv%end = momxe - call s_reconstruct_cell_boundary_values_visc_deriv( & - dq_prim_dx_qp(1)%vf(iv%beg:iv%end), & - dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, & - dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, & - id, dqL_prim_dx_n(id)%vf(iv%beg:iv%end), dqR_prim_dx_n(id)%vf(iv%beg:iv%end), & - idwbuff(1), idwbuff(2), idwbuff(3)) - if (n > 0) then - call s_reconstruct_cell_boundary_values_visc_deriv( & - dq_prim_dy_qp(1)%vf(iv%beg:iv%end), & - dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, & - dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, & - id, dqL_prim_dy_n(id)%vf(iv%beg:iv%end), dqR_prim_dy_n(id)%vf(iv%beg:iv%end), & - idwbuff(1), idwbuff(2), idwbuff(3)) - if (p > 0) then + ! Reconstruct viscous derivatives for viscosity + if (weno_Re_flux) then + iv%beg = momxb; iv%end = momxe call s_reconstruct_cell_boundary_values_visc_deriv( & - dq_prim_dz_qp(1)%vf(iv%beg:iv%end), & + dq_prim_dx_qp(1)%vf(iv%beg:iv%end), & dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, & dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, & - id, dqL_prim_dz_n(id)%vf(iv%beg:iv%end), dqR_prim_dz_n(id)%vf(iv%beg:iv%end), & + id, dqL_prim_dx_n(id)%vf(iv%beg:iv%end), dqR_prim_dx_n(id)%vf(iv%beg:iv%end), & idwbuff(1), idwbuff(2), idwbuff(3)) + if (n > 0) then + call s_reconstruct_cell_boundary_values_visc_deriv( & + dq_prim_dy_qp(1)%vf(iv%beg:iv%end), & + dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, & + dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, & + id, dqL_prim_dy_n(id)%vf(iv%beg:iv%end), dqR_prim_dy_n(id)%vf(iv%beg:iv%end), & + idwbuff(1), idwbuff(2), idwbuff(3)) + if (p > 0) then + call s_reconstruct_cell_boundary_values_visc_deriv( & + dq_prim_dz_qp(1)%vf(iv%beg:iv%end), & + dqL_rsx_vf, dqL_rsy_vf, dqL_rsz_vf, & + dqR_rsx_vf, dqR_rsy_vf, dqR_rsz_vf, & + id, dqL_prim_dz_n(id)%vf(iv%beg:iv%end), dqR_prim_dz_n(id)%vf(iv%beg:iv%end), & + idwbuff(1), idwbuff(2), idwbuff(3)) + end if + end if end if - end if - end if - call nvtxEndRange ! WENO + call nvtxEndRange ! WENO - ! Configuring Coordinate Direction Indexes ====================== - if (id == 1) then - irx%beg = -1; iry%beg = 0; irz%beg = 0 - elseif (id == 2) then - irx%beg = 0; iry%beg = -1; irz%beg = 0 - else - irx%beg = 0; iry%beg = 0; irz%beg = -1 - end if - irx%end = m; iry%end = n; irz%end = p - ! =============================================================== - - call nvtxStartRange("RHS_riemann_solver") - - ! Computing Riemann Solver Flux and Source Flux ================= - - call s_riemann_solver(qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & - dqR_prim_dx_n(id)%vf, & - dqR_prim_dy_n(id)%vf, & - dqR_prim_dz_n(id)%vf, & - qR_prim(id)%vf, & - qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & - dqL_prim_dx_n(id)%vf, & - dqL_prim_dy_n(id)%vf, & - dqL_prim_dz_n(id)%vf, & - qL_prim(id)%vf, & - q_prim_qp%vf, & - flux_n(id)%vf, & - flux_src_n(id)%vf, & - flux_gsrc_n(id)%vf, & - id, irx, iry, irz) - call nvtxEndRange - - ! Additional physics and source terms ============================== - - ! RHS addition for advection source - call nvtxStartRange("RHS_advection_source") - call s_compute_advection_source_term(id, & - rhs_vf, & - q_cons_qp, & - q_prim_qp, & - flux_src_n(id)) - call nvtxEndRange - - ! RHS additions for hypoelasticity - call nvtxStartRange("RHS_Hypoelasticity") - if (hypoelasticity) call s_compute_hypoelastic_rhs(id, & - q_prim_qp%vf, & - rhs_vf) - call nvtxEndRange - - ! RHS additions for viscosity - call nvtxStartRange("RHS_add_phys") - if (viscous .or. surface_tension) then - call s_compute_additional_physics_rhs(id, & + ! Configuring Coordinate Direction Indexes ====================== + if (id == 1) then + irx%beg = -1; iry%beg = 0; irz%beg = 0 + elseif (id == 2) then + irx%beg = 0; iry%beg = -1; irz%beg = 0 + else + irx%beg = 0; iry%beg = 0; irz%beg = -1 + end if + irx%end = m; iry%end = n; irz%end = p + ! =============================================================== + + call nvtxStartRange("RHS_riemann_solver") + + ! Computing Riemann Solver Flux and Source Flux ================= + + call s_riemann_solver(qR_rsx_vf, qR_rsy_vf, qR_rsz_vf, & + dqR_prim_dx_n(id)%vf, & + dqR_prim_dy_n(id)%vf, & + dqR_prim_dz_n(id)%vf, & + qR_prim(id)%vf, & + qL_rsx_vf, qL_rsy_vf, qL_rsz_vf, & + dqL_prim_dx_n(id)%vf, & + dqL_prim_dy_n(id)%vf, & + dqL_prim_dz_n(id)%vf, & + qL_prim(id)%vf, & + q_prim_qp%vf, & + flux_n(id)%vf, & + flux_src_n(id)%vf, & + flux_gsrc_n(id)%vf, & + id, irx, iry, irz) + call nvtxEndRange + + ! Additional physics and source terms ============================== + + ! RHS addition for advection source + call nvtxStartRange("RHS_advection_source") + call s_compute_advection_source_term(id, & + rhs_vf, & + q_cons_qp, & + q_prim_qp, & + flux_src_n(id)) + call nvtxEndRange + + ! RHS additions for hypoelasticity + call nvtxStartRange("RHS_Hypoelasticity") + if (hypoelasticity) call s_compute_hypoelastic_rhs(id, & + q_prim_qp%vf, & + rhs_vf) + call nvtxEndRange + + ! RHS additions for viscosity + call nvtxStartRange("RHS_add_phys") + if (viscous .or. surface_tension) then + call s_compute_additional_physics_rhs(id, & + q_prim_qp%vf, & + rhs_vf, & + flux_src_n(id)%vf, & + dq_prim_dx_qp(1)%vf, & + dq_prim_dy_qp(1)%vf, & + dq_prim_dz_qp(1)%vf) + end if + call nvtxEndRange + + ! RHS additions for sub-grid bubbles + call nvtxStartRange("RHS_bubbles") + if (bubbles) call s_compute_bubbles_rhs(id, & + q_prim_qp%vf) + call nvtxEndRange + + ! RHS additions for qbmm bubbles + call nvtxStartRange("RHS_qbmm") + if (qbmm) call s_compute_qbmm_rhs(id, & + q_cons_qp%vf, & q_prim_qp%vf, & rhs_vf, & - flux_src_n(id)%vf, & - dq_prim_dx_qp(1)%vf, & - dq_prim_dy_qp(1)%vf, & - dq_prim_dz_qp(1)%vf) - end if - call nvtxEndRange - - ! RHS additions for sub-grid bubbles - call nvtxStartRange("RHS_bubbles") - if (bubbles) call s_compute_bubbles_rhs(id, & - q_prim_qp%vf) - call nvtxEndRange - - ! RHS additions for qbmm bubbles - call nvtxStartRange("RHS_qbmm") - if (qbmm) call s_compute_qbmm_rhs(id, & - q_cons_qp%vf, & - q_prim_qp%vf, & - rhs_vf, & - flux_n(id)%vf, & - pb, & - rhs_pb, & - mv, & - rhs_mv) - call nvtxEndRange - ! END: Additional physics and source terms ========================= + flux_n(id)%vf, & + pb, & + rhs_pb, & + mv, & + rhs_mv) + call nvtxEndRange + ! END: Additional physics and source terms ========================= - end do - ! END: Dimensional Splitting Loop ================================= - - #:if chemistry - call nvtxStartRange("RHS_Chem_Advection") - call s_compute_chemistry_advection_flux(flux_n, rhs_vf) - call nvtxEndRange - #:endif - - if (ib) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - if (ib_markers%sf(j, k, l) /= 0) then - do i = 1, sys_size - rhs_vf(i)%sf(j, k, l) = 0d0 - end do - end if - end do end do - end do - end if - - ! Additional Physics and Source Temrs ================================== - ! Additions for acoustic_source - call nvtxStartRange("RHS_acoustic_src") - if (acoustic_source) call s_acoustic_src_calculations(q_cons_qp%vf(1:sys_size), & - q_prim_qp%vf(1:sys_size), & - t_step, & - rhs_vf) - call nvtxEndRange - - ! Add bubles source term - call nvtxStartRange("RHS_bubbles") - if (bubbles .and. (.not. adap_dt) .and. (.not. qbmm)) call s_compute_bubble_source( & - q_cons_qp%vf(1:sys_size), & - q_prim_qp%vf(1:sys_size), & - t_step, & - rhs_vf) - call nvtxEndRange - - #:if chemistry - if (chem_params%reactions) then - call nvtxStartRange("RHS_Chem_Reactions") - call s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp%vf, q_prim_qp%vf) - call nvtxEndRange - end if - #:endif + ! END: Dimensional Splitting Loop ================================= - ! END: Additional pphysics and source terms ============================ + #:if chemistry + call nvtxStartRange("RHS_Chem_Advection") + call s_compute_chemistry_advection_flux(flux_n, rhs_vf) + call nvtxEndRange + #:endif - if (run_time_info .or. probe_wrt .or. ib) then - !$acc parallel loop collapse(4) gang vector default(present) - do i = 1, sys_size - do l = idwbuff(3)%beg, idwbuff(3)%end - do k = idwbuff(2)%beg, idwbuff(2)%end - do j = idwbuff(1)%beg, idwbuff(1)%end - q_prim_vf(i)%sf(j, k, l) = q_prim_qp%vf(i)%sf(j, k, l) + if (ib) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + if (ib_markers%sf(j, k, l) /= 0) then + do i = 1, sys_size + rhs_vf(i)%sf(j, k, l) = 0d0 + end do + end if + end do end do end do - end do - end do - end if - - call cpu_time(t_finish) + end if - if (t_step >= 4) then - time_avg = (abs(t_finish - t_start) + (t_step - 4)*time_avg)/(t_step - 3) - else - time_avg = 0d0 - end if - ! ================================================================== + ! Additional Physics and Source Temrs ================================== + ! Additions for acoustic_source + call nvtxStartRange("RHS_acoustic_src") + if (acoustic_source) call s_acoustic_src_calculations(q_cons_qp%vf(1:sys_size), & + q_prim_qp%vf(1:sys_size), & + t_step, & + rhs_vf) + call nvtxEndRange - call nvtxEndRange + ! Add bubles source term + call nvtxStartRange("RHS_bubbles") + if (bubbles .and. (.not. adap_dt) .and. (.not. qbmm)) call s_compute_bubble_source( & + q_cons_qp%vf(1:sys_size), & + q_prim_qp%vf(1:sys_size), & + t_step, & + rhs_vf) + call nvtxEndRange - end subroutine s_compute_rhs + #:if chemistry + if (chem_params%reactions) then + call nvtxStartRange("RHS_Chem_Reactions") + call s_compute_chemistry_reaction_flux(rhs_vf, q_cons_qp%vf, q_prim_qp%vf) + call nvtxEndRange + end if + #:endif - subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf) + ! END: Additional pphysics and source terms ============================ - integer, intent(in) :: idir - type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - type(vector_field), intent(inout) :: q_cons_vf - type(vector_field), intent(inout) :: q_prim_vf - type(vector_field), intent(inout) :: flux_src_n_vf + if (run_time_info .or. probe_wrt .or. ib) then + !$acc parallel loop collapse(4) gang vector default(present) + do i = 1, sys_size + do l = idwbuff(3)%beg, idwbuff(3)%end + do k = idwbuff(2)%beg, idwbuff(2)%end + do j = idwbuff(1)%beg, idwbuff(1)%end + q_prim_vf(i)%sf(j, k, l) = q_prim_qp%vf(i)%sf(j, k, l) + end do + end do + end do + end do + end if - integer :: i, j, k, l, q + call cpu_time(t_finish) - if (alt_soundspeed) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - blkmod1(j, k, l) = ((gammas(1) + 1d0)*q_prim_vf%vf(E_idx)%sf(j, k, l) + & - pi_infs(1))/gammas(1) - blkmod2(j, k, l) = ((gammas(2) + 1d0)*q_prim_vf%vf(E_idx)%sf(j, k, l) + & - pi_infs(2))/gammas(2) - alpha1(j, k, l) = q_cons_vf%vf(advxb)%sf(j, k, l) + if (t_step >= 4) then + time_avg = (abs(t_finish - t_start) + (t_step - 4)*time_avg)/(t_step - 3) + else + time_avg = 0d0 + end if + ! ================================================================== - if (bubbles) then - alpha2(j, k, l) = q_cons_vf%vf(alf_idx - 1)%sf(j, k, l) - else - alpha2(j, k, l) = q_cons_vf%vf(advxe)%sf(j, k, l) - end if + call nvtxEndRange - Kterm(j, k, l) = alpha1(j, k, l)*alpha2(j, k, l)*(blkmod2(j, k, l) - blkmod1(j, k, l))/ & - (alpha1(j, k, l)*blkmod2(j, k, l) + alpha2(j, k, l)*blkmod1(j, k, l)) - end do - end do - end do - end if + end subroutine s_compute_rhs - if (idir == 1) then + subroutine s_compute_advection_source_term(idir, rhs_vf, q_cons_vf, q_prim_vf, flux_src_n_vf) - if (bc_x%beg <= -5 .and. bc_x%beg >= -13) then - call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & - flux_src_n(idir)%vf, idir, -1, irx, iry, irz) - end if + integer, intent(in) :: idir + type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf + type(vector_field), intent(inout) :: q_cons_vf + type(vector_field), intent(inout) :: q_prim_vf + type(vector_field), intent(inout) :: flux_src_n_vf - if (bc_x%end <= -5 .and. bc_x%end >= -13) then - call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & - flux_src_n(idir)%vf, idir, 1, irx, iry, irz) - end if + integer :: i, j, k, l, q - !$acc parallel loop collapse(4) gang vector default(present) - do j = 1, sys_size - do q = 0, p - do l = 0, n - do k = 0, m - rhs_vf(j)%sf(k, l, q) = 1d0/dx(k)* & - (flux_n(1)%vf(j)%sf(k - 1, l, q) & - - flux_n(1)%vf(j)%sf(k, l, q)) - end do - end do - end do - end do + if (alt_soundspeed) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + blkmod1(j, k, l) = ((gammas(1) + 1d0)*q_prim_vf%vf(E_idx)%sf(j, k, l) + & + pi_infs(1))/gammas(1) + blkmod2(j, k, l) = ((gammas(2) + 1d0)*q_prim_vf%vf(E_idx)%sf(j, k, l) + & + pi_infs(2))/gammas(2) + alpha1(j, k, l) = q_cons_vf%vf(advxb)%sf(j, k, l) + + if (bubbles) then + alpha2(j, k, l) = q_cons_vf%vf(alf_idx - 1)%sf(j, k, l) + else + alpha2(j, k, l) = q_cons_vf%vf(advxe)%sf(j, k, l) + end if - if (model_eqns == 3) then - !$acc parallel loop collapse(4) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - do i = 1, num_fluids - rhs_vf(i + intxb - 1)%sf(j, k, l) = & - rhs_vf(i + intxb - 1)%sf(j, k, l) - 1d0/dx(j)* & - q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & - q_prim_vf%vf(E_idx)%sf(j, k, l)* & - (flux_src_n(1)%vf(advxb)%sf(j, k, l) - & - flux_src_n(1)%vf(advxb)%sf(j - 1, k, l)) + Kterm(j, k, l) = alpha1(j, k, l)*alpha2(j, k, l)*(blkmod2(j, k, l) - blkmod1(j, k, l))/ & + (alpha1(j, k, l)*blkmod2(j, k, l) + alpha2(j, k, l)*blkmod1(j, k, l)) end do end do end do - end do - end if + end if - if (riemann_solver == 1) then - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe - do q = 0, p - do l = 0, n - do k = 0, m - rhs_vf(j)%sf(k, l, q) = & - rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - q_prim_vf%vf(contxe + idir)%sf(k, l, q)* & - (flux_src_n(1)%vf(j)%sf(k - 1, l, q) & - - flux_src_n(1)%vf(j)%sf(k, l, q)) + if (idir == 1) then + + if (bc_x%beg <= -5 .and. bc_x%beg >= -13) then + call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & + flux_src_n(idir)%vf, idir, -1, irx, iry, irz) + end if + + if (bc_x%end <= -5 .and. bc_x%end >= -13) then + call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & + flux_src_n(idir)%vf, idir, 1, irx, iry, irz) + end if + + !$acc parallel loop collapse(4) gang vector default(present) + do j = 1, sys_size + do q = 0, p + do l = 0, n + do k = 0, m + rhs_vf(j)%sf(k, l, q) = 1d0/dx(k)* & + (flux_n(1)%vf(j)%sf(k - 1, l, q) & + - flux_n(1)%vf(j)%sf(k, l, q)) + end do end do end do end do - end do - else - if (alt_soundspeed) then - do j = advxb, advxe - if ((j == advxe) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) - do q = 0, p - do l = 0, n - do k = 0, m - rhs_vf(j)%sf(k, l, q) = & - rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - (q_cons_vf%vf(j)%sf(k, l, q) - Kterm(k, l, q))* & - (flux_src_n(1)%vf(j)%sf(k, l, q) & - - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) + + if (model_eqns == 3) then + !$acc parallel loop collapse(4) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + do i = 1, num_fluids + rhs_vf(i + intxb - 1)%sf(j, k, l) = & + rhs_vf(i + intxb - 1)%sf(j, k, l) - 1d0/dx(j)* & + q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & + q_prim_vf%vf(E_idx)%sf(j, k, l)* & + (flux_src_n(1)%vf(advxb)%sf(j, k, l) - & + flux_src_n(1)%vf(advxb)%sf(j - 1, k, l)) end do end do end do - else if ((j == advxb) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) + end do + end if + + if (riemann_solver == 1) then + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe do q = 0, p do l = 0, n do k = 0, m rhs_vf(j)%sf(k, l, q) = & rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - (q_cons_vf%vf(j)%sf(k, l, q) + Kterm(k, l, q))* & - (flux_src_n(1)%vf(j)%sf(k, l, q) & - - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) + q_prim_vf%vf(contxe + idir)%sf(k, l, q)* & + (flux_src_n(1)%vf(j)%sf(k - 1, l, q) & + - flux_src_n(1)%vf(j)%sf(k, l, q)) end do end do end do - end if - end do - else - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe - do q = 0, p - do l = 0, n - do k = 0, m - rhs_vf(j)%sf(k, l, q) = & - rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & - q_cons_vf%vf(j)%sf(k, l, q)* & - (flux_src_n(1)%vf(j)%sf(k, l, q) & - - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) + end do + else + if (alt_soundspeed) then + do j = advxb, advxe + if ((j == advxe) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do q = 0, p + do l = 0, n + do k = 0, m + rhs_vf(j)%sf(k, l, q) = & + rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & + (q_cons_vf%vf(j)%sf(k, l, q) - Kterm(k, l, q))* & + (flux_src_n(1)%vf(j)%sf(k, l, q) & + - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) + end do + end do + end do + else if ((j == advxb) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do q = 0, p + do l = 0, n + do k = 0, m + rhs_vf(j)%sf(k, l, q) = & + rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & + (q_cons_vf%vf(j)%sf(k, l, q) + Kterm(k, l, q))* & + (flux_src_n(1)%vf(j)%sf(k, l, q) & + - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) + end do + end do + end do + end if + end do + else + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe + do q = 0, p + do l = 0, n + do k = 0, m + rhs_vf(j)%sf(k, l, q) = & + rhs_vf(j)%sf(k, l, q) + 1d0/dx(k)* & + q_cons_vf%vf(j)%sf(k, l, q)* & + (flux_src_n(1)%vf(j)%sf(k, l, q) & + - flux_src_n(1)%vf(j)%sf(k - 1, l, q)) + end do + end do end do end do - end do - end do - end if - end if - - elseif (idir == 2) then - ! RHS Contribution in y-direction =============================== - ! Applying the Riemann fluxes - - if (bc_y%beg <= -5 .and. bc_y%beg >= -13) then - call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & - flux_src_n(idir)%vf, idir, -1, irx, iry, irz) - end if + end if + end if - if (bc_y%end <= -5 .and. bc_y%end >= -13) then - call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & - flux_src_n(idir)%vf, idir, 1, irx, iry, irz) - end if + elseif (idir == 2) then + ! RHS Contribution in y-direction =============================== + ! Applying the Riemann fluxes - !$acc parallel loop collapse(4) gang vector default(present) - do j = 1, sys_size - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - (flux_n(2)%vf(j)%sf(q, k - 1, l) & - - flux_n(2)%vf(j)%sf(q, k, l)) - end do - end do - end do - end do + if (bc_y%beg <= -5 .and. bc_y%beg >= -13) then + call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & + flux_src_n(idir)%vf, idir, -1, irx, iry, irz) + end if - if (model_eqns == 3) then - !$acc parallel loop collapse(4) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - do i = 1, num_fluids - rhs_vf(i + intxb - 1)%sf(j, k, l) = & - rhs_vf(i + intxb - 1)%sf(j, k, l) - 1d0/dy(k)* & - q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & - q_prim_vf%vf(E_idx)%sf(j, k, l)* & - (flux_src_n(2)%vf(advxb)%sf(j, k, l) - & - flux_src_n(2)%vf(advxb)%sf(j, k - 1, l)) - end do - end do - end do - end do + if (bc_y%end <= -5 .and. bc_y%end >= -13) then + call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & + flux_src_n(idir)%vf, idir, 1, irx, iry, irz) + end if - if (cyl_coord) then !$acc parallel loop collapse(4) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - do i = 1, num_fluids - rhs_vf(i + intxb - 1)%sf(j, k, l) = & - rhs_vf(i + intxb - 1)%sf(j, k, l) - 5d-1/y_cc(k)* & - q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & - q_prim_vf%vf(E_idx)%sf(j, k, l)* & - (flux_src_n(2)%vf(advxb)%sf(j, k, l) + & - flux_src_n(2)%vf(advxb)%sf(j, k - 1, l)) + do j = 1, sys_size + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + (flux_n(2)%vf(j)%sf(q, k - 1, l) & + - flux_n(2)%vf(j)%sf(q, k, l)) end do end do end do end do - end if - end if - if (cyl_coord) then - !$acc parallel loop collapse(4) gang vector default(present) - do j = 1, sys_size - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) - 5d-1/y_cc(k)* & - (flux_gsrc_n(2)%vf(j)%sf(q, k - 1, l) & - + flux_gsrc_n(2)%vf(j)%sf(q, k, l)) + if (model_eqns == 3) then + !$acc parallel loop collapse(4) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + do i = 1, num_fluids + rhs_vf(i + intxb - 1)%sf(j, k, l) = & + rhs_vf(i + intxb - 1)%sf(j, k, l) - 1d0/dy(k)* & + q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & + q_prim_vf%vf(E_idx)%sf(j, k, l)* & + (flux_src_n(2)%vf(advxb)%sf(j, k, l) - & + flux_src_n(2)%vf(advxb)%sf(j, k - 1, l)) + end do + end do end do end do - end do - end do - end if - if (riemann_solver == 1) then - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - q_prim_vf%vf(contxe + idir)%sf(q, k, l)* & - (flux_src_n(2)%vf(j)%sf(q, k - 1, l) & - - flux_src_n(2)%vf(j)%sf(q, k, l)) + if (cyl_coord) then + !$acc parallel loop collapse(4) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + do i = 1, num_fluids + rhs_vf(i + intxb - 1)%sf(j, k, l) = & + rhs_vf(i + intxb - 1)%sf(j, k, l) - 5d-1/y_cc(k)* & + q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & + q_prim_vf%vf(E_idx)%sf(j, k, l)* & + (flux_src_n(2)%vf(advxb)%sf(j, k, l) + & + flux_src_n(2)%vf(advxb)%sf(j, k - 1, l)) + end do + end do + end do end do - end do - end do - end do - else + end if + end if - if (alt_soundspeed) then - do j = advxb, advxe - if ((j == advxe) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) + if (cyl_coord) then + !$acc parallel loop collapse(4) gang vector default(present) + do j = 1, sys_size do l = 0, p do k = 0, n do q = 0, m rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - (q_cons_vf%vf(j)%sf(q, k, l) - Kterm(q, k, l))* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + rhs_vf(j)%sf(q, k, l) - 5d-1/y_cc(k)* & + (flux_gsrc_n(2)%vf(j)%sf(q, k - 1, l) & + + flux_gsrc_n(2)%vf(j)%sf(q, k, l)) end do end do end do - if (cyl_coord) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) - & - (Kterm(q, k, l)/2d0/y_cc(k))* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) - end do - end do - end do - end if - else if ((j == advxb) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) + end do + end if + + if (riemann_solver == 1) then + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe do l = 0, p do k = 0, n do q = 0, m rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - (q_cons_vf%vf(j)%sf(q, k, l) + Kterm(q, k, l))* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + q_prim_vf%vf(contxe + idir)%sf(q, k, l)* & + (flux_src_n(2)%vf(j)%sf(q, k - 1, l) & + - flux_src_n(2)%vf(j)%sf(q, k, l)) end do end do end do - if (cyl_coord) then - !$acc parallel loop collapse(3) gang vector default(present) + end do + else + + if (alt_soundspeed) then + do j = advxb, advxe + if ((j == advxe) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + (q_cons_vf%vf(j)%sf(q, k, l) - Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + if (cyl_coord) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) - & + (Kterm(q, k, l)/2d0/y_cc(k))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + end if + else if ((j == advxb) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + (q_cons_vf%vf(j)%sf(q, k, l) + Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + if (cyl_coord) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + & + (Kterm(q, k, l)/2d0/y_cc(k))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + end if + end if + end do + else + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe do l = 0, p do k = 0, n do q = 0, m rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + & - (Kterm(q, k, l)/2d0/y_cc(k))* & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + q_cons_vf%vf(j)%sf(q, k, l)* & (flux_src_n(2)%vf(j)%sf(q, k, l) & - + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do - end if - end if - end do - else - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - q_cons_vf%vf(j)%sf(q, k, l)* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) - end do end do - end do - end do - end if - end if - - elseif (idir == 3) then - ! RHS Contribution in z-direction =============================== + end if + end if - ! Applying the Riemann fluxes + elseif (idir == 3) then + ! RHS Contribution in z-direction =============================== - if (bc_z%beg <= -5 .and. bc_z%beg >= -13) then - call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & - flux_src_n(idir)%vf, idir, -1, irx, iry, irz) - end if + ! Applying the Riemann fluxes - if (bc_z%end <= -5 .and. bc_z%end >= -13) then - call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & - flux_src_n(idir)%vf, idir, 1, irx, iry, irz) - end if + if (bc_z%beg <= -5 .and. bc_z%beg >= -13) then + call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & + flux_src_n(idir)%vf, idir, -1, irx, iry, irz) + end if - if (grid_geometry == 3) then ! Cylindrical Coordinates - !$acc parallel loop collapse(4) gang vector default(present) - do j = 1, sys_size - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)/y_cc(q)* & - q_prim_vf%vf(contxe + idir)%sf(l, q, k)* & - (flux_n(3)%vf(j)%sf(l, q, k - 1) & - - flux_n(3)%vf(j)%sf(l, q, k)) - end do - end do - end do - end do + if (bc_z%end <= -5 .and. bc_z%end >= -13) then + call s_cbc(q_prim_vf%vf, flux_n(idir)%vf, & + flux_src_n(idir)%vf, idir, 1, irx, iry, irz) + end if - !$acc parallel loop collapse(4) gang vector default(present) - do j = 1, sys_size - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) - 5d-1/y_cc(q)* & - (flux_gsrc_n(3)%vf(j)%sf(l, q, k - 1) & - - flux_gsrc_n(3)%vf(j)%sf(l, q, k)) + if (grid_geometry == 3) then ! Cylindrical Coordinates + !$acc parallel loop collapse(4) gang vector default(present) + do j = 1, sys_size + do k = 0, p + do q = 0, n + do l = 0, m + rhs_vf(j)%sf(l, q, k) = & + rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)/y_cc(q)* & + q_prim_vf%vf(contxe + idir)%sf(l, q, k)* & + (flux_n(3)%vf(j)%sf(l, q, k - 1) & + - flux_n(3)%vf(j)%sf(l, q, k)) + end do + end do end do end do - end do - end do - else ! Cartesian Coordinates - !$acc parallel loop collapse(4) gang vector default(present) - do j = 1, sys_size - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - (flux_n(3)%vf(j)%sf(l, q, k - 1) & - - flux_n(3)%vf(j)%sf(l, q, k)) + !$acc parallel loop collapse(4) gang vector default(present) + do j = 1, sys_size + do k = 0, p + do q = 0, n + do l = 0, m + rhs_vf(j)%sf(l, q, k) = & + rhs_vf(j)%sf(l, q, k) - 5d-1/y_cc(q)* & + (flux_gsrc_n(3)%vf(j)%sf(l, q, k - 1) & + - flux_gsrc_n(3)%vf(j)%sf(l, q, k)) + end do + end do end do end do - end do - end do - end if - if (model_eqns == 3) then - !$acc parallel loop collapse(4) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - do i = 1, num_fluids - rhs_vf(i + intxb - 1)%sf(j, k, l) = & - rhs_vf(i + intxb - 1)%sf(j, k, l) - 1d0/dz(l)* & - q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & - q_prim_vf%vf(E_idx)%sf(j, k, l)* & - (flux_src_n(3)%vf(advxb)%sf(j, k, l) - & - flux_src_n(3)%vf(advxb)%sf(j, k, l - 1)) + else ! Cartesian Coordinates + !$acc parallel loop collapse(4) gang vector default(present) + do j = 1, sys_size + do k = 0, p + do q = 0, n + do l = 0, m + rhs_vf(j)%sf(l, q, k) = & + rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & + (flux_n(3)%vf(j)%sf(l, q, k - 1) & + - flux_n(3)%vf(j)%sf(l, q, k)) + end do + end do end do end do - end do - end do - end if + end if - if (grid_geometry == 3) then - if (riemann_solver == 1) then - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe + if (model_eqns == 3) then + !$acc parallel loop collapse(4) gang vector default(present) do l = 0, p do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - q_prim_vf%vf(contxe + idir)%sf(q, k, l)* & - (flux_src_n(2)%vf(j)%sf(q, k - 1, l) & - - flux_src_n(2)%vf(j)%sf(q, k, l)) + do j = 0, m + do i = 1, num_fluids + rhs_vf(i + intxb - 1)%sf(j, k, l) = & + rhs_vf(i + intxb - 1)%sf(j, k, l) - 1d0/dz(l)* & + q_cons_vf%vf(i + advxb - 1)%sf(j, k, l)* & + q_prim_vf%vf(E_idx)%sf(j, k, l)* & + (flux_src_n(3)%vf(advxb)%sf(j, k, l) - & + flux_src_n(3)%vf(advxb)%sf(j, k, l - 1)) + end do end do end do end do - end do - else + end if - if (alt_soundspeed) then - do j = advxb, advxe - if ((j == advxe) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) + if (grid_geometry == 3) then + if (riemann_solver == 1) then + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe do l = 0, p do k = 0, n do q = 0, m rhs_vf(j)%sf(q, k, l) = & rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - (q_cons_vf%vf(j)%sf(q, k, l) - Kterm(q, k, l))* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + q_prim_vf%vf(contxe + idir)%sf(q, k, l)* & + (flux_src_n(2)%vf(j)%sf(q, k - 1, l) & + - flux_src_n(2)%vf(j)%sf(q, k, l)) end do end do end do - if (cyl_coord) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) - & - (Kterm(q, k, l)/2d0/y_cc(k))* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + else + + if (alt_soundspeed) then + do j = advxb, advxe + if ((j == advxe) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + (q_cons_vf%vf(j)%sf(q, k, l) - Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do end do end do - end do - end if - else if ((j == advxb) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - (q_cons_vf%vf(j)%sf(q, k, l) + Kterm(q, k, l))* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + if (cyl_coord) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) - & + (Kterm(q, k, l)/2d0/y_cc(k))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + end if + else if ((j == advxb) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + (q_cons_vf%vf(j)%sf(q, k, l) + Kterm(q, k, l))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do end do - end do + if (cyl_coord) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do q = 0, m + rhs_vf(j)%sf(q, k, l) = & + rhs_vf(j)%sf(q, k, l) + & + (Kterm(q, k, l)/2d0/y_cc(k))* & + (flux_src_n(2)%vf(j)%sf(q, k, l) & + + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + end do + end do + end do + end if + end if end do - if (cyl_coord) then - !$acc parallel loop collapse(3) gang vector default(present) + else + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe do l = 0, p do k = 0, n do q = 0, m rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + & - (Kterm(q, k, l)/2d0/y_cc(k))* & + rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & + q_cons_vf%vf(j)%sf(q, k, l)* & (flux_src_n(2)%vf(j)%sf(q, k, l) & - + flux_src_n(2)%vf(j)%sf(q, k - 1, l)) + - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) end do end do end do - end if + end do end if - end do + end if else - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe - do l = 0, p - do k = 0, n - do q = 0, m - rhs_vf(j)%sf(q, k, l) = & - rhs_vf(j)%sf(q, k, l) + 1d0/dy(k)* & - q_cons_vf%vf(j)%sf(q, k, l)* & - (flux_src_n(2)%vf(j)%sf(q, k, l) & - - flux_src_n(2)%vf(j)%sf(q, k - 1, l)) - end do - end do - end do - end do - end if - end if - else - if (riemann_solver == 1) then - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - q_prim_vf%vf(contxe + idir)%sf(l, q, k)* & - (flux_src_n(3)%vf(j)%sf(l, q, k - 1) & - - flux_src_n(3)%vf(j)%sf(l, q, k)) - end do - end do - end do - end do - else - if (alt_soundspeed) then - do j = advxb, advxe - if ((j == advxe) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) + if (riemann_solver == 1) then + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe do k = 0, p do q = 0, n do l = 0, m rhs_vf(j)%sf(l, q, k) = & rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - (q_cons_vf%vf(j)%sf(l, q, k) - Kterm(l, q, k))* & - (flux_src_n(3)%vf(j)%sf(l, q, k) & - - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) + q_prim_vf%vf(contxe + idir)%sf(l, q, k)* & + (flux_src_n(3)%vf(j)%sf(l, q, k - 1) & + - flux_src_n(3)%vf(j)%sf(l, q, k)) end do end do end do - else if ((j == advxb) .and. (bubbles .neqv. .true.)) then - !$acc parallel loop collapse(3) gang vector default(present) - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - (q_cons_vf%vf(j)%sf(l, q, k) + Kterm(l, q, k))* & - (flux_src_n(3)%vf(j)%sf(l, q, k) & - - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) + end do + else + if (alt_soundspeed) then + do j = advxb, advxe + if ((j == advxe) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do k = 0, p + do q = 0, n + do l = 0, m + rhs_vf(j)%sf(l, q, k) = & + rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & + (q_cons_vf%vf(j)%sf(l, q, k) - Kterm(l, q, k))* & + (flux_src_n(3)%vf(j)%sf(l, q, k) & + - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) + end do + end do end do - end do + else if ((j == advxb) .and. (bubbles .neqv. .true.)) then + !$acc parallel loop collapse(3) gang vector default(present) + do k = 0, p + do q = 0, n + do l = 0, m + rhs_vf(j)%sf(l, q, k) = & + rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & + (q_cons_vf%vf(j)%sf(l, q, k) + Kterm(l, q, k))* & + (flux_src_n(3)%vf(j)%sf(l, q, k) & + - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) + end do + end do + end do + end if end do - end if - end do - else - !$acc parallel loop collapse(4) gang vector default(present) - do j = advxb, advxe - do k = 0, p - do q = 0, n - do l = 0, m - rhs_vf(j)%sf(l, q, k) = & - rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & - q_cons_vf%vf(j)%sf(l, q, k)* & - (flux_src_n(3)%vf(j)%sf(l, q, k) & - - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) + else + !$acc parallel loop collapse(4) gang vector default(present) + do j = advxb, advxe + do k = 0, p + do q = 0, n + do l = 0, m + rhs_vf(j)%sf(l, q, k) = & + rhs_vf(j)%sf(l, q, k) + 1d0/dz(k)* & + q_cons_vf%vf(j)%sf(l, q, k)* & + (flux_src_n(3)%vf(j)%sf(l, q, k) & + - flux_src_n(3)%vf(j)%sf(l, q, k - 1)) + end do + end do end do end do - end do - end do + end if + end if end if - end if - end if - end if ! id loop + end if ! id loop - end subroutine s_compute_advection_source_term + end subroutine s_compute_advection_source_term - subroutine s_compute_additional_physics_rhs(idir, q_prim_vf, rhs_vf, flux_src_n, & - dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf) + subroutine s_compute_additional_physics_rhs(idir, q_prim_vf, rhs_vf, flux_src_n, & + dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf) - integer, intent(in) :: idir - type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf - type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf - type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n - type(scalar_field), dimension(sys_size), intent(in) :: dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf + integer, intent(in) :: idir + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + type(scalar_field), dimension(sys_size), intent(inout) :: rhs_vf + type(scalar_field), dimension(sys_size), intent(in) :: flux_src_n + type(scalar_field), dimension(sys_size), intent(in) :: dq_prim_dx_vf, dq_prim_dy_vf, dq_prim_dz_vf - integer :: i, j, k, l, q + integer :: i, j, k, l, q - if (idir == 1) then ! x-direction + if (idir == 1) then ! x-direction - if (surface_tension) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - rhs_vf(c_idx)%sf(j, k, l) = & - rhs_vf(c_idx)%sf(j, k, l) + 1d0/dx(j)* & - q_prim_vf(c_idx)%sf(j, k, l)* & - (flux_src_n(advxb)%sf(j, k, l) - & - flux_src_n(advxb)%sf(j - 1, k, l)) + if (surface_tension) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + rhs_vf(c_idx)%sf(j, k, l) = & + rhs_vf(c_idx)%sf(j, k, l) + 1d0/dx(j)* & + q_prim_vf(c_idx)%sf(j, k, l)* & + (flux_src_n(advxb)%sf(j, k, l) - & + flux_src_n(advxb)%sf(j - 1, k, l)) + end do + end do end do - end do - end do - end if + end if - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - !$acc loop seq - do i = momxb, E_idx - rhs_vf(i)%sf(j, k, l) = & - rhs_vf(i)%sf(j, k, l) + 1d0/dx(j)* & - (flux_src_n(i)%sf(j - 1, k, l) & - - flux_src_n(i)%sf(j, k, l)) + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + !$acc loop seq + do i = momxb, E_idx + rhs_vf(i)%sf(j, k, l) = & + rhs_vf(i)%sf(j, k, l) + 1d0/dx(j)* & + (flux_src_n(i)%sf(j - 1, k, l) & + - flux_src_n(i)%sf(j, k, l)) + end do + end do end do end do - end do - end do - elseif (idir == 2) then ! y-direction + elseif (idir == 2) then ! y-direction - if (surface_tension) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - rhs_vf(c_idx)%sf(j, k, l) = & - rhs_vf(c_idx)%sf(j, k, l) + 1d0/dy(k)* & - q_prim_vf(c_idx)%sf(j, k, l)* & - (flux_src_n(advxb)%sf(j, k, l) - & - flux_src_n(advxb)%sf(j, k - 1, l)) + if (surface_tension) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + rhs_vf(c_idx)%sf(j, k, l) = & + rhs_vf(c_idx)%sf(j, k, l) + 1d0/dy(k)* & + q_prim_vf(c_idx)%sf(j, k, l)* & + (flux_src_n(advxb)%sf(j, k, l) - & + flux_src_n(advxb)%sf(j, k - 1, l)) + end do + end do end do - end do - end do - end if - - if (cyl_coord .and. ((bc_y%beg == -2) .or. (bc_y%beg == -14))) then - if (viscous) then - if (p > 0) then - call s_compute_viscous_stress_tensor(q_prim_vf, & - dq_prim_dx_vf(mom_idx%beg:mom_idx%end), & - dq_prim_dy_vf(mom_idx%beg:mom_idx%end), & - dq_prim_dz_vf(mom_idx%beg:mom_idx%end), & - tau_Re_vf, & - idwbuff(1), idwbuff(2), idwbuff(3)) - else - call s_compute_viscous_stress_tensor(q_prim_vf, & - dq_prim_dx_vf(mom_idx%beg:mom_idx%end), & - dq_prim_dy_vf(mom_idx%beg:mom_idx%end), & - dq_prim_dy_vf(mom_idx%beg:mom_idx%end), & - tau_Re_vf, & - idwbuff(1), idwbuff(2), idwbuff(3)) end if - !$acc parallel loop collapse(2) gang vector default(present) - do l = 0, p - do j = 0, m - !$acc loop seq - do i = momxb, E_idx - rhs_vf(i)%sf(j, 0, l) = & - rhs_vf(i)%sf(j, 0, l) + 1d0/(y_cc(1) - y_cc(-1))* & - (tau_Re_vf(i)%sf(j, -1, l) & - - tau_Re_vf(i)%sf(j, 1, l)) + if (cyl_coord .and. ((bc_y%beg == -2) .or. (bc_y%beg == -14))) then + if (viscous) then + if (p > 0) then + call s_compute_viscous_stress_tensor(q_prim_vf, & + dq_prim_dx_vf(mom_idx%beg:mom_idx%end), & + dq_prim_dy_vf(mom_idx%beg:mom_idx%end), & + dq_prim_dz_vf(mom_idx%beg:mom_idx%end), & + tau_Re_vf, & + idwbuff(1), idwbuff(2), idwbuff(3)) + else + call s_compute_viscous_stress_tensor(q_prim_vf, & + dq_prim_dx_vf(mom_idx%beg:mom_idx%end), & + dq_prim_dy_vf(mom_idx%beg:mom_idx%end), & + dq_prim_dy_vf(mom_idx%beg:mom_idx%end), & + tau_Re_vf, & + idwbuff(1), idwbuff(2), idwbuff(3)) + end if + + !$acc parallel loop collapse(2) gang vector default(present) + do l = 0, p + do j = 0, m + !$acc loop seq + do i = momxb, E_idx + rhs_vf(i)%sf(j, 0, l) = & + rhs_vf(i)%sf(j, 0, l) + 1d0/(y_cc(1) - y_cc(-1))* & + (tau_Re_vf(i)%sf(j, -1, l) & + - tau_Re_vf(i)%sf(j, 1, l)) + end do + end do end do - end do - end do - end if + end if - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 1, n - do j = 0, m - !$acc loop seq - do i = momxb, E_idx - rhs_vf(i)%sf(j, k, l) = & - rhs_vf(i)%sf(j, k, l) + 1d0/dy(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - - flux_src_n(i)%sf(j, k, l)) + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 1, n + do j = 0, m + !$acc loop seq + do i = momxb, E_idx + rhs_vf(i)%sf(j, k, l) = & + rhs_vf(i)%sf(j, k, l) + 1d0/dy(k)* & + (flux_src_n(i)%sf(j, k - 1, l) & + - flux_src_n(i)%sf(j, k, l)) + end do + end do end do end do - end do - end do - else - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - !$acc loop seq - do i = momxb, E_idx - rhs_vf(i)%sf(j, k, l) = & - rhs_vf(i)%sf(j, k, l) + 1d0/dy(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - - flux_src_n(i)%sf(j, k, l)) + else + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + !$acc loop seq + do i = momxb, E_idx + rhs_vf(i)%sf(j, k, l) = & + rhs_vf(i)%sf(j, k, l) + 1d0/dy(k)* & + (flux_src_n(i)%sf(j, k - 1, l) & + - flux_src_n(i)%sf(j, k, l)) + end do + end do end do end do - end do - end do - end if + end if - ! Applying the geometrical viscous Riemann source fluxes calculated as average - ! of values at cell boundaries - if (cyl_coord) then - if ((bc_y%beg == -2) .or. (bc_y%beg == -14)) then + ! Applying the geometrical viscous Riemann source fluxes calculated as average + ! of values at cell boundaries + if (cyl_coord) then + if ((bc_y%beg == -2) .or. (bc_y%beg == -14)) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 1, n - do j = 0, m - !$acc loop seq - do i = momxb, E_idx - rhs_vf(i)%sf(j, k, l) = & - rhs_vf(i)%sf(j, k, l) - 5d-1/y_cc(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - + flux_src_n(i)%sf(j, k, l)) + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 1, n + do j = 0, m + !$acc loop seq + do i = momxb, E_idx + rhs_vf(i)%sf(j, k, l) = & + rhs_vf(i)%sf(j, k, l) - 5d-1/y_cc(k)* & + (flux_src_n(i)%sf(j, k - 1, l) & + + flux_src_n(i)%sf(j, k, l)) + end do + end do + end do + end do + + if (viscous) then + !$acc parallel loop collapse(2) gang vector default(present) + do l = 0, p + do j = 0, m + !$acc loop seq + do i = momxb, E_idx + rhs_vf(i)%sf(j, 0, l) = & + rhs_vf(i)%sf(j, 0, l) - 1d0/y_cc(0)* & + tau_Re_vf(i)%sf(j, 0, l) + end do + end do + end do + end if + else + + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + !$acc loop seq + do i = momxb, E_idx + rhs_vf(i)%sf(j, k, l) = & + rhs_vf(i)%sf(j, k, l) - 5d-1/y_cc(k)* & + (flux_src_n(i)%sf(j, k - 1, l) & + + flux_src_n(i)%sf(j, k, l)) + end do + end do end do end do - end do - end do - if (viscous) then - !$acc parallel loop collapse(2) gang vector default(present) + end if + end if + + elseif (idir == 3) then ! z-direction + + if (surface_tension) then + !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p - do j = 0, m - !$acc loop seq - do i = momxb, E_idx - rhs_vf(i)%sf(j, 0, l) = & - rhs_vf(i)%sf(j, 0, l) - 1d0/y_cc(0)* & - tau_Re_vf(i)%sf(j, 0, l) + do k = 0, n + do j = 0, m + rhs_vf(c_idx)%sf(j, k, l) = & + rhs_vf(c_idx)%sf(j, k, l) + 1d0/dz(l)* & + q_prim_vf(c_idx)%sf(j, k, l)* & + (flux_src_n(advxb)%sf(j, k, l) - & + flux_src_n(advxb)%sf(j, k, l - 1)) end do end do end do end if - else !$acc parallel loop collapse(3) gang vector default(present) do l = 0, p @@ -1749,325 +1784,290 @@ contains !$acc loop seq do i = momxb, E_idx rhs_vf(i)%sf(j, k, l) = & - rhs_vf(i)%sf(j, k, l) - 5d-1/y_cc(k)* & - (flux_src_n(i)%sf(j, k - 1, l) & - + flux_src_n(i)%sf(j, k, l)) + rhs_vf(i)%sf(j, k, l) + 1d0/dz(l)* & + (flux_src_n(i)%sf(j, k, l - 1) & + - flux_src_n(i)%sf(j, k, l)) end do end do end do end do - end if - end if - - elseif (idir == 3) then ! z-direction - - if (surface_tension) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - rhs_vf(c_idx)%sf(j, k, l) = & - rhs_vf(c_idx)%sf(j, k, l) + 1d0/dz(l)* & - q_prim_vf(c_idx)%sf(j, k, l)* & - (flux_src_n(advxb)%sf(j, k, l) - & - flux_src_n(advxb)%sf(j, k, l - 1)) - end do - end do - end do - end if - - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - !$acc loop seq - do i = momxb, E_idx - rhs_vf(i)%sf(j, k, l) = & - rhs_vf(i)%sf(j, k, l) + 1d0/dz(l)* & - (flux_src_n(i)%sf(j, k, l - 1) & - - flux_src_n(i)%sf(j, k, l)) - end do - end do - end do - end do - - if (grid_geometry == 3) then - !$acc parallel loop collapse(3) gang vector default(present) - do l = 0, p - do k = 0, n - do j = 0, m - rhs_vf(momxb + 1)%sf(j, k, l) = & - rhs_vf(momxb + 1)%sf(j, k, l) + 5d-1* & - (flux_src_n(momxe)%sf(j, k, l - 1) & - + flux_src_n(momxe)%sf(j, k, l)) - - rhs_vf(momxe)%sf(j, k, l) = & - rhs_vf(momxe)%sf(j, k, l) - 5d-1* & - (flux_src_n(momxb + 1)%sf(j, k, l - 1) & - + flux_src_n(momxb + 1)%sf(j, k, l)) + if (grid_geometry == 3) then + !$acc parallel loop collapse(3) gang vector default(present) + do l = 0, p + do k = 0, n + do j = 0, m + rhs_vf(momxb + 1)%sf(j, k, l) = & + rhs_vf(momxb + 1)%sf(j, k, l) + 5d-1* & + (flux_src_n(momxe)%sf(j, k, l - 1) & + + flux_src_n(momxe)%sf(j, k, l)) + + rhs_vf(momxe)%sf(j, k, l) = & + rhs_vf(momxe)%sf(j, k, l) - 5d-1* & + (flux_src_n(momxb + 1)%sf(j, k, l - 1) & + + flux_src_n(momxb + 1)%sf(j, k, l)) + end do + end do end do - end do - end do - end if - end if + end if + end if - end subroutine s_compute_additional_physics_rhs + end subroutine s_compute_additional_physics_rhs - !> The purpose of this procedure is to infinitely relax + !> The purpose of this procedure is to infinitely relax !! the pressures from the internal-energy equations to a !! unique pressure, from which the corresponding volume !! fraction of each phase are recomputed. For conservation !! purpose, this pressure is finally corrected using the !! mixture-total-energy equation. !! @param q_cons_vf Cell-average conservative variables - subroutine s_pressure_relaxation_procedure(q_cons_vf) + subroutine s_pressure_relaxation_procedure(q_cons_vf) - type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf + type(scalar_field), dimension(sys_size), intent(inout) :: q_cons_vf - !> @name Relaxed pressure, initial partial pressures, function f(p) and its partial + !> @name Relaxed pressure, initial partial pressures, function f(p) and its partial !! derivative df(p), isentropic partial density, sum of volume fractions, !! mixture density, dynamic pressure, surface energy, specific heat ratio !! function, liquid stiffness function (two variations of the last two !! ones), shear and volume Reynolds numbers and the Weber numbers - !> @{ - real(kind(0d0)) :: pres_relax - real(kind(0d0)), dimension(num_fluids) :: pres_K_init - real(kind(0d0)) :: f_pres - real(kind(0d0)) :: df_pres - real(kind(0d0)), dimension(num_fluids) :: rho_K_s - real(kind(0d0)), dimension(num_fluids) :: alpha_rho - real(kind(0d0)), dimension(num_fluids) :: alpha - real(kind(0d0)) :: sum_alpha - real(kind(0d0)) :: rho - real(kind(0d0)) :: dyn_pres - real(kind(0d0)) :: gamma - real(kind(0d0)) :: pi_inf - real(kind(0d0)), dimension(2) :: Re - - integer :: i, j, k, l, q, iter !< Generic loop iterators - integer :: relax !< Relaxation procedure determination variable - - !$acc parallel loop collapse(3) gang vector private(pres_K_init, rho_K_s, alpha_rho, alpha, Re, pres_relax) - do l = 0, p - do k = 0, n - do j = 0, m - - ! Numerical correction of the volume fractions - if (mpp_lim) then - sum_alpha = 0d0 - - !$acc loop seq - do i = 1, num_fluids - if ((q_cons_vf(i + contxb - 1)%sf(j, k, l) < 0d0) .or. & - (q_cons_vf(i + advxb - 1)%sf(j, k, l) < 0d0)) then - q_cons_vf(i + contxb - 1)%sf(j, k, l) = 0d0 - q_cons_vf(i + advxb - 1)%sf(j, k, l) = 0d0 - q_cons_vf(i + intxb - 1)%sf(j, k, l) = 0d0 - end if - - if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > 1d0) & - q_cons_vf(i + advxb - 1)%sf(j, k, l) = 1d0 - sum_alpha = sum_alpha + q_cons_vf(i + advxb - 1)%sf(j, k, l) - end do - - !$acc loop seq - do i = 1, num_fluids - q_cons_vf(i + advxb - 1)%sf(j, k, l) = q_cons_vf(i + advxb - 1)%sf(j, k, l)/sum_alpha - end do - end if - - ! Pressures relaxation procedure =================================== - - ! Is the pressure relaxation procedure necessary? - relax = 1 + !> @{ + real(kind(0d0)) :: pres_relax + real(kind(0d0)), dimension(num_fluids) :: pres_K_init + real(kind(0d0)) :: f_pres + real(kind(0d0)) :: df_pres + real(kind(0d0)), dimension(num_fluids) :: rho_K_s + real(kind(0d0)), dimension(num_fluids) :: alpha_rho + real(kind(0d0)), dimension(num_fluids) :: alpha + real(kind(0d0)) :: sum_alpha + real(kind(0d0)) :: rho + real(kind(0d0)) :: dyn_pres + real(kind(0d0)) :: gamma + real(kind(0d0)) :: pi_inf + real(kind(0d0)), dimension(2) :: Re + + integer :: i, j, k, l, q, iter !< Generic loop iterators + integer :: relax !< Relaxation procedure determination variable + + !$acc parallel loop collapse(3) gang vector private(pres_K_init, rho_K_s, alpha_rho, alpha, Re, pres_relax) + do l = 0, p + do k = 0, n + do j = 0, m - !$acc loop seq - do i = 1, num_fluids - if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > (1d0 - sgm_eps)) relax = 0 - end do + ! Numerical correction of the volume fractions + if (mpp_lim) then + sum_alpha = 0d0 - if (relax == 1) then - ! Initial state - pres_relax = 0d0 + !$acc loop seq + do i = 1, num_fluids + if ((q_cons_vf(i + contxb - 1)%sf(j, k, l) < 0d0) .or. & + (q_cons_vf(i + advxb - 1)%sf(j, k, l) < 0d0)) then + q_cons_vf(i + contxb - 1)%sf(j, k, l) = 0d0 + q_cons_vf(i + advxb - 1)%sf(j, k, l) = 0d0 + q_cons_vf(i + intxb - 1)%sf(j, k, l) = 0d0 + end if - !$acc loop seq - do i = 1, num_fluids - if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) then - pres_K_init(i) = & - (q_cons_vf(i + intxb - 1)%sf(j, k, l)/ & - q_cons_vf(i + advxb - 1)%sf(j, k, l) & - - pi_infs(i))/gammas(i) + if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > 1d0) & + q_cons_vf(i + advxb - 1)%sf(j, k, l) = 1d0 + sum_alpha = sum_alpha + q_cons_vf(i + advxb - 1)%sf(j, k, l) + end do - if (pres_K_init(i) <= -(1d0 - 1d-8)*pres_inf(i) + 1d-8) & - pres_K_init(i) = -(1d0 - 1d-8)*pres_inf(i) + 1d-8 - else - pres_K_init(i) = 0d0 + !$acc loop seq + do i = 1, num_fluids + q_cons_vf(i + advxb - 1)%sf(j, k, l) = q_cons_vf(i + advxb - 1)%sf(j, k, l)/sum_alpha + end do end if - pres_relax = pres_relax + q_cons_vf(i + advxb - 1)%sf(j, k, l)*pres_K_init(i) - end do - ! Iterative process for relaxed pressure determination - f_pres = 1d-9 - df_pres = 1d9 + ! Pressures relaxation procedure =================================== - !$acc loop seq - do i = 1, num_fluids - rho_K_s(i) = 0d0 - end do + ! Is the pressure relaxation procedure necessary? + relax = 1 - !$acc loop seq - do iter = 0, 49 + !$acc loop seq + do i = 1, num_fluids + if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > (1d0 - sgm_eps)) relax = 0 + end do - if (DABS(f_pres) > 1d-10) then - pres_relax = pres_relax - f_pres/df_pres + if (relax == 1) then + ! Initial state + pres_relax = 0d0 - ! Physical pressure + !$acc loop seq do i = 1, num_fluids - if (pres_relax <= -(1d0 - 1d-8)*pres_inf(i) + 1d-8) & - pres_relax = -(1d0 - 1d-8)*pres_inf(i) + 1d0 + if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) then + pres_K_init(i) = & + (q_cons_vf(i + intxb - 1)%sf(j, k, l)/ & + q_cons_vf(i + advxb - 1)%sf(j, k, l) & + - pi_infs(i))/gammas(i) + + if (pres_K_init(i) <= -(1d0 - 1d-8)*pres_inf(i) + 1d-8) & + pres_K_init(i) = -(1d0 - 1d-8)*pres_inf(i) + 1d-8 + else + pres_K_init(i) = 0d0 + end if + pres_relax = pres_relax + q_cons_vf(i + advxb - 1)%sf(j, k, l)*pres_K_init(i) end do - ! Newton-Raphson method - f_pres = -1d0 - df_pres = 0d0 + ! Iterative process for relaxed pressure determination + f_pres = 1d-9 + df_pres = 1d9 !$acc loop seq do i = 1, num_fluids - if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) then - rho_K_s(i) = q_cons_vf(i + contxb - 1)%sf(j, k, l)/ & - max(q_cons_vf(i + advxb - 1)%sf(j, k, l), sgm_eps) & - *((pres_relax + pres_inf(i))/(pres_K_init(i) + & - pres_inf(i)))**(1d0/gamma_min(i)) + rho_K_s(i) = 0d0 + end do - f_pres = f_pres + q_cons_vf(i + contxb - 1)%sf(j, k, l) & - /rho_K_s(i) + !$acc loop seq + do iter = 0, 49 - df_pres = df_pres - q_cons_vf(i + contxb - 1)%sf(j, k, l) & - /(gamma_min(i)*rho_K_s(i)*(pres_relax + pres_inf(i))) - end if - end do - end if + if (DABS(f_pres) > 1d-10) then + pres_relax = pres_relax - f_pres/df_pres - end do + ! Physical pressure + do i = 1, num_fluids + if (pres_relax <= -(1d0 - 1d-8)*pres_inf(i) + 1d-8) & + pres_relax = -(1d0 - 1d-8)*pres_inf(i) + 1d0 + end do - ! Cell update of the volume fraction - !$acc loop seq - do i = 1, num_fluids - if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) & - q_cons_vf(i + advxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l) & - /rho_K_s(i) - end do - end if + ! Newton-Raphson method + f_pres = -1d0 + df_pres = 0d0 - ! ================================================================== + !$acc loop seq + do i = 1, num_fluids + if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) then + rho_K_s(i) = q_cons_vf(i + contxb - 1)%sf(j, k, l)/ & + max(q_cons_vf(i + advxb - 1)%sf(j, k, l), sgm_eps) & + *((pres_relax + pres_inf(i))/(pres_K_init(i) + & + pres_inf(i)))**(1d0/gamma_min(i)) - ! Mixture-total-energy correction ================================== + f_pres = f_pres + q_cons_vf(i + contxb - 1)%sf(j, k, l) & + /rho_K_s(i) - ! The mixture-total-energy correction of the mixture pressure P is not necessary here - ! because the primitive variables are directly recovered later on by the conservative - ! variables (see s_convert_conservative_to_primitive_variables called in s_compute_rhs). - ! However, the internal-energy equations should be reset with the corresponding mixture - ! pressure from the correction. This step is carried out below. + df_pres = df_pres - q_cons_vf(i + contxb - 1)%sf(j, k, l) & + /(gamma_min(i)*rho_K_s(i)*(pres_relax + pres_inf(i))) + end if + end do + end if - !$acc loop seq - do i = 1, num_fluids - alpha_rho(i) = q_cons_vf(i)%sf(j, k, l) - alpha(i) = q_cons_vf(E_idx + i)%sf(j, k, l) - end do + end do - if (bubbles) then - rho = 0d0 - gamma = 0d0 - pi_inf = 0d0 + ! Cell update of the volume fraction + !$acc loop seq + do i = 1, num_fluids + if (q_cons_vf(i + advxb - 1)%sf(j, k, l) > sgm_eps) & + q_cons_vf(i + advxb - 1)%sf(j, k, l) = q_cons_vf(i + contxb - 1)%sf(j, k, l) & + /rho_K_s(i) + end do + end if - if (mpp_lim .and. (model_eqns == 2) .and. (num_fluids > 2)) then - !$acc loop seq - do i = 1, num_fluids - rho = rho + alpha_rho(i) - gamma = gamma + alpha(i)*gammas(i) - pi_inf = pi_inf + alpha(i)*pi_infs(i) - end do - else if ((model_eqns == 2) .and. (num_fluids > 2)) then - !$acc loop seq - do i = 1, num_fluids - 1 - rho = rho + alpha_rho(i) - gamma = gamma + alpha(i)*gammas(i) - pi_inf = pi_inf + alpha(i)*pi_infs(i) - end do - else - rho = alpha_rho(1) - gamma = gammas(1) - pi_inf = pi_infs(1) - end if - else - rho = 0d0 - gamma = 0d0 - pi_inf = 0d0 + ! ================================================================== - sum_alpha = 0d0 + ! Mixture-total-energy correction ================================== + + ! The mixture-total-energy correction of the mixture pressure P is not necessary here + ! because the primitive variables are directly recovered later on by the conservative + ! variables (see s_convert_conservative_to_primitive_variables called in s_compute_rhs). + ! However, the internal-energy equations should be reset with the corresponding mixture + ! pressure from the correction. This step is carried out below. - if (mpp_lim) then !$acc loop seq do i = 1, num_fluids - alpha_rho(i) = max(0d0, alpha_rho(i)) - alpha(i) = min(max(0d0, alpha(i)), 1d0) - sum_alpha = sum_alpha + alpha(i) + alpha_rho(i) = q_cons_vf(i)%sf(j, k, l) + alpha(i) = q_cons_vf(E_idx + i)%sf(j, k, l) end do - alpha = alpha/max(sum_alpha, sgm_eps) - - end if + if (bubbles) then + rho = 0d0 + gamma = 0d0 + pi_inf = 0d0 + + if (mpp_lim .and. (model_eqns == 2) .and. (num_fluids > 2)) then + !$acc loop seq + do i = 1, num_fluids + rho = rho + alpha_rho(i) + gamma = gamma + alpha(i)*gammas(i) + pi_inf = pi_inf + alpha(i)*pi_infs(i) + end do + else if ((model_eqns == 2) .and. (num_fluids > 2)) then + !$acc loop seq + do i = 1, num_fluids - 1 + rho = rho + alpha_rho(i) + gamma = gamma + alpha(i)*gammas(i) + pi_inf = pi_inf + alpha(i)*pi_infs(i) + end do + else + rho = alpha_rho(1) + gamma = gammas(1) + pi_inf = pi_infs(1) + end if + else + rho = 0d0 + gamma = 0d0 + pi_inf = 0d0 + + sum_alpha = 0d0 + + if (mpp_lim) then + !$acc loop seq + do i = 1, num_fluids + alpha_rho(i) = max(0d0, alpha_rho(i)) + alpha(i) = min(max(0d0, alpha(i)), 1d0) + sum_alpha = sum_alpha + alpha(i) + end do - !$acc loop seq - do i = 1, num_fluids - rho = rho + alpha_rho(i) - gamma = gamma + alpha(i)*gammas(i) - pi_inf = pi_inf + alpha(i)*pi_infs(i) - end do + alpha = alpha/max(sum_alpha, sgm_eps) - if (viscous) then - !$acc loop seq - do i = 1, 2 - Re(i) = dflt_real + end if - if (Re_size(i) > 0) Re(i) = 0d0 !$acc loop seq - do q = 1, Re_size(i) - Re(i) = alpha(Re_idx(i, q))/Res(i, q) & - + Re(i) + do i = 1, num_fluids + rho = rho + alpha_rho(i) + gamma = gamma + alpha(i)*gammas(i) + pi_inf = pi_inf + alpha(i)*pi_infs(i) end do - Re(i) = 1d0/max(Re(i), sgm_eps) + if (viscous) then + !$acc loop seq + do i = 1, 2 + Re(i) = dflt_real - end do - end if - end if + if (Re_size(i) > 0) Re(i) = 0d0 + !$acc loop seq + do q = 1, Re_size(i) + Re(i) = alpha(Re_idx(i, q))/Res(i, q) & + + Re(i) + end do - dyn_pres = 0d0 + Re(i) = 1d0/max(Re(i), sgm_eps) - !$acc loop seq - do i = momxb, momxe - dyn_pres = dyn_pres + 5d-1*q_cons_vf(i)%sf(j, k, l)* & - q_cons_vf(i)%sf(j, k, l)/max(rho, sgm_eps) - end do + end do + end if + end if + + dyn_pres = 0d0 + + !$acc loop seq + do i = momxb, momxe + dyn_pres = dyn_pres + 5d-1*q_cons_vf(i)%sf(j, k, l)* & + q_cons_vf(i)%sf(j, k, l)/max(rho, sgm_eps) + end do - pres_relax = (q_cons_vf(E_idx)%sf(j, k, l) - dyn_pres - pi_inf)/gamma + pres_relax = (q_cons_vf(E_idx)%sf(j, k, l) - dyn_pres - pi_inf)/gamma - !$acc loop seq - do i = 1, num_fluids - q_cons_vf(i + intxb - 1)%sf(j, k, l) = & - q_cons_vf(i + advxb - 1)%sf(j, k, l)* & - (gammas(i)*pres_relax + pi_infs(i)) + !$acc loop seq + do i = 1, num_fluids + q_cons_vf(i + intxb - 1)%sf(j, k, l) = & + q_cons_vf(i + advxb - 1)%sf(j, k, l)* & + (gammas(i)*pres_relax + pi_infs(i)) + end do + ! ================================================================== + end do end do - ! ================================================================== end do - end do - end do - end subroutine s_pressure_relaxation_procedure + end subroutine s_pressure_relaxation_procedure - !> The purpose of this subroutine is to WENO-reconstruct the + !> The purpose of this subroutine is to WENO-reconstruct the !! left and the right cell-boundary values, including values !! at the Gaussian quadrature points, from the cell-averaged !! variables. @@ -2077,295 +2077,295 @@ contains !! @param vR_qp Right WENO-reconstructed, cell-boundary values including !! the values at the quadrature points, of the cell-average variables !! @param norm_dir Splitting coordinate direction - subroutine s_reconstruct_cell_boundary_values(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & - norm_dir) + subroutine s_reconstruct_cell_boundary_values(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & + norm_dir) - type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf - real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z - real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z - integer, intent(in) :: norm_dir + type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf + real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z + real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z + integer, intent(in) :: norm_dir - integer :: weno_dir !< Coordinate direction of the WENO reconstruction + integer :: weno_dir !< Coordinate direction of the WENO reconstruction - integer :: i, j, k, l - ! Reconstruction in s1-direction =================================== + integer :: i, j, k, l + ! Reconstruction in s1-direction =================================== - if (norm_dir == 1) then - is1 = idwbuff(1); is2 = idwbuff(2); is3 = idwbuff(3) - weno_dir = 1; is1%beg = is1%beg + weno_polyn - is1%end = is1%end - weno_polyn + if (norm_dir == 1) then + is1 = idwbuff(1); is2 = idwbuff(2); is3 = idwbuff(3) + weno_dir = 1; is1%beg = is1%beg + weno_polyn + is1%end = is1%end - weno_polyn - elseif (norm_dir == 2) then - is1 = idwbuff(2); is2 = idwbuff(1); is3 = idwbuff(3) - weno_dir = 2; is1%beg = is1%beg + weno_polyn - is1%end = is1%end - weno_polyn + elseif (norm_dir == 2) then + is1 = idwbuff(2); is2 = idwbuff(1); is3 = idwbuff(3) + weno_dir = 2; is1%beg = is1%beg + weno_polyn + is1%end = is1%end - weno_polyn - else - is1 = idwbuff(3); is2 = idwbuff(2); is3 = idwbuff(1) - weno_dir = 3; is1%beg = is1%beg + weno_polyn - is1%end = is1%end - weno_polyn + else + is1 = idwbuff(3); is2 = idwbuff(2); is3 = idwbuff(1) + weno_dir = 3; is1%beg = is1%beg + weno_polyn + is1%end = is1%end - weno_polyn - end if + end if - if (n > 0) then - if (p > 0) then - - call s_weno(v_vf(iv%beg:iv%end), & - vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, iv%beg:iv%end), vL_z(:, :, :, iv%beg:iv%end), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, iv%beg:iv%end), vR_z(:, :, :, iv%beg:iv%end), & - norm_dir, weno_dir, & - is1, is2, is3) - else - call s_weno(v_vf(iv%beg:iv%end), & - vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, iv%beg:iv%end), vL_z(:, :, :, :), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, iv%beg:iv%end), vR_z(:, :, :, :), & - norm_dir, weno_dir, & - is1, is2, is3) - end if - else + if (n > 0) then + if (p > 0) then - call s_weno(v_vf(iv%beg:iv%end), & - vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, :), vL_z(:, :, :, :), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, :), vR_z(:, :, :, :), & - norm_dir, weno_dir, & - is1, is2, is3) - end if + call s_weno(v_vf(iv%beg:iv%end), & + vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, iv%beg:iv%end), vL_z(:, :, :, iv%beg:iv%end), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, iv%beg:iv%end), vR_z(:, :, :, iv%beg:iv%end), & + norm_dir, weno_dir, & + is1, is2, is3) + else + call s_weno(v_vf(iv%beg:iv%end), & + vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, iv%beg:iv%end), vL_z(:, :, :, :), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, iv%beg:iv%end), vR_z(:, :, :, :), & + norm_dir, weno_dir, & + is1, is2, is3) + end if + else - ! ================================================================== - end subroutine s_reconstruct_cell_boundary_values + call s_weno(v_vf(iv%beg:iv%end), & + vL_x(:, :, :, iv%beg:iv%end), vL_y(:, :, :, :), vL_z(:, :, :, :), vR_x(:, :, :, iv%beg:iv%end), vR_y(:, :, :, :), vR_z(:, :, :, :), & + norm_dir, weno_dir, & + is1, is2, is3) + end if - subroutine s_reconstruct_cell_boundary_values_first_order(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & - norm_dir) + ! ================================================================== + end subroutine s_reconstruct_cell_boundary_values - type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf - real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z - real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z - integer, intent(in) :: norm_dir + subroutine s_reconstruct_cell_boundary_values_first_order(v_vf, vL_x, vL_y, vL_z, vR_x, vR_y, vR_z, & + norm_dir) - integer :: recon_dir !< Coordinate direction of the WENO reconstruction + type(scalar_field), dimension(iv%beg:iv%end), intent(in) :: v_vf + real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vL_x, vL_y, vL_z + real(kind(0d0)), dimension(startx:, starty:, startz:, 1:), intent(inout) :: vR_x, vR_y, vR_z + integer, intent(in) :: norm_dir - integer :: i, j, k, l - ! Reconstruction in s1-direction =================================== + integer :: recon_dir !< Coordinate direction of the WENO reconstruction - if (norm_dir == 1) then - is1 = idwbuff(1); is2 = idwbuff(2); is3 = idwbuff(3) - recon_dir = 1; is1%beg = is1%beg + weno_polyn - is1%end = is1%end - weno_polyn + integer :: i, j, k, l + ! Reconstruction in s1-direction =================================== - elseif (norm_dir == 2) then - is1 = idwbuff(2); is2 = idwbuff(1); is3 = idwbuff(3) - recon_dir = 2; is1%beg = is1%beg + weno_polyn - is1%end = is1%end - weno_polyn + if (norm_dir == 1) then + is1 = idwbuff(1); is2 = idwbuff(2); is3 = idwbuff(3) + recon_dir = 1; is1%beg = is1%beg + weno_polyn + is1%end = is1%end - weno_polyn - else - is1 = idwbuff(3); is2 = idwbuff(2); is3 = idwbuff(1) - recon_dir = 3; is1%beg = is1%beg + weno_polyn - is1%end = is1%end - weno_polyn + elseif (norm_dir == 2) then + is1 = idwbuff(2); is2 = idwbuff(1); is3 = idwbuff(3) + recon_dir = 2; is1%beg = is1%beg + weno_polyn + is1%end = is1%end - weno_polyn - end if + else + is1 = idwbuff(3); is2 = idwbuff(2); is3 = idwbuff(1) + recon_dir = 3; is1%beg = is1%beg + weno_polyn + is1%end = is1%end - weno_polyn + + end if #ifndef _CRAYFTN !$acc update device(is1, is2, is3, iv) #endif - if (recon_dir == 1) then - !$acc parallel loop collapse(4) default(present) - do i = iv%beg, iv%end - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end - vL_x(j, k, l, i) = v_vf(i)%sf(j, k, l) - vR_x(j, k, l, i) = v_vf(i)%sf(j, k, l) + if (recon_dir == 1) then + !$acc parallel loop collapse(4) default(present) + do i = iv%beg, iv%end + do l = is3%beg, is3%end + do k = is2%beg, is2%end + do j = is1%beg, is1%end + vL_x(j, k, l, i) = v_vf(i)%sf(j, k, l) + vR_x(j, k, l, i) = v_vf(i)%sf(j, k, l) + end do + end do end do end do - end do - end do - !$acc end parallel loop - else if (recon_dir == 2) then - !$acc parallel loop collapse(4) default(present) - do i = iv%beg, iv%end - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end - vL_y(j, k, l, i) = v_vf(i)%sf(k, j, l) - vR_y(j, k, l, i) = v_vf(i)%sf(k, j, l) + !$acc end parallel loop + else if (recon_dir == 2) then + !$acc parallel loop collapse(4) default(present) + do i = iv%beg, iv%end + do l = is3%beg, is3%end + do k = is2%beg, is2%end + do j = is1%beg, is1%end + vL_y(j, k, l, i) = v_vf(i)%sf(k, j, l) + vR_y(j, k, l, i) = v_vf(i)%sf(k, j, l) + end do + end do end do end do - end do - end do - !$acc end parallel loop - else if (recon_dir == 3) then - !$acc parallel loop collapse(4) default(present) - do i = iv%beg, iv%end - do l = is3%beg, is3%end - do k = is2%beg, is2%end - do j = is1%beg, is1%end - vL_z(j, k, l, i) = v_vf(i)%sf(l, k, j) - vR_z(j, k, l, i) = v_vf(i)%sf(l, k, j) + !$acc end parallel loop + else if (recon_dir == 3) then + !$acc parallel loop collapse(4) default(present) + do i = iv%beg, iv%end + do l = is3%beg, is3%end + do k = is2%beg, is2%end + do j = is1%beg, is1%end + vL_z(j, k, l, i) = v_vf(i)%sf(l, k, j) + vR_z(j, k, l, i) = v_vf(i)%sf(l, k, j) + end do + end do end do end do - end do - end do - !$acc end parallel loop - end if - - end subroutine s_reconstruct_cell_boundary_values_first_order - - !> Module deallocation and/or disassociation procedures - subroutine s_finalize_rhs_module - - integer :: i, j, k, l !< Generic loop iterators - - do j = cont_idx%beg, cont_idx%end - !$acc exit data detach(q_prim_qp%vf(j)%sf) - nullify (q_prim_qp%vf(j)%sf) - end do - - do j = adv_idx%beg, adv_idx%end - !$acc exit data detach(q_prim_qp%vf(j)%sf) - nullify (q_prim_qp%vf(j)%sf) - end do + !$acc end parallel loop + end if - do j = mom_idx%beg, E_idx - @:DEALLOCATE(q_cons_qp%vf(j)%sf) - @:DEALLOCATE(q_prim_qp%vf(j)%sf) - end do + end subroutine s_reconstruct_cell_boundary_values_first_order - @:DEALLOCATE(q_cons_qp%vf, q_prim_qp%vf) - @:DEALLOCATE_GLOBAL(qL_rsx_vf, qR_rsx_vf) + !> Module deallocation and/or disassociation procedures + subroutine s_finalize_rhs_module - if (n > 0) then - @:DEALLOCATE_GLOBAL(qL_rsy_vf, qR_rsy_vf) - end if + integer :: i, j, k, l !< Generic loop iterators - if (p > 0) then - @:DEALLOCATE_GLOBAL(qL_rsz_vf, qR_rsz_vf) - end if + do j = cont_idx%beg, cont_idx%end + !$acc exit data detach(q_prim_qp%vf(j)%sf) + nullify (q_prim_qp%vf(j)%sf) + end do - if (viscous .and. weno_Re_flux) then - @:DEALLOCATE_GLOBAL(dqL_rsx_vf, dqR_rsx_vf) + do j = adv_idx%beg, adv_idx%end + !$acc exit data detach(q_prim_qp%vf(j)%sf) + nullify (q_prim_qp%vf(j)%sf) + end do - if (n > 0) then - @:DEALLOCATE_GLOBAL(dqL_rsy_vf, dqR_rsy_vf) - end if + do j = mom_idx%beg, E_idx + @:DEALLOCATE(q_cons_qp%vf(j)%sf) + @:DEALLOCATE(q_prim_qp%vf(j)%sf) + end do - if (p > 0) then - @:DEALLOCATE_GLOBAL(dqL_rsz_vf, dqR_rsz_vf) - end if - end if + @:DEALLOCATE(q_cons_qp%vf, q_prim_qp%vf) + @:DEALLOCATE_GLOBAL(qL_rsx_vf, qR_rsx_vf) - if (mpp_lim .and. bubbles) then - !$acc exit data delete(alf_sum%sf) - deallocate (alf_sum%sf) - end if + if (n > 0) then + @:DEALLOCATE_GLOBAL(qL_rsy_vf, qR_rsy_vf) + end if - if (viscous) then - do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf) - end do + if (p > 0) then + @:DEALLOCATE_GLOBAL(qL_rsz_vf, qR_rsz_vf) + end if - if (n > 0) then + if (viscous .and. weno_Re_flux) then + @:DEALLOCATE_GLOBAL(dqL_rsx_vf, dqR_rsx_vf) - do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dq_prim_dy_qp(1)%vf(l)%sf) - end do + if (n > 0) then + @:DEALLOCATE_GLOBAL(dqL_rsy_vf, dqR_rsy_vf) + end if - if (p > 0) then - do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dq_prim_dz_qp(1)%vf(l)%sf) - end do + if (p > 0) then + @:DEALLOCATE_GLOBAL(dqL_rsz_vf, dqR_rsz_vf) + end if end if - end if - - @:DEALLOCATE(dq_prim_dx_qp(1)%vf) - @:DEALLOCATE(dq_prim_dy_qp(1)%vf) - @:DEALLOCATE(dq_prim_dz_qp(1)%vf) - end if + if (mpp_lim .and. bubbles) then + !$acc exit data delete(alf_sum%sf) + deallocate (alf_sum%sf) + end if - if (viscous) then - do i = num_dims, 1, -1 if (viscous) then - do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf) - @:DEALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf) + @:DEALLOCATE(dq_prim_dx_qp(1)%vf(l)%sf) end do if (n > 0) then - do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf) - @:DEALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf) - end do - end if - if (p > 0) then do l = mom_idx%beg, mom_idx%end - @:DEALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf) - @:DEALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf) + @:DEALLOCATE(dq_prim_dy_qp(1)%vf(l)%sf) end do + + if (p > 0) then + do l = mom_idx%beg, mom_idx%end + @:DEALLOCATE(dq_prim_dz_qp(1)%vf(l)%sf) + end do + end if + end if + @:DEALLOCATE(dq_prim_dx_qp(1)%vf) + @:DEALLOCATE(dq_prim_dy_qp(1)%vf) + @:DEALLOCATE(dq_prim_dz_qp(1)%vf) end if - @:DEALLOCATE(dqL_prim_dx_n(i)%vf) - @:DEALLOCATE(dqL_prim_dy_n(i)%vf) - @:DEALLOCATE(dqL_prim_dz_n(i)%vf) - @:DEALLOCATE(dqR_prim_dx_n(i)%vf) - @:DEALLOCATE(dqR_prim_dy_n(i)%vf) - @:DEALLOCATE(dqR_prim_dz_n(i)%vf) - end do - end if + if (viscous) then + do i = num_dims, 1, -1 + if (viscous) then - @:DEALLOCATE_GLOBAL(dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n) - @:DEALLOCATE_GLOBAL(dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n) + do l = mom_idx%beg, mom_idx%end + @:DEALLOCATE(dqL_prim_dx_n(i)%vf(l)%sf) + @:DEALLOCATE(dqR_prim_dx_n(i)%vf(l)%sf) + end do - do i = num_dims, 1, -1 - if (i /= 1) then - do l = 1, sys_size - nullify (flux_n(i)%vf(l)%sf) - nullify (flux_src_n(i)%vf(l)%sf) - @:DEALLOCATE(flux_gsrc_n(i)%vf(l)%sf) - end do - else - do l = 1, sys_size - @:DEALLOCATE(flux_n(i)%vf(l)%sf) - @:DEALLOCATE(flux_gsrc_n(i)%vf(l)%sf) - end do + if (n > 0) then + do l = mom_idx%beg, mom_idx%end + @:DEALLOCATE(dqL_prim_dy_n(i)%vf(l)%sf) + @:DEALLOCATE(dqR_prim_dy_n(i)%vf(l)%sf) + end do + end if - if (viscous) then - do l = mom_idx%beg, E_idx - @:DEALLOCATE(flux_src_n(i)%vf(l)%sf) - end do - end if + if (p > 0) then + do l = mom_idx%beg, mom_idx%end + @:DEALLOCATE(dqL_prim_dz_n(i)%vf(l)%sf) + @:DEALLOCATE(dqR_prim_dz_n(i)%vf(l)%sf) + end do + end if - if (riemann_solver == 1) then - do l = adv_idx%beg + 1, adv_idx%end - @:DEALLOCATE(flux_src_n(i)%vf(l)%sf) - end do - else - do l = adv_idx%beg + 1, adv_idx%end - nullify (flux_src_n(i)%vf(l)%sf) + end if + + @:DEALLOCATE(dqL_prim_dx_n(i)%vf) + @:DEALLOCATE(dqL_prim_dy_n(i)%vf) + @:DEALLOCATE(dqL_prim_dz_n(i)%vf) + @:DEALLOCATE(dqR_prim_dx_n(i)%vf) + @:DEALLOCATE(dqR_prim_dy_n(i)%vf) + @:DEALLOCATE(dqR_prim_dz_n(i)%vf) end do end if - @:DEALLOCATE(flux_src_n(i)%vf(adv_idx%beg)%sf) - end if + @:DEALLOCATE_GLOBAL(dqL_prim_dx_n, dqL_prim_dy_n, dqL_prim_dz_n) + @:DEALLOCATE_GLOBAL(dqR_prim_dx_n, dqR_prim_dy_n, dqR_prim_dz_n) - @:DEALLOCATE(flux_n(i)%vf, flux_src_n(i)%vf, flux_gsrc_n(i)%vf) - end do + do i = num_dims, 1, -1 + if (i /= 1) then + do l = 1, sys_size + nullify (flux_n(i)%vf(l)%sf) + nullify (flux_src_n(i)%vf(l)%sf) + @:DEALLOCATE(flux_gsrc_n(i)%vf(l)%sf) + end do + else + do l = 1, sys_size + @:DEALLOCATE(flux_n(i)%vf(l)%sf) + @:DEALLOCATE(flux_gsrc_n(i)%vf(l)%sf) + end do + + if (viscous) then + do l = mom_idx%beg, E_idx + @:DEALLOCATE(flux_src_n(i)%vf(l)%sf) + end do + end if - @:DEALLOCATE_GLOBAL(flux_n, flux_src_n, flux_gsrc_n) + if (riemann_solver == 1) then + do l = adv_idx%beg + 1, adv_idx%end + @:DEALLOCATE(flux_src_n(i)%vf(l)%sf) + end do + else + do l = adv_idx%beg + 1, adv_idx%end + nullify (flux_src_n(i)%vf(l)%sf) + end do + end if - if (viscous .and. cyl_coord) then - do i = 1, num_dims - @:DEALLOCATE(tau_re_vf(cont_idx%end + i)%sf) - end do - @:DEALLOCATE(tau_re_vf(e_idx)%sf) - @:DEALLOCATE_GLOBAL(tau_re_vf) - end if + @:DEALLOCATE(flux_src_n(i)%vf(adv_idx%beg)%sf) + end if + + @:DEALLOCATE(flux_n(i)%vf, flux_src_n(i)%vf, flux_gsrc_n(i)%vf) + end do + + @:DEALLOCATE_GLOBAL(flux_n, flux_src_n, flux_gsrc_n) + + if (viscous .and. cyl_coord) then + do i = 1, num_dims + @:DEALLOCATE(tau_re_vf(cont_idx%end + i)%sf) + end do + @:DEALLOCATE(tau_re_vf(e_idx)%sf) + @:DEALLOCATE_GLOBAL(tau_re_vf) + end if - s_riemann_solver => null() - s_convert_to_mixture_variables => null() + s_riemann_solver => null() + s_convert_to_mixture_variables => null() - end subroutine s_finalize_rhs_module + end subroutine s_finalize_rhs_module -end module m_rhs + end module m_rhs