Skip to content

Fix the 6-equation phase energies at characteristic boundaries and the wall contact angle with surface tension #1954

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

Draft
wilfonba wants to merge 5 commits into
base: master
Choose a base branch
from
Draft

Fix the 6-equation phase energies at characteristic boundaries and the wall contact angle with surface tension #1954

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
17592f3
fixes for surface tension at walls and six-eqn CBCs
wilfonba Oct 7, 2026
8f8e510
precheck and tested on GPUs
wilfonba Oct 7, 2026
97035b7
update tolerance for contact line cases
wilfonba Oct 8, 2026
06bef60
Merge branch 'master' into CBCContactLineFixes
wilfonba Oct 9, 2026
f54c64c
Merge branch 'master' into CBCContactLineFixes
wilfonba Oct 14, 2026
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
12 changes: 10 additions & 2 deletions src/common/m_boundary_common.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -176,6 +176,12 @@ contains
call s_dirichlet(q_prim_vf, bc_dir, bc_loc, k, l, q_T_sf)
end select

if (surface_tension .and. (bc_code == BC_NO_SLIP_WALL .or. bc_code == BC_SLIP_WALL)) then
call s_mirror_into_ghost_cells(q_prim_vf(eqn_idx%c), bc_dir, bc_loc, k, l, sgn=1._wp, off=1)
else if (surface_tension .and. bc_code == BC_CHAR_SLIP_WALL) then
call s_mirror_into_ghost_cells(q_prim_vf(eqn_idx%c), bc_dir, bc_loc, k, l, sgn=1._wp, off=0)
end if

if (qbmm .and. (.not. polytropic) .and. present(pb_in) .and. present(mv_in) .and. (bc_code <= BC_GHOST_EXTRAP) &
& .and. .not. (bc_dir == 2 .and. bc_loc == -1 .and. bc_code == BC_AXIS)) then
call s_qbmm_extrapolation(bc_dir, bc_loc, k, l, pb_in, mv_in)
Expand Down Expand Up @@ -255,8 +261,10 @@ contains
select case (bc_code)
case (BC_PERIODIC)
call s_color_function_periodic(c_divs, bc_dir, bc_loc, k, l)
case (BC_REFLECTIVE)
call s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l)
case (BC_REFLECTIVE, BC_SLIP_WALL, BC_NO_SLIP_WALL)
call s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l, 1)
case (BC_CHAR_SLIP_WALL)
call s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l, 0)
case default
call s_color_function_ghost_cell_extrapolation(c_divs, bc_dir, bc_loc, k, l)
end select
Expand Down
116 changes: 45 additions & 71 deletions src/common/m_boundary_primitives.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1140,82 +1140,56 @@ contains

end subroutine s_color_function_periodic

!> Apply reflective boundary conditions to the color function and its divergence fields.
subroutine s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l)
!> Mirror a field into the ghost cells, times sgn: about the domain face (off = 1, a symmetry plane or ghost-cell wall) or about
!! the boundary cell's center (off = 0, where a characteristic wall sits)
subroutine s_mirror_into_ghost_cells(sf, bc_dir, bc_loc, k, l, sgn, off)

$:GPU_ROUTINE(function_name='s_mirror_into_ghost_cells', parallelism='[seq]', cray_inline=True)
type(scalar_field), intent(inout) :: sf
integer, intent(in) :: bc_dir, bc_loc, k, l, off
real(wp), intent(in) :: sgn
integer :: j

do j = 1, buff_size
if (bc_dir == 1) then
if (bc_loc == -1) then
sf%sf(-j, k, l) = sgn*sf%sf(j - off, k, l)
else
sf%sf(m + j, k, l) = sgn*sf%sf(m - j + off, k, l)
end if
else if (bc_dir == 2) then
if (bc_loc == -1) then
sf%sf(k, -j, l) = sgn*sf%sf(k, j - off, l)
else
sf%sf(k, n + j, l) = sgn*sf%sf(k, n - j + off, l)
end if
else
if (bc_loc == -1) then
sf%sf(k, l, -j) = sgn*sf%sf(k, l, j - off)
else
sf%sf(k, l, p + j) = sgn*sf%sf(k, l, p - j + off)
end if
end if
end do

