Skip to content

Immersed boundaries integration with IGR #1095

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

Open
Cowsreal wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Immersed boundaries integration with IGR #1095

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
af75f84
Integrate IBM with IGR
Dec 18, 2025
1340497
Add comments
Dec 19, 2025
b07045b
Minor fixes
Dec 19, 2025
489bc7c
Change from using jac_sf to jac
Dec 19, 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
250 changes: 196 additions & 54 deletions src/simulation/m_ibm.fpp
View file
Open in desktop
Copy link
Contributor

@qodo-code-review qodo-code-review bot Dec 19, 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 implementation for moving immersed boundaries with IGR is incomplete, as indicated by a !TEMPORARY code comment. This should be either fully implemented or disabled with a runtime check to avoid incorrect simulations. [High-level, importance: 9]

Solution Walkthrough:

Before:

! in s_ibm_correct_state
...
! !TEMPORARY, NEED TO FIX FOR MOVING IGR
if (.not. igr) then
    if (patch_ib(patch_id)%moving_ibm <= 1) then
        q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
    else
        ! Logic for moving IBM pressure calculation
        ...
    end if
end if
! When igr is true, pressure is not set for the ghost point.
...

After:

! in s_ibm_correct_state
...
! set the pressure
if (igr .and. patch_ib(patch_id)%moving_ibm > 1) then
    ! Option A: Add a runtime error to prevent invalid simulations
    call s_throw_error("Moving IBM with IGR is not yet supported.")
else
    ! Option B: Fully implement the feature
    if (patch_ib(patch_id)%moving_ibm <= 1) then
        q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
    else
        ! Correctly implement moving IBM pressure calculation for IGR
        ...
    end if
end if
...

Copy link
Contributor

@qodo-code-review qodo-code-review bot Dec 19, 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 PR currently omits pressure calculations for moving immersed boundaries when IGR is enabled, marked by a !TEMPORARY comment. This functionality should be implemented to complete the feature. [High-level, importance: 8]

Solution Walkthrough:

Before:

subroutine s_ibm_correct_state(...)
    ...
    ! !TEMPORARY, NEED TO FIX FOR MOVING IGR
    if (.not. igr) then
        if (patch_ib(patch_id)%moving_ibm <= 1) then
            ! Set pressure for static IBM
            q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
        else
            ! Set pressure for moving IBM
            ...
        end if
    end if
    ...
    ! For IGR, energy is set later from interpolated values,
    ! but without the moving boundary pressure correction.
    q_cons_vf(E_idx)%sf(j, k, l) = gamma*pres_IP + pi_inf + dyn_pres
    ...
end subroutine

After:

subroutine s_ibm_correct_state(...)
    ...
    if (.not. igr) then
        if (patch_ib(patch_id)%moving_ibm <= 1) then
            q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
        else
            ! Set pressure for moving IBM
            ...
        end if
    else ! Handle IGR case
        if (patch_ib(patch_id)%moving_ibm > 1) then
            ! Implement moving boundary pressure correction for IGR
            ! This might involve adjusting pres_IP before it's used
            ! to calculate the energy for the ghost point.
            pres_IP = pres_IP / (1._wp - ... )
        end if
    end if
    ...
    q_cons_vf(E_idx)%sf(j, k, l) = gamma*pres_IP + pi_inf + dyn_pres
    ...
end subroutine

Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,51 @@ contains

end subroutine s_populate_ib_buffers

!> Interpolates sigma from the m_igr module at all ghost points
!! @param jac: Sigma, Entropic pressure
subroutine s_interpolate_sigma_igr(jac)
real(wp), dimension(:, :, :), intent(inout) :: jac
integer :: j, k, l, r, s, t, i
integer :: j1, j2, k1, k2, l1, l2
real(wp) :: coeff, jac_IP
type(ghost_point) :: gp

! At all ghost points, use its image point to interpolate sigma
if (num_gps > 0) then
$:GPU_PARALLEL_LOOP(private='[i, j, k, l, j1, j2, k1, k2, l1, l2, r, s, t, gp, coeff, jac_IP]')
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Dec 19, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

