Skip to content

Conversation

@Cowsreal
Copy link
Contributor

@Cowsreal Cowsreal commented Dec 19, 2025

User description

User description

User description

Description

Added immersed boundaries compatibility with IGR. I deleted the relevant prohibition inside the Python case checker.

Type of change

Please delete options that are not relevant.

  • New feature (non-breaking change which adds functionality)

Scope

  • This PR comprises a set of related changes with a common goal

If you cannot check the above box, please split your PR into multiple PRs that each have a common goal.

How Has This Been Tested?

Please describe the tests that you ran to verify your changes.
Provide instructions so we can reproduce.
Please also list any relevant details for your test configuration

  • Eye norm by comparing WENO5 results to IGR5 results on 500x500 grids. (Airfoil, square, circle)
  • For the airfoil example, below is are two plots showing the pressure profile on the contour of the airfoil over time. (Demonstrate the similarity between the pressure profiles)
  • For the airfoil example, below is a plot showing the average of the difference of pressure( + rho) at each point on the contour over time. (difference of pressure corresponding to IGR5 results - WENO5 results)

Plots of mean difference between IGR5 and WENO5 on airfoil contour

Figure: This shows the average of pressure differences over the airfoil contour for all time. (i.e., (mean(p1(s) - p2(s)))(t))
igr5weno5comp1

Figure: This shows the average of density differences over the airfoil contour for all time. (i.e., (mean(rho1(s) - rho2(s)))(t))
igr5weno5comprho

Movies of pressure profile around airfoil contour

IGR5 results

presContour.mp4

WENO5 results

presContour.mp4

Movies of pressure

Airfoil, IGR5

airfoil_500_500.mp4

Airfoil, WENO5

airfoil_500_500.mp4

Square and Circle have viscosity with Re = 25000

Circle, IGR5

pres.mp4

Circle, WENO5

pres.mp4

Square, IGR5

pres.mp4

Square, WENO5

pres.mp4

Test Configuration:

  • What computers and compilers did you use to test this: Phoenix(nvhpc 24.5) + Delta(nvhpc 24.1) + Frontier(Cray 19.0)

Checklist

  • I have added comments for the new code
  • I added Doxygen docstrings to the new code
  • I ran ./mfc.sh format before committing my code
  • New and existing tests pass locally with my changes, including with GPU capability enabled (both NVIDIA hardware with NVHPC compilers and AMD hardware with CRAY compilers) and disabled
  • This PR does not introduce any repeated code (it follows the DRY principle)
  • I cannot think of a way to condense this code and reduce any introduced additional line count

If your code changes any code source files (anything in src/simulation)

To make sure the code is performing as expected on GPU devices, I have:

  • Checked that the code compiles using NVHPC compilers
  • Checked that the code compiles using CRAY compilers
  • Ran the code on either V100, A100, or H100 GPUs and ensured the new feature performed as expected (the GPU results match the CPU results)
  • Ran the code on MI200+ GPUs and ensure the new features performed as expected (the GPU results match the CPU results)

PR Type

Enhancement


Description

  • Integrate immersed boundary method with IGR solver

  • Add s_update_igr subroutine for ghost point interpolation

  • Modify state correction to handle IGR-specific variable updates

  • Remove IBM-IGR incompatibility restriction from validator


Diagram Walkthrough

flowchart LR
  IBM["IBM Module"]
  IGR["IGR Solver"]
  GhostPoints["Ghost Points"]
  ImagePoints["Image Points"]
  StateCorrection["State Correction"]
  
  IBM -- "s_update_igr" --> GhostPoints
  GhostPoints -- "interpolate from" --> ImagePoints
  ImagePoints -- "conservative vars" --> StateCorrection
  IGR -- "calls" --> StateCorrection
  StateCorrection -- "updates" --> GhostPoints
Loading

File Walkthrough

Relevant files
Enhancement
m_ibm.fpp
Integrate IBM with IGR solver capabilities                             

src/simulation/m_ibm.fpp

  • Added new s_update_igr subroutine to interpolate sigma at ghost points
    using image point data
  • Modified s_ibm_correct_state to conditionally skip interior pressure
    initialization when IGR is enabled
  • Updated ghost point state interpolation to use conservative variables
    for IGR instead of primitive variables
  • Added IGR-specific handling for continuity and volume fraction
    variables at ghost points
  • Modified pressure setting logic to skip moving IBM pressure
    calculations when IGR is active
  • Updated s_interpolate_image_point signature to accept optional
    conservative variables parameter
  • Implemented IGR interpolation path that computes pressure from
    conservative variables and mixture properties
+188/-54
m_igr.fpp
Integrate IBM updates into IGR solver loop                             

src/simulation/m_igr.fpp

  • Added use m_ibm module import for IBM integration
  • Integrated IBM ghost point updates into IGR iterative solver loop
  • Added zero-out of Jacobian at IBM interior points before calling
    s_update_igr
  • Added conditional update of jac_old for Jacobi iteration method when
    IBM is active
+28/-1   
Bug fix
case_validator.py
Remove IGR-IBM incompatibility restriction                             

toolchain/mfc/case_validator.py

  • Removed prohibition that prevented simultaneous use of IGR and IBM
  • Deleted lines that enforced incompatibility check between ib and igr
    options
+0/-3     


CodeAnt-AI Description

Support immersed boundary method (IBM) when running the IGR solver

What Changed

  • IGR now accepts and processes cells marked as immersed boundaries: ghost points inside IBM are set to zero and then populated from image-point interpolation so they receive consistent field values for IGR.
  • Image-point interpolation for IGR uses conservative variables to compute pressure, velocity and phase fractions, producing interpolated pressures and velocities from the conserved state rather than from primitive-only fields.
  • IBM interior-pressure and mixture-variable handling is skipped or adjusted when IGR is active to avoid overwriting IGR-computed values; moving-IBM pressure handling remains disabled for IGR (temporary).
  • The runtime case validator no longer blocks running IGR with the immersed boundary method, allowing previously-prohibited IB+IGR configurations.

Impact

✅ Run IGR on cases with immersed boundaries
✅ Correct pressure and velocity at immersed ghost points
✅ Fewer pre-run configuration blocks for IB+IGR

💡 Usage Guide

Checking Your Pull Request

Every time you make a pull request, our system automatically looks through it. We check for security issues, mistakes in how you're setting up your infrastructure, and common code problems. We do this to make sure your changes are solid and won't cause any trouble later.

Talking to CodeAnt AI

Got a question or need a hand with something in your pull request? You can easily get in touch with CodeAnt AI right here. Just type the following in a comment on your pull request, and replace "Your question here" with whatever you want to ask:

@codeant-ai ask: Your question here

This lets you have a chat with CodeAnt AI about your pull request, making it easier to understand and improve your code.

Example

@codeant-ai ask: Can you suggest a safer alternative to storing this secret?

Preserve Org Learnings with CodeAnt

You can record team preferences so CodeAnt AI applies them in future reviews. Reply directly to the specific CodeAnt AI suggestion (in the same thread) and replace "Your feedback here" with your input:

@codeant-ai: Your feedback here

This helps CodeAnt AI learn and adapt to your team's coding style and standards.

Example

@codeant-ai: Do not flag unused imports.

Retrigger review

Ask CodeAnt AI to review the PR again, by typing:

@codeant-ai: review

Check Your Repository Health

To analyze the health of your code repository, visit our dashboard at https://app.codeant.ai. This tool helps you identify potential issues and areas for improvement in your codebase, ensuring your repository maintains high standards of code health.

Summary by CodeRabbit

  • New Features
    • IGR now supports the immersed boundary method and uses conservative-variable image-point interpolation for assembling mixture, pressure and velocity fields.
  • Bug Fixes
    • Improved synchronization of iterative solver state for immersed-boundary markers and safer handling across single- and multi-fluid cases (non-negative densities and robust mixture updates).
  • Chores
    • Removed prior restriction preventing IGR from being used with the immersed boundary method.

