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

Closed
Cowsreal wants to merge 28 commits into from
Closed

Immersed boundaries integration with IGR #1095

Show file tree
Hide file tree
Changes from 2 commits
Commits
Show all changes
28 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
d914ba3
add igr forward facing step example
wilfonba Dec 29, 2025
3e92773
update IGR forward facing step example
wilfonba Jan 1, 2026
4be08cc
Minor fixes
Jan 2, 2026
27ea826
Added to skipped cases and README
Jan 2, 2026
845c9b5
Minor fixes and formatting
Jan 2, 2026
a8c9629
Fix inline issue on Frontier
Jan 3, 2026
01e6492
Fix inlining issue on frontier
Jan 3, 2026
e28a2c9
Remove q_cons_vf optional variable guard
Jan 3, 2026
a914653
Restore elasticity branch
Jan 3, 2026
a29bac8
Format
Jan 3, 2026
b31dd3a
Minor fixes and add guard for IGR + MIBM
Jan 3, 2026
5a6a5a8
Minor changes
Jan 4, 2026
bcce655
Merge branch 'master' into igribm
sbryngelson Jan 4, 2026
31fce2a
Minor changes regards to chemistry. Add guard inside case_validator f…
Cowsreal Jan 5, 2026
6c756a0
merge
wilfonba Jan 18, 2026
e111b08
Formatting
Jan 21, 2026
c89cb66
Merge branch 'master' into igribm
Cowsreal Jan 21, 2026
6cf419c
Merge branch 'master' into igribm
Cowsreal Jan 22, 2026
b68f694
Add 3D IGR bowshock case
Jan 24, 2026
755b013
Merge branch 'igribm' of github.com:Cowsreal/MFCMarkZhang into igribm
Jan 24, 2026
2a40bbf
Merge branch 'master' into igribm
Cowsreal Jan 24, 2026
849a030
Fix README format
Jan 24, 2026
0009ca6
Merge branch 'igribm' of github.com:Cowsreal/MFCMarkZhang into igribm
Jan 24, 2026
1b28f76
Merge branch 'master' into igribm
wilfonba Feb 9, 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
242 changes: 188 additions & 54 deletions src/simulation/m_ibm.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -150,6 +150,44 @@ contains

end subroutine s_populate_ib_buffers

subroutine s_update_igr(jac_sf)
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

$: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

!> 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 +197,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 +234,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='[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
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 +279,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
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 +380,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 +400,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 +888,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 +911,12 @@ contains
integer :: i, j, k, l, q !< Iterator variables
integer :: i1, i2, j1, j2, k1, k2 !< Iterator variables
real(wp) :: coeff
real(wp) :: alphaSum
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 +931,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 +958,96 @@ 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
alphaSum = 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)
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)
Copy link
Copy Markdown

@cursor cursor bot Jan 4, 2026

Choose a reason for hiding this comment

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

Uninitialized variables in IGR path for chemistry mode

The IGR code path in s_interpolate_image_point declares rhoYks and T (line 926-927) but never initializes them before passing to s_compute_pressure (line 1013-1016). In chemistry mode (compile-time flag), s_compute_pressure reads rhoYks to compute species fractions and uses T as a Newton solver initial guess. Since chemistry is not prohibited with IGR in the validator, using IGR + IBM with a chemistry-enabled build would result in incorrect pressure calculations due to uninitialized species fractions and potentially poor temperature solver convergence.

Fix in Cursor Fix in Web

Copy link
Copy Markdown
Contributor Author

@Cowsreal Cowsreal Jan 5, 2026

Choose a reason for hiding this comment

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

IGR is not supported with chemistry yet, so a guard needs to be added inside the validator. It's ok to leave rhoYks and T as is because they're never used inside s_compute_pressure.

Copy link
Copy Markdown
Contributor

@qodo-code-review qodo-code-review bot Jan 5, 2026

Choose a reason for hiding this comment

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

Code Implementation 🛠️

Implementation: Initialize rhoYks and T before calling s_compute_pressure in the IGR interpolation path to avoid passing uninitialized values in chemistry-enabled builds.

Suggested change
qv_K, rhoYks, pres, T, pres_mag=pres_mag)
! Ensure chemistry-related arguments are initialized even if chemistry is not supported with IGR yet.
rhoYks = 0._wp
T = 300._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)

See review comment here

Copy link
Copy Markdown
Contributor Author

@Cowsreal Cowsreal Jan 5, 2026

Choose a reason for hiding this comment

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

You're right, I should initialize them with values just for code safety anyways.

Copy link
Copy Markdown
Contributor