P2: GPU data mapping concern: jac_sf is accessed within the GPU parallel kernel but is not explicitly included in data clauses (copyin/copy/update). This could result in accessing host memory from device code or using stale data. Verify that jac_sf%sf is properly mapped to the device before this kernel executes.

Prompt for AI agents 
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_ibm.fpp, line 162:

<comment>GPU data mapping concern: `jac_sf` is accessed within the GPU parallel kernel but is not explicitly included in data clauses (copyin/copy/update). This could result in accessing host memory from device code or using stale data. Verify that `jac_sf%sf` is properly mapped to the device before this kernel executes.</comment>

<file context>
@@ -150,6 +150,44 @@ contains
+
+        ! At all ghost points, use its image point to interpolate sigma
+        if (num_gps &gt; 0) then
+            $:GPU_PARALLEL_LOOP(private=&#39;[i, j, k, l, j1, j2, k1, k2, l1, l2, r, s, t, gp, coeff, jac_IP]&#39;)
+            do i = 1, num_gps
+                jac_IP = 0._wp
</file context>

✅ Addressed in 489bc7c

do i = 1, num_gps
jac_IP = 0._wp
gp = ghost_points(i)
r = gp%loc(1)
s = gp%loc(2)
t = gp%loc(3)

j1 = gp%ip_grid(1); j2 = j1 + 1
k1 = gp%ip_grid(2); k2 = k1 + 1
l1 = gp%ip_grid(3); l2 = l1 + 1

if (p == 0) then
l1 = 0
l2 = 0
end if
Comment on lines +176 to +179
Copy link
Contributor

@qodo-code-review qodo-code-review bot Dec 19, 2025

Choose a reason for hiding this comment

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

Suggestion: In s_update_igr, import the global parameter p from m_global_parameters to fix a compilation error. [possible issue, importance: 9]

Suggested change
if (p == 0) then
l1 = 0
l2 = 0
end if
subroutine s_update_igr(jac_sf)
use m_global_parameters, only: p, wp
type(scalar_field), dimension(1), intent(inout) :: jac_sf
integer :: j, k, l, r, s, t, i
integer :: j1, j2, k1, k2, l1, l2
real(wp) :: coeff, jac_IP
type(ghost_point) :: gp
! At all ghost points, use its image point to interpolate sigma
if (num_gps > 0) then
$:GPU_PARALLEL_LOOP(private='[i, j, k, l, j1, j2, k1, k2, l1, l2, r, s, t, gp, coeff, jac_IP]')
do i = 1, num_gps
jac_IP = 0._wp
gp = ghost_points(i)
r = gp%loc(1)
s = gp%loc(2)
t = gp%loc(3)
j1 = gp%ip_grid(1); j2 = j1 + 1
k1 = gp%ip_grid(2); k2 = k1 + 1
l1 = gp%ip_grid(3); l2 = l1 + 1
if (p == 0) then
l2 = l1
end if
$:GPU_LOOP(parallelism='[seq]')
do l = l1, l2
$:GPU_LOOP(parallelism='[seq]')
do k = k1, k2
$:GPU_LOOP(parallelism='[seq]')
do j = j1, j2
coeff = gp%interp_coeffs(j - j1 + 1, k - k1 + 1, l - l1 + 1)
jac_IP = jac_IP + coeff*jac_sf(1)%sf(j, k, l)
end do
end do
end do
jac_sf(1)%sf(r, s, t) = jac_IP
end do
$:END_GPU_PARALLEL_LOOP()
end if
end subroutine s_update_igr


$:GPU_LOOP(parallelism='[seq]')
do l = l1, l2
$:GPU_LOOP(parallelism='[seq]')
do k = k1, k2
$:GPU_LOOP(parallelism='[seq]')
do j = j1, j2
coeff = gp%interp_coeffs(j - j1 + 1, k - k1 + 1, l - l1 + 1)
jac_IP = jac_IP + coeff*jac(j, k, l)
end do
end do
end do
jac(r, s, t) = jac_IP
end do
$:END_GPU_PARALLEL_LOOP()
end if
end subroutine s_interpolate_sigma_igr

!> Subroutine that updates the conservative variables at the ghost points
!! @param q_cons_vf Conservative Variables
!! @param q_prim_vf Primitive variables
Expand All @@ -159,7 +204,7 @@ contains

type(scalar_field), &
dimension(sys_size), &
intent(INOUT) :: q_cons_vf !< Primitive Variables
intent(INOUT) :: q_cons_vf !< Conservative Variables

type(scalar_field), &
dimension(sys_size), &
Expand Down Expand Up @@ -196,18 +241,20 @@ contains
type(ghost_point) :: gp
type(ghost_point) :: innerp

! set the Moving IBM interior Pressure Values
$:GPU_PARALLEL_LOOP(private='[i,j,k]', copyin='[E_idx]', collapse=3)
do l = 0, p
do k = 0, n
do j = 0, m
if (ib_markers%sf(j, k, l) /= 0) then
q_prim_vf(E_idx)%sf(j, k, l) = 1._wp
end if
if (.not. igr) then
! set the Moving IBM interior Pressure Values
$:GPU_PARALLEL_LOOP(private='[j,k,l]', copyin='[E_idx]', collapse=3)
do l = 0, p
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Dec 19, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

P2: Inconsistent loop bounds: zeroing jac for IB markers uses 0..m, 0..n, 0..p ranges, but copying to jac_old uses idwbuff(...)%beg:...%end. If idwbuff extends beyond the computational domain (e.g., includes ghost/halo cells), this may leave ghost regions of jac_old inconsistent after IB updates.

Prompt for AI agents 
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_ibm.fpp, line 240:

<comment>Inconsistent loop bounds: zeroing `jac` for IB markers uses `0..m, 0..n, 0..p` ranges, but copying to `jac_old` uses `idwbuff(...)%beg:...%end`. If `idwbuff` extends beyond the computational domain (e.g., includes ghost/halo cells), this may leave ghost regions of `jac_old` inconsistent after IB updates.</comment>

<file context>
@@ -196,18 +234,20 @@ contains
+        if (.not. igr) then
+            ! set the Moving IBM interior Pressure Values
+            $:GPU_PARALLEL_LOOP(private=&#39;[i,j,k]&#39;, copyin=&#39;[E_idx]&#39;, collapse=3)
+            do l = 0, p
+                do k = 0, n
+                    do j = 0, m
</file context>
Fix with Cubic

do k = 0, n
do j = 0, m
if (ib_markers%sf(j, k, l) /= 0) then
q_prim_vf(E_idx)%sf(j, k, l) = 1._wp
end if
end do
end do
end do
end do
$:END_GPU_PARALLEL_LOOP()
$:END_GPU_PARALLEL_LOOP()
end if

if (num_gps > 0) then
$:GPU_PARALLEL_LOOP(private='[i,physical_loc,dyn_pres,alpha_rho_IP, alpha_IP,pres_IP,vel_IP,vel_g,vel_norm_IP,r_IP, v_IP,pb_IP,mv_IP,nmom_IP,presb_IP,massv_IP,rho, gamma,pi_inf,Re_K,G_K,Gs,gp,innerp,norm,buf, radial_vector, rotation_velocity, j,k,l,q,qv_K,c_IP,nbub,patch_id]')
Expand Down Expand Up @@ -239,34 +286,52 @@ contains
call 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_in, mv_in, presb_IP, massv_IP)
else if (igr) then
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, q_cons_vf=q_cons_vf)
else
call s_interpolate_image_point(q_prim_vf, gp, &
alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP)
end if

