Skip to content

removing duplicated codes in mixture rules for EE bubbles #1071

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 11 commits into from
Dec 1, 2025
Merged

removing duplicated codes in mixture rules for EE bubbles #1071

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
5d977c3
integrating mixture routines
Nov 29, 2025
ff8ff17
Merge branch 'MFlowCode:master' into bubnorm
hyeoksu-lee Nov 29, 2025
ae827a5
update mixture rules
Nov 29, 2025
e1e953c
fix mixture rule
Nov 29, 2025
50f78b8
correction for vf of EE bubble
Nov 29, 2025
4461db2
Merge branch 'MFlowCode:master' into bubnorm
hyeoksu-lee Nov 29, 2025
1568228
add mpp_lim to mixture rule acc
Nov 29, 2025
4b7ab1c
add comments
Nov 30, 2025
262e980
Merge branch 'bubnorm' into merge
hyeoksu-lee Nov 30, 2025
570e53c
Merge pull request #24 from hyeoksu-lee/merge
hyeoksu-lee Nov 30, 2025
f882347
update bubble case validation
Nov 30, 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
400 changes: 89 additions & 311 deletions src/common/m_variables_conversion.fpp
View file
Open in desktop
Copy link
Contributor

@qodo-code-review qodo-code-review bot Nov 30, 2025

Choose a reason for hiding this comment

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

High-level Suggestion

The new s_compute_species_fraction subroutine has a confusing structure with a conditional override at the end. Refactor it to have clear, separate logic paths for each physical model to improve readability. [High-level, importance: 7]

Solution Walkthrough:

Before:

subroutine s_compute_species_fraction(...)
    if (num_fluids == 1) then
        alpha_rho_K(1) = ...
        if (igr .or. bubbles_euler) then
            alpha_K(1) = 1._wp
        else
            alpha_K(1) = ...
        end if
    else
        ! ... logic for num_fluids > 1
    end if

    if (mpp_lim) then
        ! ... apply limits and normalize alpha_K
    end if

    ! Override for a specific case
    if (num_fluids == 1 .and. bubbles_euler) then
        alpha_K(1) = q_vf(advxb)%sf(k, l, r)
    end if
end subroutine

After:

subroutine s_compute_species_fraction(...)
    if (num_fluids == 1 .and. bubbles_euler) then
        alpha_rho_K(1) = q_vf(contxb)%sf(k, l, r)
        alpha_K(1) = q_vf(advxb)%sf(k, l, r)
    else if (num_fluids == 1) then
        alpha_rho_K(1) = ...
        if (igr) then
            alpha_K(1) = 1._wp
        else
            alpha_K(1) = ...
        end if
    else
        ! ... logic for num_fluids > 1
    end if

    if (mpp_lim) then
        ! ... apply limits and normalize alpha_K
    end if
end subroutine

Large diffs are not rendered by default.

26 changes: 9 additions & 17 deletions src/simulation/m_bubbles_EE.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -242,34 +242,26 @@ contains
myalpha(ii) = q_cons_vf(advxb + ii - 1)%sf(j, k, l)
end do

myRho = 0._wp
n_tait = 0._wp
B_tait = 0._wp
if (num_fluids == 1) then
myRho = myalpha_rho(1)
n_tait = gammas(1)
B_tait = pi_infs(1)/pi_fac
else
myRho = 0._wp
n_tait = 0._wp
B_tait = 0._wp

if (mpp_lim .and. (num_fluids > 2)) then
$:GPU_LOOP(parallelism='[seq]')
do ii = 1, num_fluids
myRho = myRho + myalpha_rho(ii)
n_tait = n_tait + myalpha(ii)*gammas(ii)
B_tait = B_tait + myalpha(ii)*pi_infs(ii)
end do
else if (num_fluids > 2) then
$:GPU_LOOP(parallelism='[seq]')
do ii = 1, num_fluids - 1
myRho = myRho + myalpha_rho(ii)
n_tait = n_tait + myalpha(ii)*gammas(ii)
B_tait = B_tait + myalpha(ii)*pi_infs(ii)
B_tait = B_tait + myalpha(ii)*pi_infs(ii)/pi_fac
end do
else
myRho = myalpha_rho(1)
n_tait = gammas(1)
B_tait = pi_infs(1)/pi_fac
end if

n_tait = 1._wp/n_tait + 1._wp !make this the usual little 'gamma'
B_tait = B_tait*(n_tait - 1)/n_tait ! make this the usual pi_inf
Comment on lines 262 to 263
Copy link
Contributor

@qodo-code-review qodo-code-review bot Nov 30, 2025

Choose a reason for hiding this comment

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

Suggestion: Prevent a potential division-by-zero error in the calculation of n_tait and B_tait. Use max(n_tait, sgm_eps) in the denominator to ensure it is non-zero. [possible issue, importance: 8]

Suggested change
n_tait = 1._wp/n_tait + 1._wp !make this the usual little 'gamma'
B_tait = B_tait*(n_tait - 1)/n_tait ! make this the usual pi_inf
n_tait = 1._wp/max(n_tait, sgm_eps) + 1._wp !make this the usual little 'gamma'
B_tait = B_tait*(n_tait - 1)/max(n_tait, sgm_eps) ! make this the usual pi_inf