@qodo-code-review qodo-code-review bot Jan 5, 2026

Choose a reason for hiding this comment

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

Code Implementation 🛠️

Implementation: Initialize chemistry-related temporaries (rhoYks and T) in the IGR interpolation path before calling s_compute_pressure to avoid passing uninitialized values in chemistry-enabled builds.

Suggested change
qv_K, rhoYks, pres, T, pres_mag=pres_mag)
! Ensure chemistry-related temporaries are initialized even if chemistry is not supported with IGR yet.
rhoYks = 0._wp
T = 300._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)

See review comment here


Comment on lines +1020 to +1048
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai bot Feb 9, 2026

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Guard IGR pressure computation for model_eqns == 4.

When model_eqns == 4, rho_K, gamma_K, pi_inf_K, and qv_K are never initialized but are still used in s_compute_pressure. Either compute these quantities for this branch or fail fast.

🛠️ Suggested guard 
-                        if (model_eqns /= 4) then
+                        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)
                             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
+                        else
+                            call s_mpi_abort("IGR image-point interpolation requires mixture variables; model_eqns == 4 is not yet supported here")
                         end if

As per coding guidelines: "Call s_mpi_abort() for errors, never use stop or error stop".

🤖 Prompt for AI Agents 
In `@src/simulation/m_ibm.fpp` around lines 1020 - 1048, The branch allows
model_eqns==4 to reach the call to s_compute_pressure while rho_K, gamma_K,
pi_inf_K and qv_K may be uninitialized; fix by adding an explicit guard before
calling s_compute_pressure: if model_eqns == 4 either populate these variables
(e.g., call s_convert_species_to_mixture_variables_acc with the appropriate args
as done for other branches) or call s_mpi_abort("<descriptive message>") to fail
fast; update the logic around the use of
s_convert_species_to_mixture_variables_acc, the rho_K = max(...) line, and the
s_compute_pressure invocation to ensure the chosen branch always leaves rho_K,
gamma_K, pi_inf_K and qv_K initialized.

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 (num_fluids == 1) then
alpha_rho_IP(1) = alpha_rho_IP(1) + coeff*q_cons_vf(contxb)%sf(i, j, k)
Copy link
Copy Markdown
Contributor

@codeant-ai codeant-ai bot Jan 22, 2026

Choose a reason for hiding this comment

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

Suggestion: In the IGR interpolation branch for the single-fluid case, the code uses q_cons_vf(contxb) when accumulating alpha_rho_IP but earlier treated the single-fluid conservative density as q_cons_vf(1) — this index mismatch will pick the wrong field and produce incorrect interpolated mixture density; use a consistent index for the single-fluid conservative variable. [logic error]

Severity Level: Critical 🚨 
- ❌ IGR single-fluid interpolation uses wrong conservative field.
- ⚠️ Ghost-point mixture density interpolated incorrectly.
- ⚠️ Affects IGR-based immersed boundary corrections.
Suggested change
alpha_rho_IP(1) = alpha_rho_IP(1) + coeff*q_cons_vf(contxb)%sf(i, j, k)
alpha_rho_IP(1) = alpha_rho_IP(1) + coeff*q_cons_vf(1)%sf(i, j, k)
Steps of Reproduction ✅ 
1. s_interpolate_image_point contains an IGR branch starting at
src/simulation/m_ibm.fpp:989; inside that branch the single-fluid handling sets
alpha_rho_IP using q_cons_vf(contxb) at src/simulation/m_ibm.fpp:1044.

2. The only caller that supplies q_cons_vf is the IGR caller in s_ibm_correct_state at
src/simulation/m_ibm.fpp:307, so this path is exercised when igr = .true. and num_fluids
== 1.

3. Run an IGR-enabled single-fluid case (igr true, num_fluids == 1). During image-point
interpolation the code will read q_cons_vf(contxb)%sf(i,j,k)
(src/simulation/m_ibm.fpp:1044) instead of q_cons_vf(1)%sf(...). If contxb ≠ 1 this reads
the wrong conserved field.

4. Observe incorrect interpolated density/mixture variables and resulting incorrect
pressure/velocity at ghost points for IGR single-fluid runs; this will affect IB state
corrections and nearby solution fields.
Prompt for AI Agent 🤖 
This is a comment left during a code review.

**Path:** src/simulation/m_ibm.fpp
**Line:** 1044:1044
**Comment:**
	*Logic Error: In the IGR interpolation branch for the single-fluid case, the code uses `q_cons_vf(contxb)` when accumulating `alpha_rho_IP` but earlier treated the single-fluid conservative density as `q_cons_vf(1)` — this index mismatch will pick the wrong field and produce incorrect interpolated mixture density; use a consistent index for the single-fluid conservative variable.