✏️ Tip: You can customize this high-level summary in your review settings.


PR Type

Enhancement


Description

  • Integrate immersed boundaries (IBM) with IGR method

  • Add sigma interpolation at ghost points via s_update_igr

  • Conditionally handle IBM state updates based on IGR flag

  • Remove IGR-IBM incompatibility restriction from validator


Diagram Walkthrough

flowchart LR
  IBM["Immersed Boundary<br/>Method"]
  IGR["IGR Iterative<br/>Solver"]
  GhostPts["Ghost Points"]
  SigmaInterp["Sigma Interpolation<br/>s_update_igr"]
  StateUpdate["State Update<br/>s_ibm_correct_state"]
  
  IBM --> GhostPts
  IGR --> GhostPts
  GhostPts --> SigmaInterp
  GhostPts --> StateUpdate
  SigmaInterp --> IGR
  StateUpdate --> IBM
Loading

File Walkthrough

Relevant files
Enhancement
m_ibm.fpp
IBM-IGR integration with sigma interpolation                         

src/simulation/m_ibm.fpp

  • Added s_update_igr subroutine to interpolate sigma at ghost points
    using image point coefficients
  • Modified s_ibm_correct_state to conditionally skip interior pressure
    initialization when IGR is enabled
  • Updated s_interpolate_image_point to accept optional conservative
    variables and handle IGR-specific interpolation logic
  • Added IGR-specific state variable assignments for density and volume
    fractions at ghost points
  • Wrapped pressure setting and continuity variable updates with IGR
    conditionals
  • Fixed documentation comment for q_cons_vf parameter from "Primitive
    Variables" to "Conservative Variables"
+196/-54
m_igr.fpp
IGR solver integration with IBM ghost points                         

src/simulation/m_igr.fpp

  • Added import of IBM module (m_ibm) with ib, ib_markers, and
    s_update_igr
  • Integrated IBM ghost point sigma updates into IGR iterative solver
    loop
  • Zero out Jacobian at IBM interior points before calling s_update_igr
  • Update jac_old after sigma interpolation when using Jacobi iteration
    solver
+28/-1   
Bug fix
case_validator.py
Remove IGR-IBM incompatibility restriction                             

toolchain/mfc/case_validator.py

  • Removed prohibition that prevented IGR from being used with immersed
    boundary method
  • Deleted ib variable declaration and associated validation check
+0/-3     

@codeant-ai
Copy link

codeant-ai bot commented Dec 19, 2025

CodeAnt AI is reviewing your PR.


Thanks for using CodeAnt! 🎉

We're free for open-source projects. if you're enjoying it, help us grow by sharing.

Share on X ·
Reddit ·
LinkedIn

@coderabbitai
Copy link
Contributor

coderabbitai bot commented Dec 19, 2025

Note

Other AI code review bot(s) detected

CodeRabbit has detected other AI code review bot(s) in this pull request and will avoid duplicating their findings in the review comments. This may lead to a less comprehensive review.

Walkthrough

Adds IBM ↔ IGR integration: new subroutine s_interpolate_sigma_igr(jac) to interpolate sigma at image/ghost points using conservative variables when IGR is active, updates s_interpolate_image_point to accept q_cons_vf, calls the new subroutine from the IGR solver post-iteration, and removes the validator ban on IGR+IBM. (45 words)

Changes

Cohort / File(s) Summary
IBM module updates
src/simulation/m_ibm.fpp
Adds s_interpolate_sigma_igr(jac) to interpolate sigma at ghost/image points (GPU-parallel). Extends s_interpolate_image_point(...) signature to accept optional q_cons_vf. Introduces conservative-variable (q_cons_vf) code paths gated by igr and dual handling for single-/multi-fluid cases.
IGR solver updates
src/simulation/m_igr.fpp
Adds use m_ibm, only: ib, ib_markers, s_interpolate_sigma_igr. After the iterative solve in s_igr_iterative_solve, zeroes jac at IBM-marked cells, calls s_interpolate_sigma_igr(jac), and copies jac→jac_old for Jacobi solver when appropriate.
Validator update
toolchain/mfc/case_validator.py
Removes the prior prohibition in check_igr_simulation that disallowed using IGR with IBM, permitting IGR+IBM combinations.

Sequence Diagram(s)

sequenceDiagram
    participant Solver as IGR Solver
    participant Jac as Jacobian (jac / jac_old)
    participant IBM as m_ibm (s_interpolate_sigma_igr / s_interpolate_image_point)
    participant GPU as GPU-parallel loops

    Note over Solver: main IGR iterative solve
    Solver->>IBM: call s_interpolate_image_point(..., optional q_cons_vf)
    IBM-->>Solver: interpolated image-point data

    alt ib (IBM) active — post-iteration
        Solver->>Jac: zero jac at IB markers (parallel)
        Solver->>IBM: call s_interpolate_sigma_igr(jac)
        IBM->>GPU: perform GPU-parallel interpolation & assemble q_cons_vf / sigma
        GPU->>Jac: update jac (ghost/image points)
        alt Jacobi solver
            Solver->>Jac: copy jac → jac_old (parallel)
        end
    end
Loading

Estimated code review effort

🎯 3 (Moderate) | ⏱️ ~25 minutes

  • Review changes to s_interpolate_image_point signature and all call sites for matching argument passing.
  • Inspect GPU-parallel loops in s_interpolate_sigma_igr for indexing, race conditions, and device synchronization.
  • Verify jac zeroing and jac→jac_old copy in s_igr_iterative_solve for consistency with solver iterations and Jacobi path.
  • Check conservative vs primitive variable branching (igr guards) and multi-fluid handling (rho_K, alpha_K, alpha_rho_K, sgm_eps).

Suggested labels

Review effort 3/5

Suggested reviewers

  • sbryngelson

Poem

A rabbit nudges ghostly jac, 🐇
Image points hum, then trace the track,
Conservative whispers stitch the sea,
IBM and IGR hop in harmony,
Small code hops — a tidy snack.

Pre-merge checks and finishing touches

✅ Passed checks (3 passed)
Check name Status Explanation
Title check ✅ Passed The title clearly and concisely summarizes the main change: adding integration between immersed boundaries and the IGR solver, which is the primary objective of this PR.
Description check ✅ Passed The description includes a summary of changes, type of change, scope, test configuration and results with supporting plots/videos, and a comprehensive checklist. However, most checklist items remain unchecked, indicating testing claims may not be fully documented per template requirements.
Docstring Coverage ✅ Passed No functions found in the changed files to evaluate docstring coverage. Skipping docstring coverage check.
✨ Finishing touches
🧪 Generate unit tests (beta)
  • Create PR with unit tests
  • Post copyable unit tests in a comment

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands and usage tips.

@Cowsreal Cowsreal changed the title Igribm Immersed boundaries integration with IGR Dec 19, 2025
@qodo-code-review qodo-code-review bot changed the title Immersed boundaries integration with IGR Igribm Dec 19, 2025
@codeant-ai codeant-ai bot added the size:L This PR changes 100-499 lines, ignoring generated files label Dec 19, 2025
@codeant-ai codeant-ai bot changed the title Igribm Immersed boundaries integration with IGR Dec 19, 2025
@qodo-code-review
Copy link
Contributor

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪
🧪 No relevant tests
🔒 No security concerns identified
⚡ Recommended focus areas for review

Possible Issue