dyn_pres = 0._wp

! Set q_prim_vf params at GP so that mixture vars calculated properly
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids
q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
end do
if (igr) then
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Dec 19, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

P1: The optional q_cons_vf argument is accessed without a present(q_cons_vf) guard when igr is true. If this subroutine is ever called with igr enabled but without passing q_cons_vf, this will cause undefined behavior. Consider adding if (igr .and. .not. present(q_cons_vf)) error handling or restructuring the logic.

Prompt for AI agents 
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_ibm.fpp, line 291:

<comment>The optional `q_cons_vf` argument is accessed without a `present(q_cons_vf)` guard when `igr` is true. If this subroutine is ever called with `igr` enabled but without passing `q_cons_vf`, this will cause undefined behavior. Consider adding `if (igr .and. .not. present(q_cons_vf))` error handling or restructuring the logic.</comment>

<file context>
@@ -239,34 +279,52 @@ contains
-                    q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
-                    q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
-                end do
+                if (igr) then
+                    if (num_fluids == 1) then
+                        q_cons_vf(1)%sf(j, k, l) = alpha_rho_IP(1)
</file context>
Fix with Cubic

if (num_fluids == 1) then
q_cons_vf(1)%sf(j, k, l) = alpha_rho_IP(1)
else
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids - 1
q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
q_cons_vf(E_idx + q)%sf(j, k, l) = alpha_IP(q)
end do
q_cons_vf(num_fluids)%sf(j, k, l) = alpha_rho_IP(num_fluids)
end if
else
! Set q_prim_vf params at GP so that mixture vars calculated properly
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids
q_prim_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
q_prim_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
end do
end if

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