end subroutine s_mirror_into_ghost_cells

!> Reflect the color function's gradient and its magnitude into the ghost cells (a 90 degree contact angle): the wall-normal
!! component changes sign. off as in s_mirror_into_ghost_cells.
subroutine s_color_function_reflective(c_divs, bc_dir, bc_loc, k, l, off)

$:GPU_ROUTINE(function_name='s_color_function_reflective', parallelism='[seq]', cray_inline=True)
type(scalar_field), dimension(num_dims + 1), intent(inout) :: c_divs
integer, intent(in) :: bc_dir, bc_loc
integer, intent(in) :: k, l
integer :: j, i
integer, intent(in) :: bc_dir, bc_loc, k, l, off
integer :: i

if (bc_dir == 1) then !< x-direction
if (bc_loc == -1) then ! bc_x%beg
do i = 1, num_dims + 1
do j = 1, buff_size
if (i == bc_dir) then
c_divs(i)%sf(-j, k, l) = -c_divs(i)%sf(j - 1, k, l)
else
c_divs(i)%sf(-j, k, l) = c_divs(i)%sf(j - 1, k, l)
end if
end do
end do
else !< bc_x%end
do i = 1, num_dims + 1
do j = 1, buff_size
if (i == bc_dir) then
c_divs(i)%sf(m + j, k, l) = -c_divs(i)%sf(m - (j - 1), k, l)
else
c_divs(i)%sf(m + j, k, l) = c_divs(i)%sf(m - (j - 1), k, l)
end if
end do
end do
end if
else if (bc_dir == 2) then !< y-direction
if (bc_loc == -1) then !< bc_y%beg
do i = 1, num_dims + 1
do j = 1, buff_size
if (i == bc_dir) then
c_divs(i)%sf(k, -j, l) = -c_divs(i)%sf(k, j - 1, l)
else
c_divs(i)%sf(k, -j, l) = c_divs(i)%sf(k, j - 1, l)
end if
end do
end do
else !< bc_y%end
do i = 1, num_dims + 1
do j = 1, buff_size
if (i == bc_dir) then
c_divs(i)%sf(k, n + j, l) = -c_divs(i)%sf(k, n - (j - 1), l)
else
c_divs(i)%sf(k, n + j, l) = c_divs(i)%sf(k, n - (j - 1), l)
end if
end do
end do
end if
else if (bc_dir == 3) then !< z-direction
if (bc_loc == -1) then !< bc_z%beg
do i = 1, num_dims + 1
do j = 1, buff_size
if (i == bc_dir) then
c_divs(i)%sf(k, l, -j) = -c_divs(i)%sf(k, l, j - 1)
else
c_divs(i)%sf(k, l, -j) = c_divs(i)%sf(k, l, j - 1)
end if
end do
end do
else !< bc_z%end
do i = 1, num_dims + 1
do j = 1, buff_size
if (i == bc_dir) then
c_divs(i)%sf(k, l, p + j) = -c_divs(i)%sf(k, l, p - (j - 1))
else
c_divs(i)%sf(k, l, p + j) = c_divs(i)%sf(k, l, p - (j - 1))
end if
end do
end do
do i = 1, num_dims + 1
if (i == bc_dir) then
call s_mirror_into_ghost_cells(c_divs(i), bc_dir, bc_loc, k, l, sgn=-1._wp, off=off)
else
call s_mirror_into_ghost_cells(c_divs(i), bc_dir, bc_loc, k, l, sgn=1._wp, off=off)
end if
end if
end do

end subroutine s_color_function_reflective