In the IGR interpolation path, the computation of the last volume fraction uses an unweighted running sum of volume fractions from cell centers (not multiplied by the interpolation coefficient), which can produce incorrect interpolated alpha_IP and potentially violate bounds/conservation. This should be validated carefully because it affects mixture properties and pressure reconstruction at ghost points.

                    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

                        rho_K = max(rho_K, sgm_eps)

                        $: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 (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
                            $: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
Possible Issue

s_interpolate_image_point now relies on the optional q_cons_vf when igr is enabled, but there is no local guard ensuring the optional argument is present before dereferencing it. If any call path sets igr but forgets to pass q_cons_vf, this will be undefined behavior. Consider enforcing presence (e.g., if (igr .and. .not. present(q_cons_vf)) error stop) or refactoring the signature to avoid optional in the IGR-required path.

    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, 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

        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
        real(wp), optional, dimension(:), intent(INOUT) :: presb_IP, massv_IP

        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
        k1 = gp%ip_grid(3); k2 = k1 + 1

        if (p == 0) then
            k1 = 0
            k2 = 0
        end if

        alpha_rho_IP = 0._wp
        alpha_IP = 0._wp
        pres_IP = 0._wp
        vel_IP = 0._wp
        pres = 0._wp

        if (surface_tension) c_IP = 0._wp

        if (bubbles_euler) then
            r_IP = 0._wp
            v_IP = 0._wp
            if (.not. polytropic) then
                mv_IP = 0._wp
                pb_IP = 0._wp
            end if
        end if

        if (qbmm) then
            nmom_IP = 0._wp
            if (.not. polytropic) then
                presb_IP = 0._wp
                massv_IP = 0._wp
            end if
        end if

        $:GPU_LOOP(parallelism='[seq]')
        do i = i1, i2
            $:GPU_LOOP(parallelism='[seq]')
            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)

                    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

                        rho_K = max(rho_K, sgm_eps)

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

The ghost-point pressure handling is explicitly skipped for IGR with a comment indicating it is temporary and “NEED TO FIX FOR MOVING IGR”. This creates a behavioral gap for moving immersed boundaries when IGR is enabled (pressure may remain unset/incorrect), and the case validator now allows IGR+IBM. This should be either implemented or explicitly prohibited/validated for moving IBM cases to avoid silent wrong physics.

! set the pressure
! !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

Copy link
Contributor

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
...

Comment on lines 1014 to 1026
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
$: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
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: Reset alphaSum to zero for each interpolation stencil point to fix an incorrect volume fraction calculation. [possible issue, importance: 9]

Suggested change
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
$: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
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
alphaSum = 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)
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

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Good catch.

Comment on lines 373 to 398
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
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Update jac_old for the Gauss-Seidel solver (igr_iter_solver == 2) in addition to the Jacobi solver to ensure correct convergence checks. [possible issue, importance: 8]

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_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
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.
! This is also needed for Gauss-Seidel convergence check.
if (igr_iter_solver == 1 .or. igr_iter_solver == 2) 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
Contributor Author

Choose a reason for hiding this comment

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

Not relevant.

Comment on lines 889 to 1004
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
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
real(wp), optional, dimension(:), intent(INOUT) :: presb_IP, massv_IP
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
k1 = gp%ip_grid(3); k2 = k1 + 1
if (p == 0) then
k1 = 0
k2 = 0
end if
alpha_rho_IP = 0._wp
alpha_IP = 0._wp
pres_IP = 0._wp
vel_IP = 0._wp
pres = 0._wp
if (surface_tension) c_IP = 0._wp
if (bubbles_euler) then
r_IP = 0._wp
v_IP = 0._wp
if (.not. polytropic) then
mv_IP = 0._wp
pb_IP = 0._wp
end if
end if
if (qbmm) then
nmom_IP = 0._wp
if (.not. polytropic) then
presb_IP = 0._wp
massv_IP = 0._wp
end if
end if
$:GPU_LOOP(parallelism='[seq]')
do i = i1, i2
$:GPU_LOOP(parallelism='[seq]')
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
Contributor

Choose a reason for hiding this comment

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

Suggestion: Add a check to ensure the optional q_cons_vf argument is provided to s_interpolate_image_point when igr is enabled, and stop with an error if it is missing. [possible issue, importance: 7]

Suggested change
subroutine s_interpolate_image_point(..., presb_IP, massv_IP, q_cons_vf)
...
if (igr .and. .not. present(q_cons_vf)) then
error stop 'q_cons_vf must be provided when igr is enabled'
end if
if (igr) then
alphaSum = 0._wp
...
alpha_rho_K(1) = q_cons_vf(1)%sf(i, j, k)
...
else
...
end if

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Not relevant.

@codeant-ai
Copy link

codeant-ai bot commented Dec 19, 2025

Nitpicks 🔍

🔒 No security issues identified
⚡ Recommended areas for review

  • GPU data mapping
    The new GPU parallel kernel in s_update_igr uses jac_sf but the parallel-loop pragma does not include jac_sf (or its internal %sf) in its data clauses. This can lead to using host memory on device or stale data on the accelerator. Verify device data movement (copyin/copy/update) for jac_sf before the kernel runs.

  • Bounds Mismatch
    The new loops that zero out jac iterate j=0..m, k=0..n, l=0..p while the rest of the code initializes/updates jac and jac_old over idwbuff(...)%beg:...%end. If idwbuff does not equal 0..m / 0..n / 0..p this can cause out-of-bounds accesses or leave ghost/halo regions inconsistent. Verify indexing and use consistent buffer ranges.

  • Optional argument safety
    s_interpolate_image_point was extended with an optional q_cons_vf argument and the routine reads from it when igr is set. There is no runtime guard ensuring q_cons_vf is present when igr is true — this can lead to undefined behavior if the subroutine is ever called with igr true but without passing q_cons_vf.

  • Divergent control paths for IGR vs non-IGR
    Several assignment branches were added to treat igr differently (e.g., filling q_cons_vf vs q_prim_vf). Ensure all variables used later are set in both branches (or are guarded) to avoid reading uninitialized or inconsistent state.

  • Changed initialization behavior
    The initialization of moving-IBM interior pressure values (previously always set) is now only performed when .not. igr. Confirm this behavioral change is intended: it affects how interior q_prim_vf(E_idx) is set and might change results for mixed/transition cases.


use m_boundary_common

use m_ibm
Copy link

@codeant-ai codeant-ai 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: The PR adds a broad use m_ibm statement which imports the entire module namespace; importing all symbols increases risk of name collisions and makes it unclear which symbols are used by this file. Restrict the import to only the symbols actually required (for example ib, ib_markers, s_update_igr) to avoid accidental symbol conflicts and improve code clarity. [possible bug]

Severity Level: Critical 🚨

Suggested change
use m_ibm
use m_ibm, only: ib, ib_markers, s_update_igr
Why it matters? ⭐

Restricting module imports with "only:" is a good Fortran practice and improves clarity. In this file the newly added code uses ib, ib_markers and s_update_igr from m_ibm; switching to "use m_ibm, only: ib, ib_markers, s_update_igr" reduces namespace pollution and the risk of accidental name clashes. This is stylistic but sensible and safe, provided you list every symbol from m_ibm that this module actually uses.

Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_igr.fpp
**Line:** 18:18
**Comment:**
	*Possible Bug: The PR adds a broad `use m_ibm` statement which imports the entire module namespace; importing all symbols increases risk of name collisions and makes it unclear which symbols are used by this file. Restrict the import to only the symbols actually required (for example `ib`, `ib_markers`, `s_update_igr`) to avoid accidental symbol conflicts and improve code clarity.

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
Contributor Author

Choose a reason for hiding this comment

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

Good suggestion.

@codeant-ai
Copy link

codeant-ai bot commented Dec 19, 2025

CodeAnt AI finished reviewing your PR.

Copy link
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 3

Caution

Some comments are outside the diff and can’t be posted inline due to platform limitations.

⚠️ Outside diff range comments (1)
src/simulation/m_ibm.fpp (1)

891-1083: Add present(q_cons_vf) guard to IGR branch in s_interpolate_image_point

The if (igr) branch at line 963 dereferences q_cons_vf without checking if it is present. While the call site in s_ibm_correct_state only passes q_cons_vf= in the explicit else if (igr) branch, the mutually-exclusive if-else-if structure means that when both bubbles_euler=.true. and igr=.true. are set, the call enters the bubbles_euler branch without passing q_cons_vf, but the subroutine still executes the if (igr) block internally with an absent optional argument. The Fortran standard prohibits accessing an omitted optional argument, creating undefined behavior on lines 972–981 and 995–996.