Validate the correctness of the flagged issue. If correct, How can I resolve this? If you propose a fix, implement it and please make it concise.

Copy link
Copy Markdown
Contributor Author

@Cowsreal Cowsreal Jan 24, 2026

Choose a reason for hiding this comment

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

It's ok, if num_fluids is 1, then contxb will be 1.

Copy link
Copy Markdown
Contributor

@qodo-code-review qodo-code-review bot Jan 24, 2026

Choose a reason for hiding this comment

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

Question:

It's ok, if num_fluids is 1, then contxb will be 1.

Answer:

Based on the thread, the “logic error” depends on an invariant that the AI reviewer does not know from the diff alone: whether contxb is guaranteed to be 1 whenever num_fluids == 1. If that invariant is true, then q_cons_vf(contxb) and q_cons_vf(1) are equivalent and the suggestion is a false positive. If it’s not guaranteed (or could change in future refactors), then the suggestion is valid.

To reduce this type of false-positive “index mismatch” suggestion in similar situations, add project guidance via pr_code_suggestions.extra_instructions to require the model to avoid reporting “mismatches” when they may be enforced by invariants, and to ask for confirmation (or point to where the invariant is defined) before proposing a change.

Recommended configuration (repo .pr_agent.toml or wiki .pr_agent.toml):

[pr_code_suggestions]
extra_instructions = """\
- When flagging index/array-field mismatches (e.g., q_cons_vf(contxb) vs q_cons_vf(1)), do NOT assume it's a bug unless the diff shows the invariant is violated.
- If correctness depends on an external invariant (e.g., contxb==1 when num_fluids==1), explicitly state the assumption and ask for confirmation instead of proposing a code change.
- Prefer suggestions that are provably wrong from the local context; avoid speculative “logic error” claims.
"""

If you want an even stronger filter toward only clear bugs:

[pr_code_suggestions]
focus_only_on_problems = true

If you share where contxb is set (or a comment/doc that states the invariant), you can make the instruction more specific (e.g., “In this codebase, contxb is always 1 for single-fluid; don’t flag it”).

Relevant Sources:

alpha_IP(1) = alpha_IP(1) + coeff*1._wp
else
$: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)
alphaSum = alphaSum + 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 - alphaSum)
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

if (surface_tension) then
$: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
Comment on lines +1090 to 1092
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai bot Feb 9, 2026

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Surface-tension color function is skipped in IGR.

With igr enabled, c_IP never gets interpolated, so ghost cells will receive zero color function. Either keep primitive interpolation or interpolate from q_cons_vf(c_idx) when IGR is on.

🛠️ Suggested fix 
-                    if (surface_tension .and. .not. igr) then
-                        c_IP = c_IP + coeff*q_prim_vf(c_idx)%sf(i, j, k)
-                    end if
+                    if (surface_tension) then
+                        if (igr) then
+                            c_IP = c_IP + coeff*q_cons_vf(c_idx)%sf(i, j, k)
+                        else
+                            c_IP = c_IP + coeff*q_prim_vf(c_idx)%sf(i, j, k)
+                        end if
+                    end if
🤖 Prompt for AI Agents 
In `@src/simulation/m_ibm.fpp` around lines 1090 - 1092, The current block skips
adding the surface-tension color contribution to c_IP when igr is true, leaving
ghost cells with zero; change the logic so that if surface_tension is true you
always add coeff * <color function>, using q_prim_vf(c_idx)%sf(i,j,k) when igr
is false and using q_cons_vf(c_idx)%sf(i,j,k) (or otherwise perform the
primitive interpolation from conserved fields) when igr is true; update the
conditional around c_IP (references: surface_tension, igr, c_IP, coeff,
q_prim_vf, q_cons_vf, c_idx, sf) to ensure c_IP is interpolated in both modes.


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 +1062,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 +1076,7 @@ contains
end do
end do
end if

end if

end do
end do
end do
Expand Down
29 changes: 28 additions & 1 deletion src/simulation/m_igr.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@ module m_igr

use m_boundary_common

use m_ibm

implicit none

private; public :: s_initialize_igr_module, &
Expand Down Expand Up @@ -368,7 +370,32 @@ contains
$:END_GPU_PARALLEL_LOOP()
end if
end do