! set the pressure
if (patch_ib(patch_id)%moving_ibm <= 1) then
q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
else
q_prim_vf(E_idx)%sf(j, k, l) = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids
! Se the pressure inside a moving immersed boundary based upon the pressure of the image point. acceleration, and normal vector direction
q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) + pres_IP/(1._wp - 2._wp*abs(levelset%sf(j, k, l, patch_id)*alpha_rho_IP(q)/pres_IP)*dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, levelset_norm%sf(j, k, l, patch_id, :)))
end do
! !TEMPORARY, NEED TO FIX FOR MOVING IGR
if (.not. igr) then
if (patch_ib(patch_id)%moving_ibm <= 1) then
q_prim_vf(E_idx)%sf(j, k, l) = pres_IP
else
q_prim_vf(E_idx)%sf(j, k, l) = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids
! Se the pressure inside a moving immersed boundary based upon the pressure of the image point. acceleration, and normal vector direction
q_prim_vf(E_idx)%sf(j, k, l) = q_prim_vf(E_idx)%sf(j, k, l) + pres_IP/(1._wp - 2._wp*abs(levelset%sf(j, k, l, patch_id)*alpha_rho_IP(q)/pres_IP)*dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, levelset_norm%sf(j, k, l, patch_id, :)))
end do
end if
end if

if (model_eqns /= 4) then
Expand Down Expand Up @@ -322,12 +387,14 @@ contains
vel_g(q - momxb + 1)/2._wp
end do

! Set continuity and adv vars
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids
q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
end do
if (.not. igr) then
! Set continuity and adv vars
$:GPU_LOOP(parallelism='[seq]')
do q = 1, num_fluids
q_cons_vf(q)%sf(j, k, l) = alpha_rho_IP(q)
q_cons_vf(advxb + q - 1)%sf(j, k, l) = alpha_IP(q)
end do
end if

! Set color function
if (surface_tension) then
Expand All @@ -340,6 +407,7 @@ contains
else
q_cons_vf(E_idx)%sf(j, k, l) = gamma*pres_IP + pi_inf + dyn_pres
end if

! Set bubble vars
if (bubbles_euler .and. .not. qbmm) then
call s_comp_n_from_prim(alpha_IP(1), r_IP, nbub, weight)
Expand Down Expand Up @@ -827,11 +895,14 @@ contains
!! 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, c_IP, r_IP, v_IP, pb_IP, &
mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP)
mv_IP, nmom_IP, pb_in, mv_in, presb_IP, massv_IP, q_cons_vf)
$:GPU_ROUTINE(parallelism='[seq]')
type(scalar_field), &
dimension(sys_size), &
intent(IN) :: q_prim_vf !< Primitive Variables
type(scalar_field), optional, &
dimension(sys_size), &
intent(IN) :: q_cons_vf !< Conservative Variables

real(stp), optional, dimension(idwbuff(1)%beg:, idwbuff(2)%beg:, idwbuff(3)%beg:, 1:, 1:), intent(IN) :: pb_in, mv_in

Expand All @@ -847,6 +918,12 @@ contains
integer :: i, j, k, l, q !< Iterator variables
integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables
real(wp) :: coeff
real(wp) :: alpha_sum
real(wp) :: pres, dyn_pres, pres_mag, T
real(wp) :: rhoYks(1:num_species)
real(wp) :: rho_K, gamma_K, pi_inf_K, qv_K
real(wp), dimension(num_fluids) :: alpha_K, alpha_rho_K
real(wp), dimension(2) :: Re_K