Change line 963 from:

if (igr) then

to:

if (igr .and. present(q_cons_vf)) then

Optionally add @:ASSERT(.not. igr .or. present(q_cons_vf), 'IGR requires q_cons_vf in s_interpolate_image_point') near the top of the subroutine for early detection of unsupported configurations.

🧹 Nitpick comments (2)
src/simulation/m_igr.fpp (1)

18-18: New dependency on m_ibm looks reasonable but widens coupling

Bringing m_ibm into m_igr to access ib_markers and s_update_igr is consistent with the new IBM–IGR workflow and doesn’t introduce a circular dependency (m_ibm doesn’t use m_igr). Just be aware this further couples the IGR solver to IBM internals; if you intend s_update_igr to stay IBM-internal, consider exposing only what’s needed via an only: list on the use statement in a follow‑up.

src/simulation/m_ibm.fpp (1)

315-328: Pressure setting skipped for IGR is fine but note the TODO for moving IGR

The pressure assignment block:

! !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
    ...
  end if
end if

correctly avoids reusing primitive-variable pressure logic for IGR cases, which now build ghost states from conservative variables. The in-code TODO acknowledges that moving IBM + IGR will eventually need a dedicated treatment; for static or currently tested configurations the guard is appropriate.

I’d just suggest turning this into a more explicit FIXME (or opening a tracked issue) so the moving‑IGR pressure model doesn’t get forgotten.

📜 Review details

Configuration used: defaults

Review profile: CHILL

Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between 2902000 and 1340497.

📒 Files selected for processing (3)
  • src/simulation/m_ibm.fpp (11 hunks)
  • src/simulation/m_igr.fpp (2 hunks)
  • toolchain/mfc/case_validator.py (0 hunks)
💤 Files with no reviewable changes (1)
  • toolchain/mfc/case_validator.py
🧰 Additional context used
📓 Path-based instructions (3)
**/*.{fpp,f90}

📄 CodeRabbit inference engine (.github/copilot-instructions.md)

**/*.{fpp,f90}: Use 2-space indentation; continuation lines align beneath &
Use lower-case keywords and intrinsics (do, end subroutine, etc.)
Name modules with m_ pattern (e.g., m_transport)
Name public subroutines with s_ pattern (e.g., s_compute_flux)
Name public functions with f
_ pattern
Keep subroutine size ≤ 500 lines, helper subroutines ≤ 150 lines, functions ≤ 100 lines, files ≤ 1000 lines
Limit routine arguments to ≤ 6; use derived-type params struct if more are needed
Forbid goto statements (except in legacy code), COMMON blocks, and save globals
Every argument must have explicit intent; use dimension/allocatable/pointer as appropriate
Call s_mpi_abort() for errors, never use stop or error stop

**/*.{fpp,f90}: Indent 2 spaces; continuation lines align under &
Use lower-case keywords and intrinsics (do, end subroutine, etc.)
Name modules with m_<feature> prefix (e.g., m_transport)
Name public subroutines as s_<verb>_<noun> (e.g., s_compute_flux) and functions as f_<verb>_<noun>
Keep private helpers in the module; avoid nested procedures
Enforce size limits: subroutine ≤ 500 lines, helper ≤ 150, function ≤ 100, module/file ≤ 1000
Limit subroutines to ≤ 6 arguments; otherwise pass a derived-type 'params' struct
Avoid goto statements (except unavoidable legacy); avoid global state (COMMON, save)
Every variable must have intent(in|out|inout) specification and appropriate dimension / allocatable / pointer
Use s_mpi_abort(<msg>) for error termination instead of stop
Use !> style documentation for header comments; follow Doxygen Fortran format with !! @param and !! @return tags
Use implicit none statement in all modules
Use private declaration followed by explicit public exports in modules
Use derived types with pointers for encapsulation (e.g., pointer, dimension(:,:,:) => null())
Use pure and elemental attributes for side-effect-free functions; combine them for array ...

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
src/simulation/**/*.{fpp,f90}

📄 CodeRabbit inference engine (.github/copilot-instructions.md)

src/simulation/**/*.{fpp,f90}: Wrap tight GPU loops with !$acc parallel loop gang vector default(present) reduction(...); add collapse(n) when safe; declare loop-local variables with private(...)
Allocate large GPU arrays with managed memory or move them into persistent !$acc enter data regions at start-up
Avoid stop/error stop inside GPU device code
Ensure GPU code compiles with Cray ftn, NVIDIA nvfortran, GNU gfortran, and Intel ifx/ifort compilers

src/simulation/**/*.{fpp,f90}: Mark GPU-callable helpers with $:GPU_ROUTINE(function_name='...', parallelism='[seq]') immediately after declaration
Do not use OpenACC or OpenMP directives directly; use Fypp macros from src/common/include/parallel_macros.fpp instead
Wrap tight loops with $:GPU_PARALLEL_FOR(private='[...]', copy='[...]') macro; add collapse=n for safe nested loop merging
Declare loop-local variables with private='[...]' in GPU parallel loop macros
Allocate large arrays with managed or move them into a persistent $:GPU_ENTER_DATA(...) region at start-up
Do not place stop or error stop inside device code

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
src/**/*.fpp

📄 CodeRabbit inference engine (.cursor/rules/mfc-agent-rules.mdc)

src/**/*.fpp: Use .fpp file extension for Fypp preprocessed files; CMake transpiles them to .f90
Start module files with Fypp include for macros: #:include 'macros.fpp'
Use the fypp ASSERT macro for validating conditions: @:ASSERT(predicate, message)
Use fypp macro @:ALLOCATE(var1, var2) for device-aware allocation instead of standard Fortran allocate
Use fypp macro @:DEALLOCATE(var1, var2) for device-aware deallocation instead of standard Fortran deallocate

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
🧠 Learnings (9)
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Mark GPU-callable helpers with `$:GPU_ROUTINE(function_name='...', parallelism='[seq]')` immediately after declaration

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Wrap tight GPU loops with !$acc parallel loop gang vector default(present) reduction(...); add collapse(n) when safe; declare loop-local variables with private(...)

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Wrap tight loops with `$:GPU_PARALLEL_FOR(private='[...]', copy='[...]')` macro; add `collapse=n` for safe nested loop merging

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Ensure GPU code compiles with Cray ftn, NVIDIA nvfortran, GNU gfortran, and Intel ifx/ifort compilers

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Avoid stop/error stop inside GPU device code

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Declare loop-local variables with `private='[...]'` in GPU parallel loop macros

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Allocate large GPU arrays with managed memory or move them into persistent !$acc enter data regions at start-up

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Allocate large arrays with `managed` or move them into a persistent `$:GPU_ENTER_DATA(...)` region at start-up

Applied to files:

  • src/simulation/m_igr.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Do not use OpenACC or OpenMP directives directly; use Fypp macros from `src/common/include/parallel_macros.fpp` instead

Applied to files:

  • src/simulation/m_igr.fpp
⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (1)
  • GitHub Check: cubic · AI code reviewer
🔇 Additional comments (1)
src/simulation/m_ibm.fpp (1)

269-309: Index formulations for species volume fractions are correctly synchronized

The code correctly uses both index forms interchangeably: IGR writes to q_cons_vf(E_idx + q) while non-IGR writes to q_prim_vf(advxb + q - 1), and since advxb = E_idx + 1, both expressions reference identical array positions. The conservative variable assignments stay in sync across all branches.

Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot left a comment

Choose a reason for hiding this comment

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

4 issues found across 3 files

Prompt for AI agents (all 4 issues)

Check if these issues are valid — if so, understand the root cause of each and fix them.


<file name="src/simulation/m_ibm.fpp">

<violation number="1" location="src/simulation/m_ibm.fpp:162">
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.</violation>

<violation number="2" location="src/simulation/m_ibm.fpp:240">
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.</violation>

<violation number="3" location="src/simulation/m_ibm.fpp:291">
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.</violation>