myRho = q_prim_vf(1)%sf(j, k, l)
myP = q_prim_vf(E_idx)%sf(j, k, l)
alf = q_prim_vf(alf_idx)%sf(j, k, l)
myR = q_prim_vf(rs(q))%sf(j, k, l)
Expand Down
6 changes: 1 addition & 5 deletions src/simulation/m_bubbles_EL.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -642,11 +642,7 @@ contains
call s_get_pinf(k, q_prim_vf, 1, myPinf, cell, aux1, aux2)

! Obtain liquid density and computing speed of sound from pinf
$:GPU_LOOP(parallelism='[seq]')
do i = 1, num_fluids
myalpha_rho(i) = q_prim_vf(i)%sf(cell(1), cell(2), cell(3))
myalpha(i) = q_prim_vf(E_idx + i)%sf(cell(1), cell(2), cell(3))
end do
call s_compute_species_fraction(q_prim_vf, cell(1), cell(2), cell(3), myalpha_rho, myalpha)
call s_convert_species_to_mixture_variables_acc(myRho, gamma, pi_inf, qv, myalpha, &
myalpha_rho, Re)
call s_compute_cson_from_pinf(q_prim_vf, myPinf, cell, myRho, gamma, pi_inf, myCson)
Expand Down
6 changes: 1 addition & 5 deletions src/simulation/m_cbc.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -813,11 +813,7 @@ contains
adv_local(i) = q_prim_rs${XYZ}$_vf(0, k, r, E_idx + i)
end do

if (bubbles_euler) then
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc)
else
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc)
end if
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, adv_local, alpha_rho, Re_cbc)

$:GPU_LOOP(parallelism='[seq]')
do i = 1, contxe
Expand Down
8 changes: 3 additions & 5 deletions src/simulation/m_hyperelastic.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,11 +110,9 @@ contains
do l = 0, p
do k = 0, n
do j = 0, m
$:GPU_LOOP(parallelism='[seq]')
do i = 1, num_fluids
alpha_rho_k(i) = q_cons_vf(i)%sf(j, k, l)
alpha_k(i) = q_cons_vf(advxb + i - 1)%sf(j, k, l)
end do

call s_compute_species_fraction(q_cons_vf, j, k, l, alpha_rho_k, alpha_k)

! If in simulation, use acc mixture subroutines
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha_k, &
alpha_rho_k, Re, G_local, Gs_hyper)
Expand Down
3 changes: 0 additions & 3 deletions src/simulation/m_ibm.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,9 +249,6 @@ contains
if (elasticity) then
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
alpha_rho_IP, Re_K, G_K, Gs)
else if (bubbles_euler) then
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
alpha_rho_IP, Re_K)
else
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv_K, alpha_IP, &
alpha_rho_IP, Re_K)
Expand Down
24 changes: 1 addition & 23 deletions src/simulation/m_sim_helpers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,33 +109,11 @@ contains

integer :: i

if (igr) then
if (num_fluids == 1) then
alpha_rho(1) = q_prim_vf(contxb)%sf(j, k, l)
alpha(1) = 1._wp
else
$:GPU_LOOP(parallelism='[seq]')
do i = 1, num_fluids - 1
alpha_rho(i) = q_prim_vf(i)%sf(j, k, l)
alpha(i) = q_prim_vf(advxb + i - 1)%sf(j, k, l)
end do

alpha_rho(num_fluids) = q_prim_vf(num_fluids)%sf(j, k, l)
alpha(num_fluids) = 1._wp - sum(alpha(1:num_fluids - 1))
end if
else
$:GPU_LOOP(parallelism='[seq]')
do i = 1, num_fluids
alpha_rho(i) = q_prim_vf(i)%sf(j, k, l)
alpha(i) = q_prim_vf(advxb + i - 1)%sf(j, k, l)
end do
end if
call s_compute_species_fraction(q_prim_vf, j, k, l, alpha_rho, alpha)

if (elasticity) then
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, &
alpha_rho, Re, G_local, Gs)
elseif (bubbles_euler) then
call s_convert_species_to_mixture_variables_bubbles_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re)
else
call s_convert_species_to_mixture_variables_acc(rho, gamma, pi_inf, qv, alpha, alpha_rho, Re)
end if
Expand Down
4 changes: 2 additions & 2 deletions toolchain/mfc/case_validator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -369,8 +369,8 @@ def check_stiffened_eos(self):
if num_fluids is None:
return

# Allow one extra fluid property slot when using bubbles_euler with single-component flows
bub_fac = 1 if (bubbles_euler and num_fluids == 1) else 0
# Allow one extra fluid property slot when using bubbles_euler
bub_fac = 1 if (bubbles_euler) else 0

for i in range(1, num_fluids + 1 + bub_fac):
gamma = self.get(f'fluid_pp({i})%gamma')
Expand Down
Loading