Expand Down
18 changes: 15 additions & 3 deletions src/common/m_variables_conversion.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -967,7 +967,8 @@ contains
real(wp) :: G_K
real(wp) :: blkmod1_K, blkmod2_K, K_K
real(wp) :: T_K, mix_mol_weight, R_gas
integer :: i, j, k, l !< Generic loop iterators
real(wp) :: alpha_i, alpha_rho_i, e_i !< 6-eq phase scalars
integer :: i, j, k, l !< Generic loop iterators

is1b = is1%beg; is1e = is1%end
is2b = is2%beg; is2e = is2%end
Expand All @@ -978,8 +979,8 @@ contains
! Computing the flux variables from the primitive variables, without accounting for the contribution of either viscosity or
! capillarity
$:GPU_PARALLEL_LOOP(collapse=3, private='[alpha_rho_K, vel_K, alpha_K, Re_K, Y_K, rho_K, vel_K_sum, pres_K, E_K, gamma_K, &
& pi_inf_K, qv_K, G_K, blkmod1_K, blkmod2_K, K_K, T_K, mix_mol_weight, R_gas]', copyin='[dir_idx_in, &
& dir_flg_in, hll_u_interface_in]')
& pi_inf_K, qv_K, G_K, blkmod1_K, blkmod2_K, K_K, T_K, mix_mol_weight, R_gas, alpha_i, alpha_rho_i, &
& e_i]', copyin='[dir_idx_in, dir_flg_in, hll_u_interface_in]')
do l = is3b, is3e
do k = is2b, is2e
do j = is1b, is1e
Expand Down Expand Up @@ -1047,6 +1048,17 @@ contains
! energy flux, u(E+p)
FK_vf(j, k, l, eqn_idx%E) = vel_K(dir_idx_in(1))*(E_K + pres_K)

! 6-eq phase internal energy fluxes, u (alpha_i (Gamma_i p + Pi_i) + alpha_rho_i qv_i)
if (model_eqns == model_eqns_6eq) then
$:GPU_LOOP(parallelism='[seq]')
do i = 1, num_fluids
alpha_i = alpha_K(i)
alpha_rho_i = alpha_rho_K(i)
call s_phase_internal_energy(pres_K, alpha_i, alpha_rho_i, i, e_i)
FK_vf(j, k, l, eqn_idx%int_en%beg + i - 1) = vel_K(dir_idx_in(1))*e_i
end do
end if

! Species advection Flux, \rho*u*Y
if (chemistry) then
$:GPU_LOOP(parallelism='[seq]')
Expand Down
40 changes: 31 additions & 9 deletions src/simulation/m_cbc.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,9 +12,10 @@ module m_cbc
use m_derived_types
use m_global_parameters
use m_variables_conversion
use m_eos
use m_compute_cbc
use m_boundary_primitives, only: f_vel_ramp
use m_constants, only: riemann_solver_hll, model_eqns_gamma_law, recon_type_weno, recon_type_muscl
use m_constants, only: riemann_solver_hll, model_eqns_gamma_law, model_eqns_6eq, recon_type_weno, recon_type_muscl
use m_thermochem, only: gas_constant, get_mixture_molecular_weight, get_species_enthalpies_rt, molecular_weights, &
& get_mole_fractions
use m_thermochem_state, only: get_mixture_caloric_state
Expand Down Expand Up @@ -94,6 +95,8 @@ contains