<violation number="4" location="src/simulation/m_ibm.fpp:982">
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.</violation>
</file>

Reply to cubic to teach it or ask questions. Re-run a review with @cubic-dev-ai review this PR


! 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

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

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
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Dec 19, 2025

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

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

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

Comment on lines +982 to +983
! 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)
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Dec 19, 2025

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

Copy link
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 2

Caution

Some comments are outside the diff and can’t be posted inline due to platform limitations.

⚠️ Outside diff range comments (1)
src/simulation/m_ibm.fpp (1)

898-927: Uninitialized mixture variables when igr and elasticity are both true

The elasticity branch in the IGR path of s_interpolate_image_point remains commented out (lines 988–989), leaving rho_K, gamma_K, pi_inf_K, and qv_K uninitialized when both igr and elasticity are enabled. These uninitialized variables are then used in:

  • rho_K = max(rho_K, sgm_eps) (line 997)
  • Dynamic pressure accumulation with division by rho_K (lines 1003–1006)
  • s_compute_pressure call passing pi_inf_K, gamma_K, rho_K, qv_K (lines 1008–1012)

This produces undefined behavior and garbage values for pressure and velocity at IBM ghost points.

Enable the elasticity path and declare the missing variables G_K (scalar) and Gs (array of size num_fluids) to match the pattern used in other routines like m_sim_helpers.fpp.

🧹 Nitpick comments (2)
src/simulation/m_ibm.fpp (2)

244-257: Moving-IBM interior pressure is now explicitly disabled for IGR

Gating the “moving IBM interior pressure” initialization with if (.not. igr) makes the current limitation (moving IBM + IGR not supported) explicit and prevents that heuristic from running under IGR.

This is reasonable as an interim step; just ensure this limitation is documented or guarded elsewhere (e.g., in case setup/validator) so users don’t assume full moving-IBM support with IGR.


898-905: Add guard for optional q_cons_vf when IGR is enabled

The subroutine s_interpolate_image_point declares q_cons_vf as optional but dereferences it unconditionally inside the if (igr) block without checking present(). While the current only IGR call site (s_ibm_correct_state) always supplies q_cons_vf, this is fragile. Future callers might invoke this routine under IGR without providing q_cons_vf, causing undefined behavior.

Add a guard at the routine entry:

if (igr .and. .not. present(q_cons_vf)) then
    call s_mpi_abort('s_interpolate_image_point: q_cons_vf must be present when igr is enabled')
end if

This hardens the API boundary and prevents silent errors.

Also applies to: 970-987

📜 Review details

Configuration used: defaults

Review profile: CHILL

Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between 1340497 and b07045b.

📒 Files selected for processing (2)
  • src/simulation/m_ibm.fpp (11 hunks)
  • src/simulation/m_igr.fpp (2 hunks)
🧰 Additional context used
📓 Path-based instructions (3)
**/*.{fpp,f90}

📄 CodeRabbit inference engine (.github/copilot-instructions.md)

**/*.{fpp,f90}: Use 2-space indentation; continuation lines align beneath &
Use lower-case keywords and intrinsics (do, end subroutine, etc.)
Name modules with m_ pattern (e.g., m_transport)
Name public subroutines with s_ pattern (e.g., s_compute_flux)
Name public functions with f
_ pattern
Keep subroutine size ≤ 500 lines, helper subroutines ≤ 150 lines, functions ≤ 100 lines, files ≤ 1000 lines
Limit routine arguments to ≤ 6; use derived-type params struct if more are needed
Forbid goto statements (except in legacy code), COMMON blocks, and save globals
Every argument must have explicit intent; use dimension/allocatable/pointer as appropriate
Call s_mpi_abort() for errors, never use stop or error stop

**/*.{fpp,f90}: Indent 2 spaces; continuation lines align under &
Use lower-case keywords and intrinsics (do, end subroutine, etc.)
Name modules with m_<feature> prefix (e.g., m_transport)
Name public subroutines as s_<verb>_<noun> (e.g., s_compute_flux) and functions as f_<verb>_<noun>
Keep private helpers in the module; avoid nested procedures
Enforce size limits: subroutine ≤ 500 lines, helper ≤ 150, function ≤ 100, module/file ≤ 1000
Limit subroutines to ≤ 6 arguments; otherwise pass a derived-type 'params' struct
Avoid goto statements (except unavoidable legacy); avoid global state (COMMON, save)
Every variable must have intent(in|out|inout) specification and appropriate dimension / allocatable / pointer
Use s_mpi_abort(<msg>) for error termination instead of stop
Use !> style documentation for header comments; follow Doxygen Fortran format with !! @param and !! @return tags
Use implicit none statement in all modules
Use private declaration followed by explicit public exports in modules
Use derived types with pointers for encapsulation (e.g., pointer, dimension(:,:,:) => null())
Use pure and elemental attributes for side-effect-free functions; combine them for array ...

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
src/simulation/**/*.{fpp,f90}

📄 CodeRabbit inference engine (.github/copilot-instructions.md)

src/simulation/**/*.{fpp,f90}: Wrap tight GPU loops with !$acc parallel loop gang vector default(present) reduction(...); add collapse(n) when safe; declare loop-local variables with private(...)
Allocate large GPU arrays with managed memory or move them into persistent !$acc enter data regions at start-up
Avoid stop/error stop inside GPU device code
Ensure GPU code compiles with Cray ftn, NVIDIA nvfortran, GNU gfortran, and Intel ifx/ifort compilers

src/simulation/**/*.{fpp,f90}: Mark GPU-callable helpers with $:GPU_ROUTINE(function_name='...', parallelism='[seq]') immediately after declaration
Do not use OpenACC or OpenMP directives directly; use Fypp macros from src/common/include/parallel_macros.fpp instead
Wrap tight loops with $:GPU_PARALLEL_FOR(private='[...]', copy='[...]') macro; add collapse=n for safe nested loop merging
Declare loop-local variables with private='[...]' in GPU parallel loop macros
Allocate large arrays with managed or move them into a persistent $:GPU_ENTER_DATA(...) region at start-up
Do not place stop or error stop inside device code

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
src/**/*.fpp

📄 CodeRabbit inference engine (.cursor/rules/mfc-agent-rules.mdc)

