Skip to content

Make surface tension compatible with the 5-eqn model and immersed boundaries #821

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Apr 26, 2025
Merged

Make surface tension compatible with the 5-eqn model and immersed boundaries #821

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
70bb126
5eqn and IBM ST
wilfonba Apr 20, 2025
84c938a
format
wilfonba Apr 20, 2025
9309268
longer frontier runner
wilfonba Apr 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .github/workflows/frontier/submit.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@ sbatch <<EOT
#SBATCH -A CFD154 # charge account
#SBATCH -N 1 # Number of nodes required
#SBATCH -n 8 # Number of cores required
#SBATCH -t 01:00:00 # Duration of the job (Ex: 15 mins)
#SBATCH -t 01:30:00 # Duration of the job (Ex: 15 mins)
#SBATCH -o$job_slug.out # Combined output and error messages file
#SBATCH -q debug # Use debug QOS - only one job per user allowed in queue!
#SBATCH -W # Do not exit until the submitted job terminates.
Expand Down
4 changes: 2 additions & 2 deletions src/common/m_checker_common.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -316,8 +316,8 @@ contains
@:PROHIBIT(.not. f_is_default(sigma) .and. .not. surface_tension, &
"sigma is set but surface_tension is not enabled")

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

@:PROHIBIT(surface_tension .and. num_fluids /= 2, &
"The surface tension model requires num_fluids=2")
Expand Down
4 changes: 2 additions & 2 deletions src/simulation/m_global_parameters.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,7 +296,7 @@ module m_global_parameters
!! conditions data to march the solution in the physical computational domain
!! to the next time-step.

!$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)
!$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)

integer :: shear_num !! Number of shear stress components
integer, dimension(3) :: shear_indices !<
Expand Down Expand Up @@ -1212,7 +1212,7 @@ contains
chemxb = species_idx%beg
chemxe = species_idx%end

!$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)
!$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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can we just do eqn_idx%<whatever> where <whatever> -> {color,energy,alpha,...} please open a PR for it.

Copy link
Collaborator

@Malmahrouqi3 Malmahrouqi3 Jun 5, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just the variables at line 230?

eqn_idx%E_idx                                   !< Index of energy equation
eqn_idx%n_idx                                   !< Index of number density
eqn_idx%alf_idx                                 !< Index of void fraction
eqn_idx%gamma_idx                         !< Index of specific heat ratio func. eqn.
eqn_idx%pi_inf_idx                             !< Index of liquid stiffness func. eqn.
eqn_idx%b_size                                  !< Number of elements in the symmetric b tensor, plus one
eqn_idx%tensor_size                          !< Number of elements in the full tensor plus one
eqn_idx%c_idx                                    !< Index of color function
eqn_idx%damage_idx                         !< Index of damage state variable (D) for continuum

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

i meant all of the equations @mohdsaid497566 . you will see there are many other _idx variables (i think?)

!$acc update device(b_size, xibeg, xiend, tensor_size)

!$acc update device(species_idx)
Expand Down
32 changes: 27 additions & 5 deletions src/simulation/m_ibm.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,7 @@

real(wp) :: pres_IP, coeff
real(wp), dimension(3) :: vel_IP, vel_norm_IP
real(wp) :: c_IP
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't each IP be a derived type with components alpha_rho,alpha,pressure,vel,c,r,v,pb,mv,nmom,...?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

added "ip" where type(point_data) :: ip

type :: point_data
real(wp), dimension(:), allocatable :: alpha_rho  !< Partial densities
real(wp), dimension(:), allocatable :: alpha      !< Volume fractions
real(wp) :: pressure                              !< Pressure
real(wp), dimension(3) :: vel                     !< Velocity
real(wp) :: c                                     !< Color function (for surface tension)
real(wp), dimension(:), allocatable :: r          !< Bubble radii
real(wp), dimension(:), allocatable :: v          !< Bubble radial velocities
real(wp), dimension(:), allocatable :: pb         !< Bubble pressures
real(wp), dimension(:), allocatable :: mv         !< Mass of vapor
real(wp), dimension(:), allocatable :: nmom       !< Moments for QBMM
real(wp), dimension(:), allocatable :: presb      !< Node pressures for bubbles
real(wp), dimension(:), allocatable :: massv      !< Node masses for bubbles
end type point_data

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nice, assuming it works of course

real(wp), dimension(num_fluids) :: alpha_rho_IP, alpha_IP
real(wp), dimension(nb) :: r_IP, v_IP, pb_IP, mv_IP
real(wp), dimension(nb*nmom) :: nmom_IP
Expand Down Expand Up @@ -170,19 +171,19 @@
!Interpolate primitive variables at image point associated w/ GP
if (bubbles_euler .and. .not. qbmm) then
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
r_IP, v_IP, pb_IP, mv_IP)
else if (qbmm .and. polytropic) then
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, &
r_IP, v_IP, pb_IP, mv_IP, nmom_IP)
else if (qbmm .and. .not. polytropic) then
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, &
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)
else
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP)
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP)
end if