i1 = gp%ip_grid(1); i2 = i1 + 1
j1 = gp%ip_grid(2); j2 = j1 + 1
Expand All @@ -861,6 +938,7 @@ contains
alpha_IP = 0._wp
pres_IP = 0._wp
vel_IP = 0._wp
pres = 0._wp

if (surface_tension) c_IP = 0._wp

Expand All @@ -887,31 +965,97 @@ contains
do j = j1, j2
$:GPU_LOOP(parallelism='[seq]')
do k = k1, k2

coeff = gp%interp_coeffs(i - i1 + 1, j - j1 + 1, k - k1 + 1)

pres_IP = pres_IP + coeff* &
q_prim_vf(E_idx)%sf(i, j, k)
if (igr) then
! For IGR, we will need to perform operations on
! the conservative variables instead
Comment on lines +970 to +972
Copy link
Contributor

@qodo-code-review qodo-code-review bot Dec 19, 2025

Choose a reason for hiding this comment

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

Suggestion: In s_interpolate_image_point, add a check to ensure q_cons_vf is present before using it inside the if (igr) block, as it is an optional argument. [possible issue, importance: 8]

Suggested change
if (igr) then
! For IGR, we will need to perform operations on
! the conservative variables instead
if (igr) then
if (.not. present(q_cons_vf)) then
stop 'q_cons_vf is required for IGR interpolation'
end if
! For IGR, we will need to perform operations on
! the conservative variables instead

alpha_sum = 0._wp
dyn_pres = 0._wp
if (num_fluids == 1) then
alpha_rho_K(1) = q_cons_vf(1)%sf(i, j, k)
alpha_K(1) = 1._wp
else
$:GPU_LOOP(parallelism='[seq]')
do l = 1, num_fluids - 1
alpha_rho_K(l) = q_cons_vf(l)%sf(i, j, k)
alpha_K(l) = q_cons_vf(E_idx + l)%sf(i, j, k)
end do
alpha_rho_K(num_fluids) = q_cons_vf(num_fluids)%sf(i, j, k)
alpha_K(num_fluids) = 1._wp - sum(alpha_K(1:num_fluids - 1))
end if
if (model_eqns /= 4) then
if (elasticity) then
! call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, &
! alpha_rho_K, Re_K, G_K, Gs_vc)
Comment on lines +989 to +990
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Dec 19, 2025 •
edited
Loading

Choose a reason for hiding this comment

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

P1: Variables rho_K, gamma_K, pi_inf_K, and qv_K will be uninitialized when both elasticity and igr are true. The elasticity case is commented out, but these variables are still used in subsequent calculations (max(rho_K, sgm_eps), division by rho_K, and s_compute_pressure call). Either uncomment and fix the elasticity call, or add a runtime check to prevent this code path.

Prompt for AI agents 
Check if this issue is valid — if so, understand the root cause and fix it. At src/simulation/m_ibm.fpp, line 982:

<comment>Variables `rho_K`, `gamma_K`, `pi_inf_K`, and `qv_K` will be uninitialized when both `elasticity` and `igr` are true. The elasticity case is commented out, but these variables are still used in subsequent calculations (`max(rho_K, sgm_eps)`, division by `rho_K`, and `s_compute_pressure` call). Either uncomment and fix the elasticity call, or add a runtime check to prevent this code path.</comment>

<file context>
@@ -887,31 +958,96 @@ contains
+                        end if
+                        if (model_eqns /= 4) then
+                            if (elasticity) then
+!                            call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, &amp;
+!                                                                            alpha_rho_K, Re_K, G_K, Gs_vc)
+                            else
</file context>
Suggested change
! call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, &
! alpha_rho_K, Re_K, G_K, Gs_vc)
call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, alpha_K, &
alpha_rho_K, Re_K, G_K, Gs)
Fix with Cubic

else
call s_convert_species_to_mixture_variables_acc(rho_K, gamma_K, pi_inf_K, qv_K, &
alpha_K, alpha_rho_K, Re_K)
end if
end if