src/**/*.fpp: Use .fpp file extension for Fypp preprocessed files; CMake transpiles them to .f90
Start module files with Fypp include for macros: #:include 'macros.fpp'
Use the fypp ASSERT macro for validating conditions: @:ASSERT(predicate, message)
Use fypp macro @:ALLOCATE(var1, var2) for device-aware allocation instead of standard Fortran allocate
Use fypp macro @:DEALLOCATE(var1, var2) for device-aware deallocation instead of standard Fortran deallocate

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
🧠 Learnings (10)
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Ensure GPU code compiles with Cray ftn, NVIDIA nvfortran, GNU gfortran, and Intel ifx/ifort compilers

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Wrap tight GPU loops with !$acc parallel loop gang vector default(present) reduction(...); add collapse(n) when safe; declare loop-local variables with private(...)

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Wrap tight loops with `$:GPU_PARALLEL_FOR(private='[...]', copy='[...]')` macro; add `collapse=n` for safe nested loop merging

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Mark GPU-callable helpers with `$:GPU_ROUTINE(function_name='...', parallelism='[seq]')` immediately after declaration

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Do not use OpenACC or OpenMP directives directly; use Fypp macros from `src/common/include/parallel_macros.fpp` instead

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Declare loop-local variables with `private='[...]'` in GPU parallel loop macros

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to **/*.{fpp,f90} : Use `private` declaration followed by explicit `public` exports in modules

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Avoid stop/error stop inside GPU device code

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Allocate large arrays with `managed` or move them into a persistent `$:GPU_ENTER_DATA(...)` region at start-up

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Allocate large GPU arrays with managed memory or move them into persistent !$acc enter data regions at start-up

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (1)
  • GitHub Check: cubic · AI code reviewer
🔇 Additional comments (2)
src/simulation/m_igr.fpp (1)

373-398: IBM post-processing of jac and jac_old looks consistent

Zeroing jac inside IBM-marked cells over 0:m,0:n,0:p, then calling s_update_igr(jac_sf) and finally re-syncing jac_old for Jacobi (igr_iter_solver == 1) matches the allocation/usage patterns of jac/jac_old and the IBM ghost-point interpolation.

No further issues spotted in this block.

src/simulation/m_ibm.fpp (1)

289-307: IGR-specific conservative updates in s_ibm_correct_state look coherent

For igr runs:

  • The image-point interpolation now fills alpha_rho_IP/alpha_IP via the IGR path in s_interpolate_image_point.
  • You correctly map those into q_cons_vf (single- and multi-fluid) at the ghost points, instead of going through q_prim_vf.
  • The later .not. igr block that resets continuity and advective variables is skipped, so the IGR-consistent values aren’t overwritten.

This matches the new conservative-variable-based design; no issues spotted.

Also applies to: 390-397

Comment on lines 153 to 196
!> Interpolates sigma from the m_igr module at all ghost points
!! @param jac_sf Sigma, Entropic pressure
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

if (p == 0) then
l1 = 0
l2 = 0
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
Copy link
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Tighten 2D handling in s_update_igr to avoid reading an unused ip_grid(3)

The new 2D guard (if (p == 0) then l1 = 0; l2 = 0) correctly prevents looping over an undefined k-stencil, but gp%ip_grid(3) is still read unconditionally before being overridden. In 2D runs, that component is never set by s_compute_image_points, so this is technically an uninitialized read even though the value is not used to form loop bounds.

You can avoid touching ip_grid(3) entirely when p == 0:

Suggested change
-        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
+        j1 = gp%ip_grid(1); j2 = j1 + 1
+        k1 = gp%ip_grid(2); k2 = k1 + 1
+        if (p == 0) then
+            l1 = 0
+            l2 = 0
+        else
+            l1 = gp%ip_grid(3)
+            l2 = l1 + 1
+        end if
🤖 Prompt for AI Agents
In src/simulation/m_ibm.fpp around lines 153-196, the code reads gp%ip_grid(3)
unconditionally which is uninitialized in 2D (p==0); change the order so
ip_grid(1) and ip_grid(2) are read as now, but only read gp%ip_grid(3) and set
l1,l2 from it when p /= 0; when p == 0 set l1 = 0 and l2 = 0 without touching
ip_grid(3). Ensure any uses of l1/l2 (loops and interp_coeffs indexing) still
work by guaranteeing l1/l2 are defined in both branches.


use m_boundary_common

use m_ibm, only: ib, ib_markers, s_update_igr
Copy link
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🔴 Critical

Fix use m_ibm ONLY-list: ib is not exported here and likely conflicts with m_global_parameters

m_ibm (as shown) does not declare an ib symbol, while m_igr already has ib via use m_global_parameters. Importing ib from m_ibm will either fail to compile or create an ambiguous symbol if ib is also exported from m_global_parameters.

Rely on ib from m_global_parameters and restrict the m_ibm import to the IBM-specific data and routine:

Suggested change
-    use m_ibm, only: ib, ib_markers, s_update_igr
+    use m_ibm, only: ib_markers, s_update_igr
🤖 Prompt for AI Agents
In src/simulation/m_igr.fpp around line 18, the "use m_ibm, only: ib,
ib_markers, s_update_igr" clause is wrong because m_ibm does not export "ib"
(and "ib" is already provided by m_global_parameters); remove "ib" from the
only-list so you only import the IBM-specific symbols (e.g., "ib_markers,
s_update_igr") and rely on the existing "ib" from m_global_parameters to avoid
compile errors or symbol ambiguity.

@codecov
Copy link

codecov bot commented Dec 19, 2025

Codecov Report

❌ Patch coverage is 12.50000% with 84 lines in your changes missing coverage. Please review.
✅ Project coverage is 43.84%. Comparing base (8ae42aa) to head (489bc7c).
⚠️ Report is 3 commits behind head on master.

Files with missing lines Patch % Lines
src/simulation/m_ibm.fpp 14.28% 68 Missing and 4 partials ⚠️
src/simulation/m_igr.fpp 0.00% 11 Missing and 1 partial ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master    #1095      +/-   ##
==========================================
- Coverage   44.08%   43.84%   -0.25%     
==========================================
  Files          71       71              
  Lines       20332    20364      +32     
  Branches     1981     1991      +10     
==========================================
- Hits         8963     8928      -35     
- Misses      10236    10296      +60     
- Partials     1133     1140       +7     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

!> Interpolates sigma from the m_igr module at all ghost points
!! @param jac_sf Sigma, Entropic pressure
subroutine s_update_igr(jac_sf)
type(scalar_field), dimension(1), intent(inout) :: jac_sf
Copy link
Contributor

Choose a reason for hiding this comment

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

The frontier error might go away if you replace jac_sf here with the non-pointer variable jac from m_igr

end do
end do
$:END_GPU_PARALLEL_LOOP()
call s_update_igr(jac_sf)
Copy link
Contributor

Choose a reason for hiding this comment

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

See previous comment. Passing by jac rather than jac_sf is probably preferred. I only made the jac_sf variable so that I could reuse MPI routines elsewhere in this module.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

ok, got it

@wilfonba
Copy link
Contributor

@Cowsreal It would be nice to see a convergence result for IGR vs WENO. I.e., run the supersonic flow over a cylinder case at a moderately high resolution with WENO5, then run IGR at three or four resolutions approaching that moderately high resolution, and compare the surface pressure along the same axis. Depending on how automated/quick this is, it might be good to see the same for the airfoil (or do the airfoil instead. Since they're NACA profiles, you could compute the pressure lift/drag and compare it to published values).

@Cowsreal Cowsreal closed this Dec 19, 2025
@Cowsreal
Copy link
Contributor Author

@Cowsreal It would be nice to see a convergence result for IGR vs WENO. I.e., run the supersonic flow over a cylinder case at a moderately high resolution with WENO5, then run IGR at three or four resolutions approaching that moderately high resolution, and compare the surface pressure along the same axis. Depending on how automated/quick this is, it might be good to see the same for the airfoil (or do the airfoil instead. Since they're NACA profiles, you could compute the pressure lift/drag and compare it to published values).

Yes, ok that would be cool to find a real study and compare it. I'll try to do that.

@Cowsreal Cowsreal reopened this Dec 19, 2025
@qodo-code-review
Copy link
Contributor

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

⏱️ Estimated effort to review: 4 🔵🔵🔵🔵⚪
🧪 No relevant tests
🔒 No security concerns identified
⚡ Recommended focus areas for review

Possible Issue

The new IGR branch in s_interpolate_image_point relies on the optional q_cons_vf argument but does not guard its use with present(q_cons_vf). If igr is enabled and any call site forgets to pass q_cons_vf, this will dereference an absent optional argument at runtime. Consider enforcing the argument for IGR (separate interface/subroutine) or adding a present() check with a clear error path.

    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, 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

        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
        real(wp), optional, dimension(:), intent(INOUT) :: presb_IP, massv_IP

        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
        k1 = gp%ip_grid(3); k2 = k1 + 1

        if (p == 0) then
            k1 = 0
            k2 = 0
        end if

        alpha_rho_IP = 0._wp
        alpha_IP = 0._wp
        pres_IP = 0._wp
        vel_IP = 0._wp
        pres = 0._wp

        if (surface_tension) c_IP = 0._wp

        if (bubbles_euler) then
            r_IP = 0._wp
            v_IP = 0._wp
            if (.not. polytropic) then
                mv_IP = 0._wp
                pb_IP = 0._wp
            end if
        end if

        if (qbmm) then
            nmom_IP = 0._wp
            if (.not. polytropic) then
                presb_IP = 0._wp
                massv_IP = 0._wp
            end if
        end if

        $:GPU_LOOP(parallelism='[seq]')
        do i = i1, i2
            $:GPU_LOOP(parallelism='[seq]')
            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)

                    if (igr) then
                        ! 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)
                            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

                        rho_K = max(rho_K, sgm_eps)

                        $: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 (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
Known Limitation

The ghost-point pressure update is explicitly skipped for IGR (.not. igr) with a comment indicating moving-IGR is not handled. This creates a high risk of silently incorrect results for cases with igr enabled and moving immersed boundaries, since ghost-point pressure/energy consistency may break. At minimum, consider failing fast (validation/runtime) when igr and moving IBM are both enabled until implemented, or implement the missing path.

! set the pressure
! !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
Algorithm Risk

The solver now zeroes jac at IBM interior points and then calls s_update_igr to overwrite ghost points via image-point interpolation. This changes the iteration state mid-solve and (for Jacobi) conditionally re-syncs jac_old. Please double-check that: (1) the interior/ghost classification via ib_markers matches ghost_points locations, (2) jac_old consistency is preserved for both solver modes, and (3) MPI halo/BC updates aren’t required between the zeroing and interpolation (otherwise boundary-adjacent ghost points may use stale jac_sf values).

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

Copy link
Contributor

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

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

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

Comment on lines 373 to 398
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
Copy link
Contributor

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, when updating jac_old for the Jacobi solver, use loop bounds 0:m, 0:n, 0:p to match the bounds used when modifying jac, ensuring consistency. [general, importance: 6]

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_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
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
!$acc parallel loop collapse(3) private(j,k,l) copyin(jac) copyout(jac_old)
do l = 0, p
do k = 0, n
do j = 0, m
jac_old(j, k, l) = jac(j, k, l)
end do
end do
end do
end if
end if

!> Interpolates sigma from the m_igr module at all ghost points
!! @param jac_sf Sigma, Entropic pressure
subroutine s_update_igr(jac_sf)
type(scalar_field), dimension(1), intent(inout) :: jac_sf
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, declare the dummy argument jac_sf as an assumed-shape array (dimension(:)) instead of a fixed-size array (dimension(1)). [possible issue, importance: 7]

Suggested change
type(scalar_field), dimension(1), intent(inout) :: jac_sf
type(scalar_field), dimension(:), intent(inout) :: jac_sf

Comment on lines +970 to +972
if (igr) then
! For IGR, we will need to perform operations on
! the conservative variables instead
Copy link
Contributor

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

@codeant-ai
Copy link

codeant-ai bot commented Dec 19, 2025

CodeAnt AI is running Incremental review


Thanks for using CodeAnt! 🎉

We're free for open-source projects. if you're enjoying it, help us grow by sharing.

Share on X ·
Reddit ·
LinkedIn

Copy link
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 0

♻️ Duplicate comments (3)
src/simulation/m_ibm.fpp (2)

153-196: Avoid reading unused ip_grid(3) in 2D in s_interpolate_sigma_igr

In 2D (p == 0), gp%ip_grid(3) is not populated by s_compute_image_points, yet it’s still read before being overridden:

l1 = gp%ip_grid(3); l2 = l1 + 1
if (p == 0) then
    l1 = 0
    l2 = 0
end if

Restructure to only access ip_grid(3) when p /= 0:

Proposed change
-                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
+                j1 = gp%ip_grid(1); j2 = j1 + 1
+                k1 = gp%ip_grid(2); k2 = k1 + 1
+                if (p == 0) then
+                    l1 = 0
+                    l2 = 0
+                else
+                    l1 = gp%ip_grid(3)
+                    l2 = l1 + 1
+                end if

This keeps the current behavior while eliminating an unnecessary read of a potentially undefined component in 2D.


898-905: IGR path in s_interpolate_image_point: uninitialized mixture vars with elasticity, and brittle use of optional q_cons_vf

Two issues in the new IGR branch need tightening:

  1. Elasticity + IGR leaves mixture variables uninitialized (real bug)
    In the igr block:

    if (model_eqns /= 4) then
        if (elasticity) then

! call s_convert_species_to_mixture_variables_acc(...)
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

rho_K = max(rho_K, sgm_eps)
...
call s_compute_pressure(..., pi_inf_K, gamma_K, rho_K, qv_K, ...)


When `elasticity` is true, the conversion call is commented out, so `rho_K`, `gamma_K`, `pi_inf_K`, and `qv_K` are never set before use. This is undefined behavior for any elastic IGR case.

You likely want to mirror the non‑IGR path and actually call the elastic variant here, e.g.:

<details>
<summary>Suggested fix for elasticity branch</summary>

```diff
-                        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
+                        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
+                        end if

(Adjust the final arguments to match the actual elastic variant’s signature used elsewhere in this module.)

  1. Optional q_cons_vf is dereferenced unconditionally when igr is true
    The igr branch assumes q_cons_vf is present:

    if (igr) then
        ...
        alpha_rho_K(1) = q_cons_vf(1)%sf(i, j, k)
        ...
        q_cons_vf(E_idx + l)%sf(i, j, k)
        ...
    end if

    but q_cons_vf is declared optional. Right now, s_ibm_correct_state passes q_cons_vf only in the else if (igr) branch; other call sites (bubbles/qbmm) do not. If IGR is ever combined with those modes, this routine will execute with igr == .true. and an absent q_cons_vf, causing undefined behavior.

    To make this robust:

    • Either make q_cons_vf a required argument (simplest and safest given all current usages are in IBM/IGR), or
    • Add an explicit guard and treat missing q_cons_vf as a configuration error, e.g. check at the caller level (host side) that igr implies q_cons_vf is supplied before entering the GPU loop.

    In any case, relying on a global igr flag with an optional argument here is fragile and easy to break with future call sites.

Also applies to: 921-927, 970-1053

src/simulation/m_igr.fpp (1)

18-18: Fix use m_ibm ONLY-list to avoid invalid/ambiguous ib import

m_ibm does not export ib, and this module already gets ib from m_global_parameters. Keeping ib in the ONLY-list will either fail compilation or create ambiguity. Restrict the import to IBM-specific symbols:

Proposed change
-    use m_ibm, only: ib, ib_markers, s_interpolate_sigma_igr
+    use m_ibm, only: ib_markers, s_interpolate_sigma_igr
🧹 Nitpick comments (1)
src/simulation/m_ibm.fpp (1)

289-309: IGR conservative updates at ghost points look consistent; moving‑IB+IGR remains a known TODO

  • The new else if (igr) branch routes image-point interpolation through the conservative q_cons_vf path and writes back alpha_rho_IP/alpha_IP into q_cons_vf at ghost points. This avoids reusing the primitive-variable path and is coherent with how q_cons_vf(E_idx) is formed later.
  • The .not. igr guards around primitive pressure assignment and continuity/adv updates prevent double-writing in the IGR case.

The !TEMPORARY, NEED TO FIX FOR MOVING IGR comment still signals that moving immersed boundaries with IGR are not fully supported; that’s an existing limitation rather than a new regression, but worth keeping documented in user-facing docs or validators.

Also applies to: 390-397, 322-335

📜 Review details

Configuration used: defaults

Review profile: CHILL

Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between b07045b and 489bc7c.

📒 Files selected for processing (2)
  • src/simulation/m_ibm.fpp (11 hunks)
  • src/simulation/m_igr.fpp (2 hunks)
🧰 Additional context used
📓 Path-based instructions (3)
**/*.{fpp,f90}