if (ib) then
$:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
do l = 0, p
do k = 0, n
do j = 0, m
if (ib_markers%sf(j, k, l) /= 0) then
jac(j, k, l) = 0._wp
end if
end do
end do
end do
$:END_GPU_PARALLEL_LOOP()
call s_update_igr(jac_sf)
! If using Jacobi Iter, we need to update jac_old again
if (igr_iter_solver == 1) then
$:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
do l = idwbuff(3)%beg, idwbuff(3)%end
do k = idwbuff(2)%beg, idwbuff(2)%end
do j = idwbuff(1)%beg, idwbuff(1)%end
jac_old(j, k, l) = jac(j, k, l)
end do
end do
end do
$:END_GPU_PARALLEL_LOOP()
end if
end if
Comment on lines +366 to +391
Copy link
Copy Markdown
Contributor

@qodo-code-review qodo-code-review bot Jan 7, 2026

Choose a reason for hiding this comment

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

Suggestion: In s_igr_iterative_solve, call s_interpolate_sigma_igr with jac_old as the source for interpolation. This ensures ghost-point values are calculated from the field before it is zeroed-out within the immersed boundary, preventing corrupted stencil values. [possible issue, importance: 9]

Suggested change
if (ib) then
$:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
do l = 0, p
do k = 0, n
do j = 0, m
if (ib_markers%sf(j, k, l) /= 0) then
jac(j, k, l) = 0._wp
end if
end do
end do
end do
$:END_GPU_PARALLEL_LOOP()
call s_interpolate_sigma_igr(jac)
! If using Jacobi Iter, we need to update jac_old again
if (igr_iter_solver == 1) then
$:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
do l = idwbuff(3)%beg, idwbuff(3)%end
do k = idwbuff(2)%beg, idwbuff(2)%end
do j = idwbuff(1)%beg, idwbuff(1)%end
jac_old(j, k, l) = jac(j, k, l)
end do
end do
end do
$:END_GPU_PARALLEL_LOOP()
end if
end if
if (ib) then
$:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
do l = 0, p
do k = 0, n
do j = 0, m
if (ib_markers%sf(j, k, l) /= 0) then
jac(j, k, l) = 0._wp
end if
end do
end do
end do
$:END_GPU_PARALLEL_LOOP()
call s_interpolate_sigma_igr(jac, jac_old)
if (igr_iter_solver == 1) then
$:GPU_PARALLEL_LOOP(private='[j,k,l]', collapse=3)
do l = idwbuff(3)%beg, idwbuff(3)%end
do k = idwbuff(2)%beg, idwbuff(2)%end
do j = idwbuff(1)%beg, idwbuff(1)%end
jac_old(j, k, l) = jac(j, k, l)
end do
end do
end do
$:END_GPU_PARALLEL_LOOP()
end if
end if

Copy link
Copy Markdown
Contributor Author

@Cowsreal Cowsreal Jan 8, 2026

Choose a reason for hiding this comment

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

Good point. Although this won't work for Gauss Siedel because jac_old is not allocated in that case to space space. However, now that you bring this up, the rest of ibm module also suffers from this, so perhaps we can try using atomic writes.

Copy link
Copy Markdown
Contributor

@qodo-code-review qodo-code-review bot Jan 8, 2026

Choose a reason for hiding this comment

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

Question:

Good point. Although this won't work for Gauss Siedel because jac_old is not allocated in that case to space space. However, now that you bring this up, the rest of ibm module also suffers from this, so perhaps we can try using atomic writes.

Answer:

Add project-specific guidance so /improve avoids suggestions that rely on arrays/buffers that may not exist in some solver modes (e.g., Gauss-Seidel), and instead proposes safe alternatives (conditional allocation, in-place compute ordering, atomics, or per-solver branches).

[pr_code_suggestions]
extra_instructions = """\
For solver/algorithm variants (e.g., Gauss-Seidel vs others), do NOT suggest using a work array unless you can confirm it is allocated/available in that mode.
- If a suggestion would require an unallocated array (like jac_old in Gauss-Seidel), propose one of:
  (a) conditional allocation with clear guards,
  (b) reordering/in-place computation that preserves required values,
  (c) atomic writes / race-safe parallel updates (GPU/OMP) when relevant,
  (d) per-solver branching (igr_iter_solver dependent) with minimal code duplication.
- When proposing changes to GPU-parallel loops, explicitly call out any race conditions and only recommend atomics when correctness requires it.
- Prefer solutions that are correct for all solver modes; if not possible, state the mode limitation explicitly.
"""

If you want even fewer “speculative” refactors in performance-sensitive GPU code, you can also keep suggestions focused on correctness issues:

[pr_code_suggestions]
focus_only_on_problems = true

Relevant Sources:

end subroutine s_igr_iterative_solve

subroutine s_igr_sigma_x(q_cons_vf, rhs_vf)
Expand Down
Loading