$:GPU_LOOP(parallelism='[seq]')
do q = momxb, momxe
vel_IP(q + 1 - momxb) = vel_IP(q + 1 - momxb) + coeff* &
q_prim_vf(q)%sf(i, j, k)
end do
rho_K = max(rho_K, sgm_eps)

$:GPU_LOOP(parallelism='[seq]')
do l = contxb, contxe
alpha_rho_IP(l) = alpha_rho_IP(l) + coeff* &
q_prim_vf(l)%sf(i, j, k)
alpha_IP(l) = alpha_IP(l) + coeff* &
q_prim_vf(advxb + l - 1)%sf(i, j, k)
end do
$:GPU_LOOP(parallelism='[seq]')
do l = momxb, momxe
if (model_eqns /= 4) then
dyn_pres = dyn_pres + 5.e-1_wp*q_cons_vf(l)%sf(i, j, k) &
*q_cons_vf(l)%sf(i, j, k)/rho_K
end if
end do
pres_mag = 0._wp

call s_compute_pressure(q_cons_vf(E_idx)%sf(i, j, k), &
q_cons_vf(alf_idx)%sf(i, j, k), &
dyn_pres, pi_inf_K, gamma_K, rho_K, &
qv_K, rhoYks, pres, T, pres_mag=pres_mag)

pres_IP = pres_IP + coeff*pres

$:GPU_LOOP(parallelism='[seq]')
do q = momxb, momxe
vel_IP(q + 1 - momxb) = vel_IP(q + 1 - momxb) + coeff* &
q_cons_vf(q)%sf(i, j, k)/rho_K
end do

if (surface_tension) then
if (num_fluids == 1) then
alpha_rho_IP(1) = alpha_rho_IP(1) + coeff*q_cons_vf(contxb)%sf(i, j, k)
alpha_IP(1) = alpha_IP(1) + coeff*1._wp
else
alpha_sum = 0._wp
$:GPU_LOOP(parallelism='[seq]')
do l = 1, num_fluids - 1
alpha_rho_IP(l) = alpha_rho_IP(l) + coeff*q_cons_vf(l)%sf(i, j, k)
alpha_IP(l) = alpha_IP(l) + coeff*q_cons_vf(E_idx + l)%sf(i, j, k)
alpha_sum = alpha_sum + q_cons_vf(E_idx + l)%sf(i, j, k)
end do
alpha_rho_IP(num_fluids) = alpha_rho_IP(num_fluids) + coeff*q_cons_vf(num_fluids)%sf(i, j, k)
alpha_IP(num_fluids) = alpha_IP(num_fluids) + coeff*(1._wp - alpha_sum)
end if
else
pres_IP = pres_IP + coeff* &
q_prim_vf(E_idx)%sf(i, j, k)

$:GPU_LOOP(parallelism='[seq]')
do q = momxb, momxe
vel_IP(q + 1 - momxb) = vel_IP(q + 1 - momxb) + coeff* &
q_prim_vf(q)%sf(i, j, k)
end do

$:GPU_LOOP(parallelism='[seq]')
do l = contxb, contxe
alpha_rho_IP(l) = alpha_rho_IP(l) + coeff* &
q_prim_vf(l)%sf(i, j, k)
alpha_IP(l) = alpha_IP(l) + coeff* &
q_prim_vf(advxb + l - 1)%sf(i, j, k)
end do
end if

if (surface_tension .and. .not. igr) then
c_IP = c_IP + coeff*q_prim_vf(c_idx)%sf(i, j, k)
end if

if (bubbles_euler .and. .not. qbmm) then
if (bubbles_euler .and. .not. qbmm .and. .not. igr) then
$:GPU_LOOP(parallelism='[seq]')
do l = 1, nb
if (polytropic) then
Expand All @@ -926,7 +1070,7 @@ contains
end do
end if

if (qbmm) then
if (qbmm .and. .not. igr) then
do l = 1, nb*nmom
nmom_IP(l) = nmom_IP(l) + coeff*q_prim_vf(bubxb - 1 + l)%sf(i, j, k)
end do
Expand All @@ -940,9 +1084,7 @@ contains
end do
end do
end if

end if

end do
end do
end do
Expand Down
Loading
Loading