📄 CodeRabbit inference engine (.github/copilot-instructions.md)

**/*.{fpp,f90}: Use 2-space indentation; continuation lines align beneath &
Use lower-case keywords and intrinsics (do, end subroutine, etc.)
Name modules with m_ pattern (e.g., m_transport)
Name public subroutines with s_ pattern (e.g., s_compute_flux)
Name public functions with f
_ pattern
Keep subroutine size ≤ 500 lines, helper subroutines ≤ 150 lines, functions ≤ 100 lines, files ≤ 1000 lines
Limit routine arguments to ≤ 6; use derived-type params struct if more are needed
Forbid goto statements (except in legacy code), COMMON blocks, and save globals
Every argument must have explicit intent; use dimension/allocatable/pointer as appropriate
Call s_mpi_abort() for errors, never use stop or error stop

**/*.{fpp,f90}: Indent 2 spaces; continuation lines align under &
Use lower-case keywords and intrinsics (do, end subroutine, etc.)
Name modules with m_<feature> prefix (e.g., m_transport)
Name public subroutines as s_<verb>_<noun> (e.g., s_compute_flux) and functions as f_<verb>_<noun>
Keep private helpers in the module; avoid nested procedures
Enforce size limits: subroutine ≤ 500 lines, helper ≤ 150, function ≤ 100, module/file ≤ 1000
Limit subroutines to ≤ 6 arguments; otherwise pass a derived-type 'params' struct
Avoid goto statements (except unavoidable legacy); avoid global state (COMMON, save)
Every variable must have intent(in|out|inout) specification and appropriate dimension / allocatable / pointer
Use s_mpi_abort(<msg>) for error termination instead of stop
Use !> style documentation for header comments; follow Doxygen Fortran format with !! @param and !! @return tags
Use implicit none statement in all modules
Use private declaration followed by explicit public exports in modules
Use derived types with pointers for encapsulation (e.g., pointer, dimension(:,:,:) => null())
Use pure and elemental attributes for side-effect-free functions; combine them for array ...

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
src/simulation/**/*.{fpp,f90}

