Skip to content

Commit 16e51c9

Browse files
wilfonbawilfonba
authored
Make surface tension compatible with the 5-eqn model and immersed boundaries (#821)
1 parent be77ab2 commit 16e51c9

File tree

5 files changed

+42
-11
lines changed

5 files changed

+42
-11
lines changed

‎.github/workflows/frontier/submit.sh

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@ sbatch <<EOT
2121
#SBATCH -A CFD154 # charge account
2222
#SBATCH -N 1 # Number of nodes required
2323
#SBATCH -n 8 # Number of cores required
24-
#SBATCH -t 01:00:00 # Duration of the job (Ex: 15 mins)
24+
#SBATCH -t 01:30:00 # Duration of the job (Ex: 15 mins)
2525
#SBATCH -o$job_slug.out # Combined output and error messages file
2626
#SBATCH -q debug # Use debug QOS - only one job per user allowed in queue!
2727
#SBATCH -W # Do not exit until the submitted job terminates.

‎src/common/m_checker_common.fpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -316,8 +316,8 @@ contains
316316
@:PROHIBIT(.not. f_is_default(sigma) .and. .not. surface_tension, &
317317
"sigma is set but surface_tension is not enabled")
318318

319-
@:PROHIBIT(surface_tension .and. model_eqns /= 3, &
320-
"The surface tension model requires model_eqns=3")
319+
@:PROHIBIT(surface_tension .and. (model_eqns /= 3 .and. model_eqns /=2), &
320+
"The surface tension model requires model_eqns=3 or model_eqns=2")
321321

322322
@:PROHIBIT(surface_tension .and. num_fluids /= 2, &
323323
"The surface tension model requires num_fluids=2")

‎src/simulation/m_global_parameters.fpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -296,7 +296,7 @@ module m_global_parameters
296296
!! conditions data to march the solution in the physical computational domain
297297
!! to the next time-step.
298298

299-
!$acc declare create(sys_size, buff_size, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, b_size, tensor_size, xi_idx, species_idx, B_idx)
299+
!$acc declare create(sys_size, buff_size, E_idx, gamma_idx, pi_inf_idx, alf_idx, n_idx, stress_idx, b_size, tensor_size, xi_idx, species_idx, B_idx, c_idx)
300300

301301
integer :: shear_num !! Number of shear stress components
302302
integer, dimension(3) :: shear_indices !<
@@ -1212,7 +1212,7 @@ contains
12121212
chemxb = species_idx%beg
12131213
chemxe = species_idx%end
12141214
1215-
!$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe)
1215+
!$acc update device(momxb, momxe, advxb, advxe, contxb, contxe, bubxb, bubxe, intxb, intxe, sys_size, buff_size, E_idx, alf_idx, n_idx, adv_n, adap_dt, pi_fac, strxb, strxe, chemxb, chemxe, c_idx)
12161216
!$acc update device(b_size, xibeg, xiend, tensor_size)
12171217
12181218
!$acc update device(species_idx)

‎src/simulation/m_ibm.fpp

Lines changed: 27 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -135,6 +135,7 @@ contains
135135

136136
real(wp) :: pres_IP, coeff
137137
real(wp), dimension(3) :: vel_IP, vel_norm_IP
138+
real(wp) :: c_IP
138139
real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP
139140
real(wp), dimension(nb) :: r_IP, v_IP, pb_IP, mv_IP
140141
real(wp), dimension(nb*nmom) :: nmom_IP
@@ -170,19 +171,19 @@ contains
170171
!Interpolate primitive variables at image point associated w/ GP
171172
if (bubbles_euler .and. .not. qbmm) then
172173
call s_interpolate_image_point(q_prim_vf, gp, &
173-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
174+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
174175
r_IP, v_IP, pb_IP, mv_IP)
175176
else if (qbmm .and. polytropic) then
176177
call s_interpolate_image_point(q_prim_vf, gp, &
177-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
178+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
178179
r_IP, v_IP, pb_IP, mv_IP, nmom_IP)
179180
else if (qbmm .and. .not. polytropic) then
180181
call s_interpolate_image_point(q_prim_vf, gp, &
181-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
182+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
182183
r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP)
183184
else
184185
call s_interpolate_image_point(q_prim_vf, gp, &
185-
alpha_rho_IP, alpha_IP, pres_IP, vel_IP)
186+
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP)
186187
end if
187188