if (chemistry) then
flux_cbc_index = sys_size
else if (model_eqns == model_eqns_6eq) then
flux_cbc_index = eqn_idx%int_en%end ! the phase internal energies follow the waves too
else
flux_cbc_index = eqn_idx%adv%end
end if
Expand Down Expand Up @@ -500,18 +503,20 @@ contains
real(wp), dimension(${NUM_SPECIES}$) :: Ys, h_k, dYs_dt, dYs_ds, Xs, Gamma_i, Cp_i
real(wp), dimension(2) :: Re_cbc
real(wp), dimension(3) :: lambda
real(wp) :: rho !< Cell averaged density
real(wp) :: pres !< Cell averaged pressure
real(wp) :: E !< Cell averaged energy
real(wp) :: gamma !< Cell averaged specific heat ratio
real(wp) :: pi_inf !< Cell averaged liquid stiffness
real(wp) :: qv !< Cell averaged fluid reference energy
real(wp) :: rho !< Cell averaged density
real(wp) :: pres !< Cell averaged pressure
real(wp) :: E !< Cell averaged energy
real(wp) :: gamma !< Cell averaged specific heat ratio
real(wp) :: pi_inf !< Cell averaged liquid stiffness
real(wp) :: qv !< Cell averaged fluid reference energy
real(wp) :: c
real(wp) :: Ma
real(wp) :: T, sum_Enthalpies
real(wp) :: Cv, Cp, e_mix, Mw, R_gas
real(wp) :: vel_K_sum, vel_dv_dt_sum
integer :: i, j, k, r !< Generic loop iterators
real(wp) :: alpha_rho_k, alpha_k !< state of phase k (6-eq)
real(wp) :: rho_k, gamma_k, pi_inf_k, dpi_k, dgamma_k !< coefficients of phase k (6-eq)
integer :: i, j, k, r !< Generic loop iterators
! Reshaping of inputted data and association of the FD and PI coefficients, or CBC coefficients, respectively, hinging on
! selected CBC coordinate direction

Expand Down Expand Up @@ -598,7 +603,8 @@ contains
& dalpha_rho_ds, dpres_ds, dvel_dt, dadv_dt, dalpha_rho_dt, L, lambda, Ys, dYs_dt, dYs_ds, &
& h_k, Cp_i, Gamma_i, Xs, drho_dt, dpres_dt, dpi_inf_dt, dqv_dt, dgamma_dt, rho, pres, E, &
& gamma, pi_inf, qv, c, Ma, T, sum_Enthalpies, Cv, Cp, e_mix, Mw, R_gas, vel_K_sum, &
& vel_dv_dt_sum, i, j, ramp]', copyin='[dir_idx]')
& vel_dv_dt_sum, i, j, ramp, alpha_rho_k, alpha_k, rho_k, gamma_k, pi_inf_k, dpi_k, &
& dgamma_k]', copyin='[dir_idx]')
do r = is3%beg, is3%end
do k = is2%beg, is2%end
! Ramp factor for a smoothly starting inflow, evaluated here from mytime rather than
Expand Down Expand Up @@ -878,6 +884,22 @@ contains
& + rho*vel_dv_dt_sum + 5.e-1_wp*drho_dt*vel_K_sum)
end if

! 6-eq: each phase's energy alpha (Gamma p + Pi) + alpha_rho qv changes at the same rates, so the phase
! energies keep summing to the mixture's. The boundary face's interface velocity is the interior one, so
! the alpha p div(u) source adds nothing here.
if (model_eqns == model_eqns_6eq) then
$:GPU_LOOP(parallelism='[seq]')
do i = 1, num_fluids
alpha_rho_k = alpha_rho(i)
alpha_k = adv_local(i)
call s_phase_coefficients(alpha_rho_k, alpha_k, i, rho_k, gamma_k, pi_inf_k, dpi_k, dgamma_k)
flux_rs${XYZ}$_vf_l(-1, k, r, eqn_idx%int_en%beg + i - 1) = flux_rs${XYZ}$_vf_l(0, k, r, &
& eqn_idx%int_en%beg + i - 1) + ds(0)*(gamma_k*(adv_local(i)*dpres_dt &
& + pres*dadv_dt(i)) + pi_inf_k*dadv_dt(i) + (pres*dgamma_k + dpi_k) &
& *(dalpha_rho_dt(i) - rho_k*dadv_dt(i)) + qvs(i)*dalpha_rho_dt(i))
end do
end if

! Only HLL Method 1 uses per-fluid alpha source traces. HLL Method 2 carries a shared interface velocity and
! must follow the same CBC representation as HLLC.
if (riemann_solver == riemann_solver_hll .and. .not. hll_u_interface) then
Expand Down
Loading
Loading