📄 CodeRabbit inference engine (.github/copilot-instructions.md)

src/simulation/**/*.{fpp,f90}: Wrap tight GPU loops with !$acc parallel loop gang vector default(present) reduction(...); add collapse(n) when safe; declare loop-local variables with private(...)
Allocate large GPU arrays with managed memory or move them into persistent !$acc enter data regions at start-up
Avoid stop/error stop inside GPU device code
Ensure GPU code compiles with Cray ftn, NVIDIA nvfortran, GNU gfortran, and Intel ifx/ifort compilers

src/simulation/**/*.{fpp,f90}: Mark GPU-callable helpers with $:GPU_ROUTINE(function_name='...', parallelism='[seq]') immediately after declaration
Do not use OpenACC or OpenMP directives directly; use Fypp macros from src/common/include/parallel_macros.fpp instead
Wrap tight loops with $:GPU_PARALLEL_FOR(private='[...]', copy='[...]') macro; add collapse=n for safe nested loop merging
Declare loop-local variables with private='[...]' in GPU parallel loop macros
Allocate large arrays with managed or move them into a persistent $:GPU_ENTER_DATA(...) region at start-up
Do not place stop or error stop inside device code

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
src/**/*.fpp

📄 CodeRabbit inference engine (.cursor/rules/mfc-agent-rules.mdc)

src/**/*.fpp: Use .fpp file extension for Fypp preprocessed files; CMake transpiles them to .f90
Start module files with Fypp include for macros: #:include 'macros.fpp'
Use the fypp ASSERT macro for validating conditions: @:ASSERT(predicate, message)
Use fypp macro @:ALLOCATE(var1, var2) for device-aware allocation instead of standard Fortran allocate
Use fypp macro @:DEALLOCATE(var1, var2) for device-aware deallocation instead of standard Fortran deallocate

Files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
🧠 Learnings (12)
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Ensure GPU code compiles with Cray ftn, NVIDIA nvfortran, GNU gfortran, and Intel ifx/ifort compilers

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Mark GPU-callable helpers with `$:GPU_ROUTINE(function_name='...', parallelism='[seq]')` immediately after declaration

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Wrap tight GPU loops with !$acc parallel loop gang vector default(present) reduction(...); add collapse(n) when safe; declare loop-local variables with private(...)

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Wrap tight loops with `$:GPU_PARALLEL_FOR(private='[...]', copy='[...]')` macro; add `collapse=n` for safe nested loop merging

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Declare loop-local variables with `private='[...]'` in GPU parallel loop macros

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Allocate large arrays with `managed` or move them into a persistent `$:GPU_ENTER_DATA(...)` region at start-up

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Do not use OpenACC or OpenMP directives directly; use Fypp macros from `src/common/include/parallel_macros.fpp` instead

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Avoid stop/error stop inside GPU device code

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:16.713Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .github/copilot-instructions.md:0-0
Timestamp: 2025-11-24T21:50:16.713Z
Learning: Applies to src/simulation/**/*.{fpp,f90} : Allocate large GPU arrays with managed memory or move them into persistent !$acc enter data regions at start-up

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to **/*.{fpp,f90} : Use `private` declaration followed by explicit `public` exports in modules

Applied to files:

  • src/simulation/m_igr.fpp
  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to **/*.{fpp,f90} : Avoid `goto` statements (except unavoidable legacy); avoid global state (`COMMON`, `save`)

Applied to files:

  • src/simulation/m_ibm.fpp
📚 Learning: 2025-11-24T21:50:46.909Z
Learnt from: CR
Repo: MFlowCode/MFC PR: 0
File: .cursor/rules/mfc-agent-rules.mdc:0-0
Timestamp: 2025-11-24T21:50:46.909Z
Learning: Applies to **/*.{fpp,f90} : Use `wp` (working precision) parameter from `m_precision_select` instead of hardcoded precision like `real*8`

Applied to files:

  • src/simulation/m_ibm.fpp
⏰ Context from checks skipped due to timeout of 90000ms. You can increase the timeout in your CodeRabbit configuration to a maximum of 15 minutes (900000ms). (17)
  • GitHub Check: Code Cleanliness Check
  • GitHub Check: Coverage Test on CodeCov
  • GitHub Check: Github (ubuntu, no-mpi, single, no-debug, false)
  • GitHub Check: Github (macos, mpi, no-debug, false)
  • GitHub Check: Oak Ridge | Frontier (gpu-acc)
  • GitHub Check: Github (macos, mpi, debug, false)
  • GitHub Check: Github (ubuntu, mpi, no-debug, false)
  • GitHub Check: Github (ubuntu, mpi, debug, false)
  • GitHub Check: Oak Ridge | Frontier (cpu)
  • GitHub Check: Oak Ridge | Frontier (gpu-omp)
  • GitHub Check: Github (ubuntu, mpi, no-debug, true)
  • GitHub Check: Georgia Tech | Phoenix (gpu-omp)
  • GitHub Check: Georgia Tech | Phoenix (gpu-acc)
  • GitHub Check: Github (ubuntu, mpi, debug, true)
  • GitHub Check: Georgia Tech | Phoenix (cpu)
  • GitHub Check: cubic · AI code reviewer
  • GitHub Check: Build & Publish
🔇 Additional comments (2)
src/simulation/m_igr.fpp (1)

373-398: IBM–IGR post-solve integration looks consistent

Zeroing jac at IB-marked cells, then calling s_interpolate_sigma_igr(jac) and refreshing jac_old for Jacobi in GPU-parallel loops is consistent with existing patterns (correct privates, collapse, and domains). This is a reasonable place in the iteration loop to apply the IBM correction.

Please sanity-check a 2D IBM+IGR case (p = 0) to confirm s_interpolate_sigma_igr gives expected ghost-point sigma after this block.

src/simulation/m_ibm.fpp (1)

244-257: GPU loop for non‑IGR interior pressure initialization is now well-formed

The .not. igr GPU loop correctly marks IB cells in q_prim_vf(E_idx) and now lists all loop indices (j,k,l) as private, matching the project’s GPU macro pattern and avoiding races.

@sbryngelson
Copy link
Member

the results make look the same but we need to be sure (1) it's taking the IGR code path (not WENO). are the results identical? and (2) we need to see it working for a 'real' case. This case has a discontinuity in the initial condition at the airfoil location that causes spurious artifacts, it seems. we should see real vortex shedding

@sbryngelson
Copy link
Member

the AI comments seem relevant, or at least several of them, as well

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Review effort 4/5 size:L This PR changes 100-499 lines, ignoring generated files

Development

Successfully merging this pull request may close these issues.

3 participants