188189
dyn_pres = 0._wp
@@ -194,6 +195,10 @@ contains
194195
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
195196
end do
196197

198+
if (surface_tension) then
199+
q_prim_vf(c_idx)%sf(j, k, l) = c_IP
200+
end if
201+
197202
if (model_eqns /= 4) then
198203
! If in simulation, use acc mixture subroutines
199204
if (elasticity) then
@@ -234,6 +239,11 @@ contains
234239
q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
235240
end do
236241

242+
! Set color function
243+
if (surface_tension) then
244+
q_cons_vf(c_idx)%sf(j, k, l) = c_IP
245+
end if
246+
237247
! Set Energy
238248
if (bubbles_euler) then
239249
q_cons_vf(E_idx)%sf(j, k, l) = (1 - alpha_IP(1))*(gamma*pres_IP + pi_inf + dyn_pres)
@@ -309,6 +319,10 @@ contains
309319
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
310320
end do
311321

322+
if (surface_tension) then
323+
q_prim_vf(c_idx)%sf(j, k, l) = c_IP
324+
end if
325+
312326
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
313327
alpha_rho_IP, Re_K, j, k, l)
314328

@@ -725,16 +739,18 @@ contains
725739

726740
!> Function that uses the interpolation coefficients and the current state
727741
!! at the cell centers in order to estimate the state at the image point
728-
subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP)
742+
subroutine s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, mv_IP, nmom_IP, pb, mv, presb_IP, massv_IP)
729743
!$acc routine seq
730744
type(scalar_field), &
731745
dimension(sys_size), &
732746
intent(IN) :: q_prim_vf !< Primitive Variables
747+
733748
real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(INOUT) :: pb, mv
734749

735750
type(ghost_point), intent(IN) :: gp
736751
real(wp), intent(INOUT) :: pres_IP
737752
real(wp), dimension(3), intent(INOUT) :: vel_IP
753+
real(wp), intent(INOUT) :: c_IP
738754
real(wp), dimension(num_fluids), intent(INOUT) :: alpha_IP, alpha_rho_IP
739755
real(wp), optional, dimension(:), intent(INOUT) :: r_IP, v_IP, pb_IP, mv_IP
740756
real(wp), optional, dimension(:), intent(INOUT) :: nmom_IP
@@ -758,6 +774,8 @@ contains
758774
pres_IP = 0._wp
759775
vel_IP = 0._wp
760776

777+
if (surface_tension) c_IP = 0._wp
778+
761779
if (bubbles_euler) then
762780
r_IP = 0._wp
763781
v_IP = 0._wp
@@ -801,6 +819,10 @@ contains
801819
q_prim_vf(advxb + l - 1)%sf(i, j, k)
802820
end do
803821

822+
if (surface_tension) then
823+
c_IP = c_IP + coeff*q_prim_vf(c_idx)%sf(i, j, k)
824+
end if
825+
804826
if (bubbles_euler .and. .not. qbmm) then
805827
!$acc loop seq
806828
do l = 1, nb

‎src/simulation/m_riemann_solvers.fpp

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1650,7 +1650,7 @@ contains
16501650
end do
16511651
end if
16521652
1653-
! SURFACE TENSION FLUX. need to check
1653+
! COLOR FUNCTION FLUX
16541654
if (surface_tension) then
16551655
flux_rs${XYZ}$_vf(j, k, l, c_idx) = &
16561656
(xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, c_idx) + &
@@ -2831,6 +2831,15 @@ contains
28312831
s_P*(xi_R - 1._wp))
28322832
end do
28332833
2834+
! COLOR FUNCTION FLUX
2835+
if (surface_tension) then
2836+
flux_rs${XYZ}$_vf(j, k, l, c_idx) = &
2837+
xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, c_idx) &
2838+
*(vel_L(idx1) + s_M*(xi_L - 1._wp)) &
2839+
+ xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, c_idx) &
2840+
*(vel_R(idx1) + s_P*(xi_R - 1._wp))
2841+
end if
2842+
28342843
! REFERENCE MAP FLUX.
28352844
if (hyperelasticity) then
28362845
!$acc loop seq

0 commit comments

Comments
 (0)