dyn_pres = 0._wp
Expand All @@ -194,6 +195,10 @@
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
end do

if (surface_tension) then
q_prim_vf(c_idx)%sf(j, k, l) = c_IP

Check warning on line 199 in src/simulation/m_ibm.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_ibm.fpp#L199 

Added line #L199 was not covered by tests
end if

if (model_eqns /= 4) then
! If in simulation, use acc mixture subroutines
if (elasticity) then
Expand Down Expand Up @@ -234,6 +239,11 @@
q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
end do

! Set color function
if (surface_tension) then
q_cons_vf(c_idx)%sf(j, k, l) = c_IP

Check warning on line 244 in src/simulation/m_ibm.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_ibm.fpp#L244 

Added line #L244 was not covered by tests
end if

! Set Energy
if (bubbles_euler) then
q_cons_vf(E_idx)%sf(j, k, l) = (1 - alpha_IP(1))*(gamma*pres_IP + pi_inf + dyn_pres)
Expand Down Expand Up @@ -309,6 +319,10 @@
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
end do

if (surface_tension) then
q_prim_vf(c_idx)%sf(j, k, l) = c_IP

Check warning on line 323 in src/simulation/m_ibm.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_ibm.fpp#L322-L323 

Added lines #L322 - L323 were not covered by tests
end if

call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
alpha_rho_IP, Re_K, j, k, l)

Expand Down Expand Up @@ -725,16 +739,18 @@

!> Function that uses the interpolation coefficients and the current state
!! at the cell centers in order to estimate the state at the image point
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)
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)

Check warning on line 742 in src/simulation/m_ibm.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_ibm.fpp#L742 

Added line #L742 was not covered by tests
!$acc routine seq
type(scalar_field), &
dimension(sys_size), &
intent(IN) :: q_prim_vf !< Primitive Variables

real(wp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(INOUT) :: pb, mv

type(ghost_point), intent(IN) :: gp
real(wp), intent(INOUT) :: pres_IP
real(wp), dimension(3), intent(INOUT) :: vel_IP
real(wp), intent(INOUT) :: c_IP
real(wp), dimension(num_fluids), intent(INOUT) :: alpha_IP, alpha_rho_IP
real(wp), optional, dimension(:), intent(INOUT) :: r_IP, v_IP, pb_IP, mv_IP
real(wp), optional, dimension(:), intent(INOUT) :: nmom_IP
Expand All @@ -758,6 +774,8 @@
pres_IP = 0._wp
vel_IP = 0._wp

if (surface_tension) c_IP = 0._wp
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should c_IP be optional? if surface tension is optional, just like bubbles_euler

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

good to have I guess even though it pertains to only surface tension. I added it to the ghost point variables for now.


if (bubbles_euler) then
r_IP = 0._wp
v_IP = 0._wp
Expand Down Expand Up @@ -801,6 +819,10 @@
q_prim_vf(advxb + l - 1)%sf(i, j, k)
end do

if (surface_tension) then
c_IP = c_IP + coeff*q_prim_vf(c_idx)%sf(i, j, k)

Check warning on line 823 in src/simulation/m_ibm.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_ibm.fpp#L823 

Added line #L823 was not covered by tests
end if

if (bubbles_euler .and. .not. qbmm) then
!$acc loop seq
do l = 1, nb
Expand Down
11 changes: 10 additions & 1 deletion src/simulation/m_riemann_solvers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1650,7 +1650,7 @@
end do
end if

! SURFACE TENSION FLUX. need to check
! COLOR FUNCTION FLUX
if (surface_tension) then
flux_rs${XYZ}$_vf(j, k, l, c_idx) = &
(xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, c_idx) + &
Expand Down Expand Up @@ -2831,6 +2831,15 @@
s_P*(xi_R - 1._wp))
end do

! COLOR FUNCTION FLUX
if (surface_tension) then
flux_rs${XYZ}$_vf(j, k, l, c_idx) = &
xi_M*qL_prim_rs${XYZ}$_vf(j, k, l, c_idx) &

Check warning on line 2837 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L2836-L2837 

Added lines #L2836 - L2837 were not covered by tests
*(vel_L(idx1) + s_M*(xi_L - 1._wp)) &
+ xi_P*qR_prim_rs${XYZ}$_vf(j + 1, k, l, c_idx) &
*(vel_R(idx1) + s_P*(xi_R - 1._wp))

Check warning on line 2840 in src/simulation/m_riemann_solvers.fpp

View check run for this annotation

Codecov / codecov/patch

src/simulation/m_riemann_solvers.fpp#L2839-L2840 

Added lines #L2839 - L2840 were not covered by tests
end if

! REFERENCE MAP FLUX.
if (hyperelasticity) then
!$acc loop seq
Expand Down
Loading