Skip to content

Conversation

@danieljvickers
Copy link
Member

@danieljvickers danieljvickers commented Jan 15, 2026

User description

User description

Description

Adds a numeric method of computing the centroid of IBs where the center of mass is not well posed. In particular, I added hooks for STL and airfiol IBs.

This PR also includes support for parameterized IB motion. It allows one to analytically describe the velocity and rotation rate of any IB experiencing 1-way coupling.

Both of these features have been tested with IB's the move and rotate using the new analytic parameterization. I have also tested with an airfoil rotating about the correct center of mass. Videos available in my slides.

Fixes #1105
Fixes #1106

Type of change

Please delete options that are not relevant.

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)

Scope

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

That goal was enabling the ability to run a pitching airfoil case, which oscillates in time.

How Has This Been Tested?

  • I ran a pitching airfoil case which oscillates in time at a fixed frequency about the new correct center of mass
  • I ran a case of a cylinder moving in a circular pattern to "stir" a fluid

The best configuration can be summarized with this test case file for the pitching airfoil:

import json
import math

Mu = 1.84e-05
gam_a = 1.4
gam_b = 1.1

total_time = 5.0 
num_time_steps = 2000
num_saves = 500
dt = total_time / float(num_time_steps)
steps_per_save = int(num_time_steps / num_saves)

# Configuring case dictionary
print(
    json.dumps(
        {
            # Logistics
            "run_time_info": "T",
            # Computational Domain Parameters
            # axial direction
            "x_domain%beg": 0.0e00,
            "x_domain%end": 6.0,
            # r direction
            "y_domain%beg": 0.0e00,
            "y_domain%end": 3.0,
            "cyl_coord": "F",
            "m": 400,
            "n": 200,
            "p": 0,
            "dt": dt,
            "t_step_start": 0,
            "t_step_stop": num_time_steps,
            "t_step_save": steps_per_save,
            # Simulation Algorithm Parameters
            # Only one patches are necessary, the air tube
            "num_patches": 1,
            # Use the 5 equation model
            "model_eqns": 2,
            # 6 equations model does not need the K \div(u) term
            "alt_soundspeed": "F",
            # One fluids: air
            "num_fluids": 2,
            # time step
            "mpp_lim": "F",
            # Correct errors when computing speed of sound
            "mixture_err": "T",
            # Use TVD RK3 for time marching
            "time_stepper": 3,
            # Reconstruct the primitive variables to minimize spurious
            # Use WENO5
            "weno_order": 5,
            "weno_eps": 1.0e-16,
            "weno_Re_flux": "T",
            "weno_avg": "T",
            "avg_state": 2,
            "viscous": "T",
            # Use the mapped WENO weights to maintain monotinicity
            "mapped_weno": "T",
            "null_weights": "F",
            "mp_weno": "T",
            # Use the HLLC  Riemann solver
            "riemann_solver": 2,
            "wave_speeds": 1,
            # We use reflective boundary conditions at octant edges and
            # non-reflective boundary conditions at the domain edges
            "bc_x%beg": -3,
            "bc_x%end": -3,
            "bc_y%beg": -3,
            "bc_y%end": -3,
            # Set IB to True and add 1 patch
            "ib": "T",
            "num_ibs": 1,
            # Formatted Database Files Structure Parameters
            "format": 1,
            "precision": 2,
            "prim_vars_wrt": "T",
            "E_wrt": "T",
            "parallel_io": "T",
            # Patch: Constant Tube filled with air
            # Specify the cylindrical air tube grid geometry
            "patch_icpp(1)%geometry": 3,
            "patch_icpp(1)%x_centroid": 3.0,
            # Uniform medium density, centroid is at the center of the domain
            "patch_icpp(1)%y_centroid": 1.5,
            "patch_icpp(1)%length_x": 6.0,
            "patch_icpp(1)%length_y": 3.0,
            # Specify the patch primitive variables
            "patch_icpp(1)%vel(1)": 3.0e00,
            "patch_icpp(1)%vel(2)": 0.0e00,
            "patch_icpp(1)%pres": 1.0e00,
            "patch_icpp(1)%alpha_rho(1)": 0.8e00,
            "patch_icpp(1)%alpha(1)": 0.8e00,
            "patch_icpp(1)%alpha_rho(2)": 0.2e00,
            "patch_icpp(1)%alpha(2)": 0.2e00,
            # Patch: Airfoil Immersed Boundary
            "patch_ib(1)%geometry": 4,
            "patch_ib(1)%x_centroid": 1.5,
            "patch_ib(1)%y_centroid": 1.5,
            "patch_ib(1)%c": 1.0,
            "patch_ib(1)%t": 0.15,
            "patch_ib(1)%p": 0.4,
            "patch_ib(1)%m": 0.02,
            "patch_ib(1)%angular_vel(3)": "-10.0 * 2.0 * pi * sin(2.0 * pi * t) * pi / 180.",
            "patch_ib(1)%angles(3)": -5.0 * math.pi / 180.,
            "patch_ib(1)%moving_ibm": 1,
            "patch_ib(1)%slip": "F",
            # Fluids Physical Parameters
            # Use the same stiffness as the air bubble
            "fluid_pp(1)%gamma": 1.0e00 / (gam_a - 1.0e00),  # 2.50(Not 1.40)
            "fluid_pp(1)%pi_inf": 0,
            "fluid_pp(1)%Re(1)": 1000000,
            "fluid_pp(2)%gamma": 1.0e00 / (gam_b - 1.0e00),  # 2.50(Not 1.40)
            "fluid_pp(2)%pi_inf": 0,
            "fluid_pp(2)%Re(1)": 1000000,
        }
    )
)

Checklist

  • I have added comments for the new code
  • I added Doxygen docstrings to the new code
  • I have made corresponding changes to the documentation (docs/)
  • I have added regression tests to the test suite so that people can verify in the future that the feature is behaving as expected
  • I have added example cases in examples/ that demonstrate my new feature performing as expected.
    They run to completion and demonstrate "interesting physics"
  • 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)
  • Enclosed the new feature via nvtx ranges so that they can be identified in profiles
  • Ran a Nsight Systems profile using ./mfc.sh run XXXX --gpu -t simulation --nsys, and have attached the output file (.nsys-rep) and plain text results to this PR
  • Ran a Rocprof Systems profile using ./mfc.sh run XXXX --gpu -t simulation --rsys --hip-trace, and have attached the output file and plain text results to this PR.
  • Ran my code using various numbers of different GPUs (1, 2, and 8, for example) in parallel and made sure that the results scale similarly to what happens if you run without the new code/feature

PR Type

Enhancement

Description

  • Add ellipse immersed boundary patch geometry with levelset computation

  • Implement centroid offset calculation for complex IB shapes (airfoils, STL)

  • Support analytical parameterization of moving IB velocities and rotation rates

  • Enable 1-way coupling with analytically-defined IB motion in time-stepping

Diagram Walkthrough

flowchart LR
  A["IB Patch Definition"] -->|"geometry=6"| B["Ellipse Levelset"]
  A -->|"centroid_offset"| C["Centroid Computation"]
  C -->|"airfoil/STL"| D["Center of Mass Calculation"]
  D -->|"offset vector"| E["Levelset Adjustment"]
  F["Analytical Parameters"] -->|"vel, angular_vel"| G["1-Way Coupling"]
  G -->|"time-stepping"| H["IB Motion Update"]
Loading

File Walkthrough

Relevant files
Enhancement
9 files
m_compute_levelset.fpp
Add ellipse levelset computation and centroid offset         
+60/-1   
m_derived_types.fpp
Add centroid offset field to IB patch type                             
+1/-0     
m_ib_patches.fpp
Add ellipse IB patch initialization and marker setup         
+47/-0   
m_check_ib_patches.fpp
Add validation for ellipse IB patch geometry                         
+2/-0     
m_global_parameters.fpp
Initialize centroid offset to zero for all patches             
+1/-0     
m_ibm.fpp
Add centroid offset computation and 1-way coupling support
+56/-1   
m_time_steppers.fpp
Integrate analytical velocity parameterization for 1-way coupling
+5/-1     
case.py
Add code generation for analytical moving IB velocities   
+91/-3   
case_dicts.py
Mark IB velocity and angular velocity as analytically-parameterizable
+2/-4     
Documentation
1 files
case.py
Add example case configuration for ellipse IB simulation 
+97/-0   

CodeAnt-AI Description

Analytic moving immersed-boundary velocities and correct airfoil centroid handling

What Changed

  • Input can now provide analytic formulas for an IB patch's translational and rotational velocities so one-way-coupled bodies move according to user expressions (time- and space-dependent).
  • Airfoil/STL-like patch centroids are computed from the actual immersed-marker distribution and the code applies a centroid offset so geometry is positioned and rotated about its true center of mass.
  • Generated case files now accept and emit analytic moving-IB definitions and will ignore those parameters in numeric input where appropriate.

Impact

✅ Accurate pitching/rotating airfoils about their true center
✅ Can run prescribed analytic IB motions (e.g., oscillating/prescribed trajectories)
✅ Fewer setup mismatches for STL/airfoil geometries

💡 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

Release Notes

  • New Features
    • Added support for ellipse-shaped immersed boundaries in simulations.
    • Enhanced moving immersed boundary capabilities with analytical solution pathways for improved simulation flexibility.
    • Improved immersed boundary positioning with centroid offset adjustments for better accuracy with complex geometries.

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

@codeant-ai
Copy link
Contributor

codeant-ai bot commented Jan 15, 2026

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 Jan 15, 2026

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

Walkthrough

Adds centroid offset field to track true center of mass for immersed boundaries; implements ellipse-based levelset computation; enables parametrized one-way MIBM motion via analytical expressions; fixes airfoil rotations to apply about true centroid instead of tip.

Changes

Cohort / File(s) Summary
Data Type Extensions
src/common/m_derived_types.fpp
Added centroid_offset(1:3) field to ib_patch_parameters type to store center-of-mass offset from cell center for odd-shaped immersed boundaries.
Ellipse Levelset Computation
src/common/m_compute_levelset.fpp
Exposed new s_ellipse_levelset subroutine for ellipse-based level set computation via quadratic form solving; applied centroid offset adjustments to existing s_airfoil_levelset and s_3D_airfoil_levelset paths. Note: contains duplicate subroutine definition.
IB Patch Handling
src/common/m_ib_patches.fpp
Added centroid offset subtraction in s_ib_airfoil and s_ib_3D_airfoil; introduced geometry code 6 dispatcher routing to new internal s_ib_ellipse subroutine for GPU-parallel ellipse membership marking.
IBM Initialization
src/pre_process/m_global_parameters.fpp
Initialize centroid_offset to zero in default immersed boundary patch configuration.
Centroid Computation & Setup
src/simulation/m_ibm.fpp
Added new internal s_compute_centroid_offset subroutine performing center-of-mass aggregation with MPI reductions for geometry codes 4, 5, 11, 12; integrated into s_ibm_setup loop. Note: contains duplicate subroutine insertion.
Time-Stepping Control
src/simulation/m_time_steppers.fpp
Restructured moving IBM conditional logic to distinguish one-way coupling (calls mib_analytical() for analytic velocities) from two-way coupling (retains force/torque computation).
Macro Definitions
src/common/include/case.fpp
Added mib_analytical() macro definition for moving immersed boundary analytical motion.
Toolchain Case Generation
toolchain/mfc/case.py
Introduced MIBM_ANALYTIC_VARS constant; refactored analytic FPP generation into separate __get_analytic_ic_fpp and __get_analytic_mib_fpp methods; added __is_mib_analytical detection; expanded get_inp to recognize MIB analytic parameters.
Patch Parameter Configuration
toolchain/mfc/run/case_dicts.py
Changed patch_ib velocity and angular velocity parameters to use ParamType.REAL.analytic(); removed velocity/angular_vel generation in PRE_PROCESS section.

Estimated code review effort

🎯 4 (Complex) | ⏱️ ~60 minutes

Suggested labels

size:XL

Suggested reviewers

  • sbryngelson
  • wilfonba

Poem

🐰 A rabbit hops through centroid dreams,
Where ellipses spin and IB streams,
No longer tips that twist and turn,
But true centers where physics learn,
One-way waltz with analytic grace! 🌀

🚥 Pre-merge checks | ✅ 1 | ❌ 4
❌ Failed checks (2 warnings, 2 inconclusive)
Check name Status Explanation Resolution
Docstring Coverage ⚠️ Warning Docstring coverage is 20.00% which is insufficient. The required threshold is 80.00%. Write docstrings for the functions missing them to satisfy the coverage threshold.
Description check ⚠️ Warning PR description is comprehensive but significantly incomplete against the template. Most checklist items are unchecked despite reported testing. Check completed items in the checklist (GPU testing, DRY principle, code comments/docstrings). Confirm documentation updates in docs/ and example cases in examples/. Verify ./mfc.sh format was run.
Linked Issues check ❓ Inconclusive While the PR implementation addresses core requirements from #1105 and #1106 (analytic velocity parameterization and centroid offset computation), the summary indicates duplicate implementations of s_ellipse_levelset exist, and incomplete test/GPU validation checklist items raise questions about full requirement satisfaction. Verify that duplicate s_ellipse_levelset definitions are intentional or remove redundancy; confirm all GPU, test, and documentation requirements from the PR checklist are met before merge.
Out of Scope Changes check ❓ Inconclusive Most changes align with requirements (ellipse IB support, centroid offset, analytic motion), but the addition of ellipse geometry (geometry code 6) extends beyond the airfoil/STL centroid fix and introduces new IB shape support not explicitly mentioned in the original issue requirements. Clarify whether ellipse geometry addition is an intentional enhancement or incidental to the core requirements; document rationale if intentional scope expansion.
✅ Passed checks (1 passed)
Check name Status Explanation
Title check ✅ Passed The title 'Analytic mibm velocities and airfoil centroid' directly summarizes the two main changes: support for analytical parameterization of moving IB velocities and correction of airfoil centroid positioning.

✏️ Tip: You can configure your own custom pre-merge checks in the settings.

✨ Finishing touches
  • 📝 Generate docstrings

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.

@codeant-ai codeant-ai bot added the size:L This PR changes 100-499 lines, ignoring generated files label Jan 15, 2026
@qodo-code-review
Copy link
Contributor

PR Reviewer Guide 🔍

Here are some key observations to aid the review process:

🎫 Ticket compliance analysis 🔶

1106 - Partially compliant

Compliant requirements:

  • Fix the "centroid" used for airfoils so it is the true centroid/center of mass rather than the airfoil tip.
  • Ensure rotations in 1-way coupling occur about the center of mass.

Non-compliant requirements:

  • Ensure torque calculations in 2-way coupling use the correct centroid/center of mass reference.

Requires further human verification:

  • Works for 2D and 3D airfoils across all environments.
  • Ensure torque calculations in 2-way coupling use the correct centroid/center of mass reference.

1105 - Partially compliant

Compliant requirements:

  • Support parametrizing 1-way MIBM motion (velocity and rotation rate) as a function of time in the case file.
  • Provide a toolchain/codegen addition that injects analytic functions into a case file (similar to analytic patches/ICs).
  • Apply the prescribed motion inside the time-step update logic using a conditional statement.

Non-compliant requirements:

  • ""

Requires further human verification:

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

Possible Issue

s_compute_centroid_offset appears to be called during s_ibm_setup before ib_markers%sf is allocated/populated (and before the IB patch marker routines run), yet it iterates over ib_markers%sf(i,j,k) to find cells. This risks invalid memory access or always computing a zero-count COM depending on initialization order.

! do all set up for moving immersed boundaries
moving_immersed_boundary_flag = .false.
do i = 1, num_ibs
    if (patch_ib(i)%moving_ibm /= 0) then
        moving_immersed_boundary_flag = .true.
    end if
    call s_update_ib_rotation_matrix(i)
    call s_compute_centroid_offset(i)
end do
$:GPU_ENTER_DATA(copyin='[patch_ib]')

! Allocating the patch identities bookkeeping variable
allocate (patch_id_fp(0:m, 0:n, 0:p))

$:GPU_UPDATE(device='[ib_markers%sf]')
$:GPU_UPDATE(device='[levelset%sf]')
$:GPU_UPDATE(device='[levelset_norm%sf]')

! Get neighboring IB variables from other processors
call s_populate_ib_buffers()

$:GPU_UPDATE(host='[ib_markers%sf]')

! find the number of ghost points and set them to be the maximum total across ranks
call s_find_num_ghost_points(num_gps, num_inner_gps)
call s_mpi_allreduce_integer_sum(num_gps, max_num_gps)
call s_mpi_allreduce_integer_sum(num_inner_gps, max_num_inner_gps)

! set the size of the ghost point arrays to be the amount of points total, plus a factor of 2 buffer
$:GPU_UPDATE(device='[num_gps, num_inner_gps]')
@:ALLOCATE(ghost_points(1:int((max_num_gps + max_num_inner_gps) * 2.0)))
@:ALLOCATE(inner_points(1:int((max_num_gps + max_num_inner_gps) * 2.0)))

$:GPU_ENTER_DATA(copyin='[ghost_points,inner_points]')

call s_find_ghost_points(ghost_points, inner_points)
$:GPU_UPDATE(device='[ghost_points, inner_points]')

call s_compute_image_points(ghost_points, levelset, levelset_norm)
$:GPU_UPDATE(device='[ghost_points]')
Possible Issue

The MPI reduction logic computes global sums into center_of_mass and global cell count into num_cells, but then divides by num_cells_local (local count) instead of the global count. This will yield incorrect centroids in parallel runs. Also verify 2D handling: the accumulation uses [x_cc(i), y_cc(j), 0] but still loops over k; in 2D this may overweight if p spans more than a single plane.

!> Computes the center of mass for IB patch types where we are unable to determine their center of mass analytically.
!> These patches include things like NACA airfoils and STL models
subroutine s_compute_centroid_offset(ib_marker)

    integer, intent(in) :: ib_marker

    integer :: i, j, k, num_cells, num_cells_local
    real(wp), dimension(1:3) :: center_of_mass

    ! Offset only needs to be computes for specific geometries
    if (patch_ib(ib_marker)%geometry == 4 .or. &
        patch_ib(ib_marker)%geometry == 5 .or. &
        patch_ib(ib_marker)%geometry == 11 .or. &
        patch_ib(ib_marker)%geometry == 12) then

        center_of_mass = [0._wp, 0._wp, 0._wp]
        num_cells_local = 0

        ! get the summed mass distribution and number of cells to divide by
        do i = 0, m
            do j = 0, n
                do k = 0, p
                    if (ib_markers%sf(i, j, k) == ib_marker) then
                      num_cells_local = num_cells_local + 1
                        center_of_mass = center_of_mass + [x_cc(i), y_cc(j), 0._wp]
                        if (num_dims == 3) center_of_mass(3) = center_of_mass(3) + z_cc(k)
                    end if
                end do
            end do
        end do

        ! reduce the mass contribution over all MPI ranks and compute COM
        ! print *, "Before reduction ", center_of_mass, num_cells_local
        call s_mpi_allreduce_sum(center_of_mass(1), center_of_mass(1))
        call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2))
        call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3))
        call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
        center_of_mass =  center_of_mass / real(num_cells_local, wp)

        ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid
        patch_ib(ib_marker)%centroid_offset = [patch_ib(ib_marker)%x_centroid, patch_ib(ib_marker)%y_centroid, patch_ib(ib_marker)%z_centroid] &
            - center_of_mass
        patch_ib(ib_marker)%x_centroid = center_of_mass(1)
        patch_ib(ib_marker)%y_centroid = center_of_mass(2)
        patch_ib(ib_marker)%z_centroid = center_of_mass(3)

        ! rotate the centroid offset back into the local coords of the IB
        patch_ib(ib_marker)%centroid_offset = matmul(patch_ib(ib_marker)%rotation_matrix_inverse, patch_ib(ib_marker)%centroid_offset)
    else
        patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp]
    end if
Possible Issue

s_ellipse_levelset only assigns levelset/levelset_norm when a point is inside the ellipse; outside points are left untouched, which can leave uninitialized/stale values in the levelset fields and break downstream logic that assumes the full field is defined. Consider explicitly setting an outside value and normal (or at least levelset) for all cells.

subroutine s_ellipse_levelset(ib_patch_id, levelset, levelset_norm)

    type(levelset_field), intent(INOUT), optional :: levelset
    type(levelset_norm_field), intent(INOUT), optional :: levelset_norm

    integer, intent(in) :: ib_patch_id
    real(wp) :: ellipse_coeffs(2) ! a and b in the ellipse equation
    real(wp) :: quadratic_coeffs(3) ! A, B, C in the quadratic equation to compute levelset

    real(wp) :: length_x, length_y
    real(wp), dimension(1:3) :: xy_local, normal_vector !< x and y coordinates in local IB frame
    real(wp), dimension(2) :: center !< x and y coordinates in local IB frame
    real(wp), dimension(1:3, 1:3) :: rotation, inverse_rotation

    integer :: i, j, k !< Loop index variables
    integer :: idx !< Shortest path direction indicator

    length_x = patch_ib(ib_patch_id)%length_x
    length_y = patch_ib(ib_patch_id)%length_y
    center(1) = patch_ib(ib_patch_id)%x_centroid
    center(2) = patch_ib(ib_patch_id)%y_centroid
    inverse_rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix_inverse(:, :)
    rotation(:, :) = patch_ib(ib_patch_id)%rotation_matrix(:, :)

    ellipse_coeffs(1) = 0.5_wp * length_x
    ellipse_coeffs(2) = 0.5_wp * length_y

    $:GPU_PARALLEL_LOOP(private='[i,j,k,idx,quadratic_coeffs,xy_local,normal_vector]', &
              & copyin='[ib_patch_id,center,ellipse_coeffs,inverse_rotation,rotation]', collapse=2)
    do i = 0, m
        do j = 0, n
            xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp]
            xy_local = matmul(inverse_rotation, xy_local)

            ! we will get NaNs in the levelset if we compute this outside the ellipse
            if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then

                normal_vector = xy_local
                normal_vector(2) = normal_vector(2) * (ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transofmration method
                normal_vector = normal_vector / sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
                levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame

                ! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3
                quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2
                quadratic_coeffs(2) = 2._wp * ((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2)))
                quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp

                ! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
                levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1)
            end if
        end do
    end do
    $:END_GPU_PARALLEL_LOOP()

end subroutine s_ellipse_levelset

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 toolchain code for analytic IB motion incorrectly allows spatial variables (x, y, z) in expressions for velocity and angular velocity. This should be refactored to only permit time (t) as a variable, as rigid body motion is not spatially dependent. [High-level, importance: 9]

Solution Walkthrough:

Before:

# toolchain/mfc/case.py
def __get_analytic_mib_fpp(self, print: bool) -> str:
    # ...
    def rhs_replace(match):
        return {
            'x': 'x_cc(i)', 'y': 'y_cc(j)', 'z': 'z_cc(k)',
            't': 'mytime',
            'r': f'patch_ib({pid})%radius',
            'e' : f'{math.e}',
        }.get(match.group(), match.group())
    # ...
    for attribute, expr in items:
        rhs = re.sub(r"[a-zA-Z]+", rhs_replace, expr)
        lines.append(f"        {attribute} = {rhs}")
    # ...

After:

# toolchain/mfc/case.py
def __get_analytic_mib_fpp(self, print: bool) -> str:
    # ...
    def rhs_replace(match):
        return {
            't': 'mytime',
            'r': f'patch_ib({pid})%radius',
            'e' : f'{math.e}',
        }.get(match.group(), match.group())
    # ...
    for attribute, expr in items:
        rhs = re.sub(r"[a-zA-Z]+", rhs_replace, expr)
        lines.append(f"        {attribute} = {rhs}")
    # ...

call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2))
call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3))
call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
center_of_mass = center_of_mass / real(num_cells_local, wp)
Copy link
Contributor

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

Choose a reason for hiding this comment

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

Suggestion: In s_compute_centroid_offset, divide center_of_mass by the global cell count num_cells instead of the local count num_cells_local to correctly calculate the global center of mass after MPI reduction. [possible issue, importance: 9]

Suggested change
center_of_mass = center_of_mass / real(num_cells_local, wp)
center_of_mass = center_of_mass / real(num_cells, wp)

Comment on lines 376 to 391
! we will get NaNs in the levelset if we compute this outside the ellipse
if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then

normal_vector = xy_local
normal_vector(2) = normal_vector(2) * (ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transofmration method
normal_vector = normal_vector / sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame

! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3
quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2
quadratic_coeffs(2) = 2._wp * ((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2)))
quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp

! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1)
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_ellipse_levelset, compute the level set for all points, not just those inside the ellipse, to ensure the signed distance function is correctly defined everywhere. Use an algebraic distance approximation for points outside the ellipse. [possible issue, importance: 8]

Suggested change
! we will get NaNs in the levelset if we compute this outside the ellipse
if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
normal_vector = xy_local
normal_vector(2) = normal_vector(2) * (ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transofmration method
normal_vector = normal_vector / sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame
! use the normal vector to set up the quadratic equation for the levelset, using A, B, and C in indices 1, 2, and 3
quadratic_coeffs(1) = (normal_vector(1)/ellipse_coeffs(1))**2 + (normal_vector(2)/ellipse_coeffs(2))**2
quadratic_coeffs(2) = 2._wp * ((xy_local(1)*normal_vector(1)/(ellipse_coeffs(1)**2)) + (xy_local(2)*normal_vector(2)/(ellipse_coeffs(2)**2)))
quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp
! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1)
end if
! algebraic distance approximation
levelset%sf(i, j, 0, ib_patch_id) = sqrt(ellipse_coeffs(1) * ellipse_coeffs(2)) * &
(sqrt((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2) - 1._wp)
! we still need to compute the normal vector for the ghost cell method
if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
normal_vector = xy_local
normal_vector(2) = normal_vector(2) * (ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp ! get the normal direction via the coordinate transofmration method
normal_vector = normal_vector / sqrt(dot_product(normal_vector, normal_vector)) ! normalize the vector
levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector) ! save after rotating the vector to the global frame
end if

quadratic_coeffs(3) = (xy_local(1)/ellipse_coeffs(1))**2._wp + (xy_local(2)/ellipse_coeffs(2))**2._wp - 1._wp

! compute the levelset with the quadratic equation [ -B + sqrt(B^2 - 4AC) ] / 2A
levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: The quadratic root is computed with the wrong sign and denominator: the existing expression evaluates to (B - sqrt(...)) / (2A) instead of the correct (-B + sqrt(...)) / (2A), producing a levelset with the wrong sign; change the expression to use the standard quadratic formula denominator (2*A) with the correct numerator sign. [logic error]

Severity Level: Critical 🚨
- ❌ s_ellipse_levelset produces wrong ellipse geometry.
- ⚠️ Simulations using ellipse IB patches misrepresent surfaces.
- ⚠️ Postprocessing (levelset-dependent diagnostics) inaccurate.
Suggested change
levelset%sf(i, j, 0, ib_patch_id) = -0.5_wp * (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / quadratic_coeffs(1)
levelset%sf(i, j, 0, ib_patch_id) = (-quadratic_coeffs(2) + sqrt(quadratic_coeffs(2)**2._wp - 4._wp * quadratic_coeffs(1)*quadratic_coeffs(3))) / (2._wp * quadratic_coeffs(1))
Steps of Reproduction ✅
1. Compile the PR code so module m_compute_levelset is built (file:
src/common/m_compute_levelset.fpp). This produces the subroutine s_ellipse_levelset
defined starting at line ~342.

2. Arrange to call s_ellipse_levelset(ib_patch_id, levelset, levelset_norm) directly (unit
test or existing levelset driver that dispatches ellipse patches). When a patch_ib entry
has non-zero length_x and length_y the routine runs the GPU loop (lines 371-393) and
evaluates the ellipse branch at lines ~377-391.

3. For any grid cell where the condition ((xy_local(1)/ellipse_coeffs(1))**2 +
(xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) is true (see lines ~377-388), the code
computes quadratic_coeffs and executes the assignment at line 390 (the existing_code
line). This is the concrete execution point for the reported issue.

4. Observe resulting levelset values: compare levelset%sf(i,j,0,ib_patch_id) produced at
line 390 against the analytic nearest-surface distance computed via (-B +
sqrt(B^2-4AC))/(2A). The current expression at line 390 simplifies to (B - sqrt(...))/(2A)
(wrong sign/ordering), producing an incorrect sign/branch for the root. Replacing the
expression with the improved code yields the expected analytic root.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/m_compute_levelset.fpp
**Line:** 390:390
**Comment:**
	*Logic Error: The quadratic root is computed with the wrong sign and denominator: the existing expression evaluates to (B - sqrt(...)) / (2A) instead of the correct (-B + sqrt(...)) / (2A), producing a levelset with the wrong sign; change the expression to use the standard quadratic formula denominator (2*A) with the correct numerator sign.

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.

do i = 0, m
xy_local = [x_cc(i) - center(1), y_cc(j) - center(2), 0._wp] ! get coordinate frame centered on IB
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinates
xy_local = xy_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: The centroid offset is subtracted from xy_local after rotating the global coordinate into the IB frame, but patch_ib(patch_id)%centroid_offset is likely defined in the global/model frame and must be rotated into the same frame before subtraction — otherwise the offset is applied in the wrong coordinate system causing misplaced airfoil geometry. Rotate the offset into the IB/local frame (or subtract the offset before rotation) so both vectors are in the same coordinate basis. [logic error]

Severity Level: Critical 🚨
- ❌ 2D airfoil placement in s_ib_airfoil mispositioned.
- ⚠️ Pitching-airfoil tests produce incorrect geometry mapping.
Suggested change
xy_local = xy_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
xy_local = xy_local - matmul(inverse_rotation, patch_ib(patch_id)%centroid_offset) ! rotate centroid_offset into IB/local frame before subtracting
Steps of Reproduction ✅
1. Prepare a 2D case that includes an airfoil patch (patch_ib(... )%geometry == 4). The 2D
branch in s_apply_ib_patches calls s_ib_airfoil (src/common/m_ib_patches.fpp — the 2D IB
patch logic that eventually invokes s_ib_airfoil).

2. Ensure the patch has a non-zero centroid_offset and a non-trivial rotation (set
patch_ib(patch_id)%centroid_offset and non-zero patch_ib(patch_id)%angles so
rotation_matrix_inverse != identity). Run to the phase where IB markers are populated so
s_ib_airfoil executes.

3. Inside s_ib_airfoil, execution computes xy_local and rotates it at the lines shown in
the diff (see src/common/m_ib_patches.fpp:279-281 where xy_local is computed and rotated,
then the subtraction at line 282). Because the code subtracts
patch_ib(patch_id)%centroid_offset without rotating it, the offset is applied in the
already-rotated/local frame incorrectly.

4. Observe incorrect IB marking: cell membership checks that follow (the if on xy_local(1)
at src/common/m_ib_patches.fpp:284 and subsequent assignments to ib_markers_sf) will mark
the airfoil in the wrong location or shape. This reproduces reliably when centroid_offset
is set and rotation is non-zero; the fix (rotate centroid_offset into local frame) aligns
coordinate bases and corrects marker placement.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/m_ib_patches.fpp
**Line:** 282:282
**Comment:**
	*Logic Error: The centroid offset is subtracted from `xy_local` after rotating the global coordinate into the IB frame, but `patch_ib(patch_id)%centroid_offset` is likely defined in the global/model frame and must be rotated into the same frame before subtraction — otherwise the offset is applied in the wrong coordinate system causing misplaced airfoil geometry. Rotate the offset into the IB/local frame (or subtract the offset before rotation) so both vectors are in the same coordinate basis.

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.

do i = 0, m
xyz_local = [x_cc(i) - center(1), y_cc(j) - center(2), z_cc(l) - center(3)] ! get coordinate frame centered on IB
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordinates
xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: The 3D airfoil routine subtracts patch_ib(patch_id)%centroid_offset from xyz_local after applying inverse_rotation; if centroid_offset is expressed in the global/model frame this subtraction is done in the wrong basis. Rotate the centroid offset into the IB/local frame with the same inverse_rotation before subtracting to ensure the offset is applied consistently in the same coordinate system. [logic error]

Severity Level: Critical 🚨
- ❌ 3D airfoil marker placement incorrect in s_ib_3D_airfoil.
- ⚠️ 3D pitching/rotating airfoil simulations show wrong geometry.
Suggested change
xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
xyz_local = xyz_local - matmul(inverse_rotation, patch_ib(patch_id)%centroid_offset) ! rotate centroid_offset into IB/local frame before subtracting
Steps of Reproduction ✅
1. Prepare a 3D case with an airfoil patch (patch_ib(... )%geometry == 11) so
s_apply_ib_patches calls s_ib_3D_airfoil during marker generation.

2. Set patch_ib(patch_id)%centroid_offset to a non-zero vector and set non-zero rotation
angles so rotation_matrix_inverse is non-trivial; run to the IB marker generation step.

3. In s_ib_3D_airfoil the code computes xyz_local and applies inverse_rotation at
src/common/m_ib_patches.fpp:435-436, then subtracts patch_ib(patch_id)%centroid_offset at
line 437 without rotating that offset into the local frame.

4. As a result the z/y/x cell-inclusion tests (for example the z_min/z_max check at the
surrounding block) and the later assignments to ib_markers_sf(i,j,l) are done with an
inconsistent offset and the 3D airfoil appears shifted. Reproduces when centroid_offset ≠
0 and rotation ≠ identity.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/m_ib_patches.fpp
**Line:** 437:437
**Comment:**
	*Logic Error: The 3D airfoil routine subtracts `patch_ib(patch_id)%centroid_offset` from `xyz_local` after applying `inverse_rotation`; if `centroid_offset` is expressed in the global/model frame this subtraction is done in the wrong basis. Rotate the centroid offset into the IB/local frame with the same `inverse_rotation` before subtracting to ensure the offset is applied consistently in the same coordinate system.

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.

patch_ib(i)%geometry == 12) then
call s_check_model_ib_patch_geometry(i)
else if (patch_ib(i)%geometry == 6) then
print *, "Ellipse Patch"
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Missing validation hook: if ellipse geometry (6) is intended to be supported, the branch should invoke a dedicated validation routine (like other geometries do) rather than just printing; call the ellipse check subroutine so ellipse parameters are validated consistently with other patch types. [possible bug]

Severity Level: Critical 🚨
- ❌ Build fails at m_check_ib_patches.fpp:65 blocking compilation.
- ❌ CI and developer builds will be broken until function implemented.
- ⚠️ Prevents running any preprocessing or simulations.
Suggested change
print *, "Ellipse Patch"
call s_check_ellipse_ib_patch_geometry(i)
Steps of Reproduction ✅
1. Apply the proposed change that replaces the print with a call to

   s_check_ellipse_ib_patch_geometry in src/pre_process/m_check_ib_patches.fpp

   at lines 65-66 (i.e., insert the call in place of the print).

2. Run the project's build (compile all Fortran sources; e.g., invoke the project's

   usual build target that compiles src/pre_process/m_check_ib_patches.fpp).

3. The Fortran compiler / linker will fail because there is no definition of

   s_check_ellipse_ib_patch_geometry in the current codebase (a search of

   src/pre_process/m_check_ib_patches.fpp and the file's subroutines shows no such

   symbol). The error will reference the call site at

   src/pre_process/m_check_ib_patches.fpp:65 (undefined reference / missing procedure).

4. Result: build fails immediately; unless s_check_ellipse_ib_patch_geometry is added

   elsewhere, this change is a compile-time break, not a runtime validation improvement.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/pre_process/m_check_ib_patches.fpp
**Line:** 66:66
**Comment:**
	*Possible Bug: Missing validation hook: if ellipse geometry (6) is intended to be supported, the branch should invoke a dedicated validation routine (like other geometries do) rather than just printing; call the ellipse check subroutine so ellipse parameters are validated consistently with other patch types.

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.

call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2))
call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3))
call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
center_of_mass = center_of_mass / real(num_cells_local, wp)
Copy link
Contributor

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

Choose a reason for hiding this comment

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

Suggestion: Logic error: the code divides the reduced (global) sum of coordinates by the local rank cell count (num_cells_local) instead of the global cell count (num_cells), which yields an incorrect centroid and can cause a division-by-zero when the local rank has zero cells; change the divisor to the globally reduced num_cells and guard against zero. [logic error]

Severity Level: Critical 🚨
- ❌ Centroid miscomputed for airfoil/STL IBs (s_compute_centroid_offset)
- ❌ Possible runtime crash on ranks with zero local IB cells
- ⚠️ Affects s_ibm_setup and subsequent IB motion calculations
- ⚠️ Breaks pitching-airfoil test cases (user-reported)
Suggested change
center_of_mass = center_of_mass / real(num_cells_local, wp)
if (num_cells == 0) then
! No cells globally for this IB -- leave offsets at zero and exit safely
patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp]
return
end if
center_of_mass = center_of_mass / real(num_cells, wp)
Steps of Reproduction ✅
1. Build and run the code with MPI (>=2 ranks) and a case that registers an IB of geometry
type handled here (geometry == 4/5/11/12). The initialization path runs s_ibm_setup
(src/simulation/m_ibm.fpp: around lines 93-105) which calls s_compute_centroid_offset.

2. Ensure the IB lies entirely within one MPI rank's subdomain such that other ranks have
zero cells for that IB. During s_compute_centroid_offset (subroutine start at
src/simulation/m_ibm.fpp: ~1086) each rank counts local cells into num_cells_local and
sums coordinates into center_of_mass.

3. The code calls s_mpi_allreduce_integer_sum(num_cells_local, num_cells) (line 1120)
which returns global num_cells, but then divides by the local num_cells_local at line
1121. If the local rank had num_cells_local == 0, this division triggers a divide-by-zero
on that rank causing a crash; otherwise the centroid is computed incorrectly (using local
count instead of global).

4. Observe either a floating-point exception / program abort on ranks with zero local
cells or an incorrect centroid value being written to patch_ib(...)%x_centroid etc.
Confirm by inspecting outputs written after s_compute_centroid_offset or by adding a debug
print near lines 1115-1121.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_ibm.fpp
**Line:** 1121:1121
**Comment:**
	*Logic Error: Logic error: the code divides the reduced (global) sum of coordinates by the local rank cell count (`num_cells_local`) instead of the global cell count (`num_cells`), which yields an incorrect centroid and can cause a division-by-zero when the local rank has zero cells; change the divisor to the globally reduced `num_cells` and guard against zero.

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.

Comment on lines +369 to +371
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic()
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Type mismatch: assigning the result of ParamType.REAL.analytic() (a plain dict) directly into SIMULATION mixes raw dicts with ParamType Enum members. Later the code builds _properties = {k: v.value for k, v in ALL.items()} which assumes every value is an Enum with a .value attribute; a raw dict will raise AttributeError. Replace the two analytic assignments with the Enum member ParamType.REAL to keep dictionary values consistent and avoid the runtime AttributeError. [type error]

Severity Level: Critical 🚨
- ❌ Module import fails building schema and validator.
- ❌ get_validator() unusable; schema compilation blocked.
- ⚠️ Any tool consuming CASE options fails initialization.
Suggested change
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic()
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
Steps of Reproduction ✅
1. Open the module file toolchain/mfc/run/case_dicts.py and locate the
ParamType.analytic() implementation (function at toolchain/mfc/run/case_dicts.py:17-22).
The analytic() returns a plain dict for ParamType.REAL.

2. Observe the SIMULATION assignments inside the ib loop at the hunk lines where the PR
added analytic variants (see toolchain/mfc/run/case_dicts.py lines 369 and 371): those two
lines assign the raw dicts returned from ParamType.REAL.analytic() into SIMULATION.

3. Import the module (e.g., python -c "import toolchain.mfc.run.case_dicts as c") or call
get_validator() from toolchain/mfc/run/case_dicts.py. Module import triggers building of
ALL and then the comprehension _properties = { k: v.value for k, v in ALL.items() } (the
_properties line is defined immediately after ALL is assembled in the same file).

4. At that comprehension, Python evaluates v.value for every value in ALL. For the two
entries assigned above, v is a plain dict and has no .value attribute; this raises
AttributeError: 'dict' object has no attribute 'value'. The import or get_validator() call
fails with that exception.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** toolchain/mfc/run/case_dicts.py
**Line:** 369:371
**Comment:**
	*Type Error: Type mismatch: assigning the result of `ParamType.REAL.analytic()` (a plain dict) directly into `SIMULATION` mixes raw dicts with `ParamType` Enum members. Later the code builds `_properties = {k: v.value for k, v in ALL.items()}` which assumes every value is an Enum with a `.value` attribute; a raw dict will raise AttributeError. Replace the two analytic assignments with the Enum member `ParamType.REAL` to keep dictionary values consistent and avoid the runtime AttributeError.

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.

@codeant-ai
Copy link
Contributor

codeant-ai bot commented Jan 15, 2026

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

Caution

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

⚠️ Outside diff range comments (2)
src/pre_process/m_check_ib_patches.fpp (1)

65-70: Add validation for ellipse geometry and update error message to include geometry 6.

The ellipse patch (geometry == 6) branch lacks parameter validation that exists for all other geometry types. Additionally, the error message on line 70 must be updated to include geometry 6 in the list of valid geometries to match the explicit handling at line 65.

🔧 Suggested fixes
  1. Update the error message to include geometry 6:
                 call s_prohibit_abort("Invalid IB patch", &
                                       "patch_ib("//trim(iStr)//")%geometry must be "// &
-                                      "2-4, 8-10, 11 or 12.")
+                                      "2-6, 8-10, 11 or 12.")
  1. Add a validation subroutine for ellipse geometry and call it:
else if (patch_ib(i)%geometry == 6) then
    call s_check_ellipse_ib_patch_geometry(i)

Then implement s_check_ellipse_ib_patch_geometry to validate centroid and radii parameters similarly to the circle patch validation.

src/simulation/m_time_steppers.fpp (1)

629-643: The @:mib_analytical() macro at line 631 is undefined and will cause a preprocessing error.

The macro is invoked for 1-way coupling (moving_ibm == 1) but is not defined in case.fpp, macros.fpp, or any other .fpp file in the codebase. Only an empty analytical() macro (without the mib_ prefix) exists in case.fpp.

The macro must be properly defined to specify how analytic velocities are applied to the IB patches in 1-way coupling mode. Add the implementation to case.fpp or the appropriate macro file.

🤖 Fix all issues with AI agents
In `@examples/2D_ibm_ellipse/case.py`:
- Around line 84-90: Update the misleading comment "Cylinder Immersed Boundary"
to describe the actual geometry being configured (an ellipse) where the patch
keys like "patch_ib(1)%geometry" are set to 6; change the comment near the block
that sets "patch_ib(1)%geometry", "patch_ib(1)%x_centroid",
"patch_ib(1)%y_centroid", "patch_ib(1)%length_x", "patch_ib(1)%length_y", and
"patch_ib(1)%slip" so it accurately states "Ellipse Immersed Boundary" (or
equivalent) to match geometry=6.

In `@src/common/m_compute_levelset.fpp`:
- Line 380: The inline comment on the line assigning normal_vector(2) contains a
typo "transofmration"; update the comment text to "transformation" in the same
statement so it reads "... via the coordinate transformation method" while
leaving the code (normal_vector(2) = normal_vector(2) *
(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp) and surrounding whitespace
unchanged.

In `@src/common/m_ib_patches.fpp`:
- Around line 789-792: The comment in the m_ib_patches module that reads
"Checking whether the rectangle covers..." is incorrect for the ellipse-handling
subroutine; update the comment to say "ellipse" (and adjust surrounding phrasing
if needed, e.g., "ellipse covers a particular cell") so it accurately reflects
the geometry handled by the subroutine in this file.

In `@src/simulation/m_ibm.fpp`:
- Around line 1115-1122: The center_of_mass calculation divides by the local
count; after calling s_mpi_allreduce_sum(center_of_mass) and
s_mpi_allreduce_integer_sum(num_cells_local, num_cells) replace the division by
real(num_cells_local, wp) with real(num_cells, wp) so the globally reduced
center_of_mass is normalized by the global cell count (ensure you use the
reduced integer num_cells, not num_cells_local, in the division and keep the wp
real-kind conversion).
🧹 Nitpick comments (1)
toolchain/mfc/case.py (1)

223-256: Consider binding pid explicitly to avoid closure capture issue.

The static analyzer flagged that rhs_replace captures the loop variable pid by reference. While this is safe here since the function is called immediately within the same iteration, it's a potential maintenance hazard if the code is refactored.

♻️ Optional fix using default argument binding
         for pid, items in ib_patches.items():
 
             # function that defines how we will replace variable names with
             # values from the case file
-            def rhs_replace(match):
+            def rhs_replace(match, pid=pid):
                 return {
                     'x': 'x_cc(i)', 'y': 'y_cc(j)', 'z': 'z_cc(k)',
                     't': 'mytime',
 
                     'r':     f'patch_ib({pid})%radius',
 
                     'e' : f'{math.e}',
                 }.get(match.group(), match.group())
📜 Review details

Configuration used: defaults

Review profile: CHILL

Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between 31c938f and 7bd28cd.

📒 Files selected for processing (10)
  • examples/2D_ibm_ellipse/case.py
  • src/common/m_compute_levelset.fpp
  • src/common/m_derived_types.fpp
  • src/common/m_ib_patches.fpp
  • src/pre_process/m_check_ib_patches.fpp
  • src/pre_process/m_global_parameters.fpp
  • src/simulation/m_ibm.fpp
  • src/simulation/m_time_steppers.fpp
  • toolchain/mfc/case.py
  • toolchain/mfc/run/case_dicts.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/common/m_derived_types.fpp
  • src/pre_process/m_check_ib_patches.fpp
  • src/simulation/m_ibm.fpp
  • src/simulation/m_time_steppers.fpp
  • src/common/m_compute_levelset.fpp
  • src/common/m_ib_patches.fpp
  • src/pre_process/m_global_parameters.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/common/m_derived_types.fpp
  • src/pre_process/m_check_ib_patches.fpp
  • src/simulation/m_ibm.fpp
  • src/simulation/m_time_steppers.fpp
  • src/common/m_compute_levelset.fpp
  • src/common/m_ib_patches.fpp
  • src/pre_process/m_global_parameters.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_ibm.fpp
  • src/simulation/m_time_steppers.fpp
🧠 Learnings (7)
📚 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_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_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_ibm.fpp
  • src/simulation/m_time_steppers.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_time_steppers.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_time_steppers.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_time_steppers.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/common/m_compute_levelset.fpp
🧬 Code graph analysis (3)
src/simulation/m_ibm.fpp (1)
toolchain/mfc/test/cases.py (1)
  • alter_ib (377-422)
examples/2D_ibm_ellipse/case.py (1)
toolchain/mfc/test/cases.py (8)
  • alter_ib (377-422)
  • alter_3d (322-365)
  • alter_2d (293-320)
  • alter_hypoelasticity (590-652)
  • alter_num_fluids (224-291)
  • alter_mixlayer_perturb (686-717)
  • alter_viscosity (788-839)
  • alter_bc_patches (908-955)
toolchain/mfc/case.py (1)
toolchain/mfc/common.py (1)
  • is_number (190-201)
🪛 GitHub Actions: Spell Check
src/common/m_compute_levelset.fpp

[error] 380-380: Typos found. 'transofmration' should be 'transformation'.

🪛 Ruff (0.14.11)
toolchain/mfc/case.py

232-232: Function definition does not bind loop variable pid

(B023)

⏰ 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). (11)
  • GitHub Check: Code Cleanliness Check
  • GitHub Check: Coverage Test on CodeCov
  • GitHub Check: Github (ubuntu, no-mpi, single, no-debug, false)
  • GitHub Check: Github (ubuntu, mpi, no-debug, false)
  • GitHub Check: Github (macos, mpi, no-debug, false)
  • GitHub Check: Github (macos, mpi, debug, false)
  • GitHub Check: Github (ubuntu, mpi, no-debug, true)
  • GitHub Check: Github (ubuntu, mpi, debug, false)
  • GitHub Check: Github (ubuntu, mpi, debug, true)
  • GitHub Check: cubic · AI code reviewer
  • GitHub Check: Build & Publish
🔇 Additional comments (20)
toolchain/mfc/run/case_dicts.py (1)

368-371: LGTM! Analytic type support for IB velocity parameters.

The change to ParamType.REAL.analytic() for vel and angular_vel correctly enables time-parameterized motion expressions for moving immersed boundaries, consistent with how other analytic initial conditions are handled elsewhere in the codebase.

toolchain/mfc/case.py (5)

15-18: LGTM! Well-structured constant for MIB analytic variables.

The MIBM_ANALYTIC_VARS list correctly enumerates all velocity and angular velocity components that can be analytically defined for moving immersed boundaries.


101-111: LGTM! Proper regex escaping for variable matching.

Good use of re.escape(variable) to handle parentheses in variable names like vel(1). The method correctly mirrors the pattern established by __is_ic_analytical.


350-358: LGTM! Clean refactoring for FPP generation.

The separation of __get_analytic_ic_fpp and __get_analytic_mib_fpp with proper concatenation in __get_sim_fpp ensures both IC and MIB analytic macros are available during simulation compilation.


54-57: LGTM! Proper handling of MIB analytic expressions in input generation.

Correctly extends the analytic placeholder logic to MIB parameters, ensuring string expressions are replaced with 0d0 in the .inp file while the actual expressions are handled via generated Fypp macros.


258-265: No action needed — mytime is properly available in scope.

The generated #:def mib_analytical() macro correctly maps the time variable t to mytime, which is defined as a module-level variable in m_global_parameters.fpp and is initialized and maintained throughout the simulation lifecycle. It is accessible wherever the macro is expanded.

src/common/m_derived_types.fpp (1)

299-299: LGTM! Well-documented centroid offset field.

The centroid_offset field correctly addresses issue #1106 by allowing proper tracking of the center of mass offset for asymmetric immersed boundaries like airfoils. The field type and dimensioning are appropriate.

src/pre_process/m_global_parameters.fpp (1)

565-565: LGTM! Proper default initialization for centroid_offset.

The initialization to 0._wp is consistent with other array initializations in the IB patch defaults and ensures a clean starting state before the actual centroid offset is computed during simulation setup.

examples/2D_ibm_ellipse/case.py (1)

1-97: LGTM - New 2D ellipse example case.

The case file correctly demonstrates the new ellipse IB patch (geometry=6) with appropriate domain parameters, boundary conditions, and fluid properties for viscous flow simulation.

src/simulation/m_time_steppers.fpp (1)

5-6: LGTM - Include case.fpp for analytic MIB macros.

Adding the include enables the @:mib_analytical() macro for prescribing analytic velocities in 1-way MIBM coupling.

src/simulation/m_ibm.fpp (2)

96-104: LGTM - Added centroid offset computation in IBM setup.

The call to s_compute_centroid_offset(i) is correctly placed after the rotation matrix update and before GPU data transfer, enabling proper centroid calculations for airfoil and STL geometries.


1103-1113: Use the more direct if (p > 0) check for consistency with the rest of the codebase. Line 1109 should be if (p > 0) instead of if (num_dims == 3), matching the pattern used elsewhere (e.g., lines 224–227) and avoiding reliance on a global parameter that should inherently correspond to the grid dimensions.

src/common/m_compute_levelset.fpp (4)

26-27: LGTM - Added ellipse levelset to public exports.

The new s_ellipse_levelset subroutine is correctly added to the module's public interface.


97-97: LGTM - Centroid offset adjustment for 2D airfoil levelset.

Correctly subtracts centroid_offset from local coordinates after rotation, aligning the airfoil's geometric center with its true center of mass for proper levelset computation.


193-193: LGTM - Centroid offset adjustment for 3D airfoil levelset.

Consistent with the 2D implementation, correctly applies centroid offset in the local IB frame.


342-396: Resolve the sign discrepancy between the formula comment and implementation on line 390.

The comment on line 389 states to compute [ -B + sqrt(B^2 - 4AC) ] / 2A (the plus root), but the code on line 390 implements -0.5_wp * (-quadratic_coeffs(2) + sqrt(...)) / quadratic_coeffs(1), which mathematically evaluates to (B - sqrt(...)) / (2A) (the minus root with a sign flip). Verify whether this negation is intentional (e.g., for inward-facing signed distance) and update the comment to reflect the actual formula being computed.

Regarding the other points: centroid_offset is correctly omitted since ellipses have a well-defined geometric center (unlike airfoils which require explicit offset correction, as noted in s_airfoil_levelset). Minor style: replace **2._wp exponents on lines 380, 387 with **2 per coding guidelines.

src/common/m_ib_patches.fpp (4)

122-125: LGTM - Added ellipse geometry dispatch in 2D IB patch application.

The new geometry type 6 for ellipse is correctly routed to both the marker function (s_ib_ellipse) and levelset computation (s_ellipse_levelset).


282-282: LGTM - Centroid offset for 2D airfoil patch markers.

Correctly applies centroid offset in the local IB frame, consistent with the levelset computation in m_compute_levelset.fpp.


437-437: LGTM - Centroid offset for 3D airfoil patch markers.

Consistent implementation with the 2D version.


771-810: LGTM - New ellipse IB patch subroutine.

The implementation correctly:

  • Extracts ellipse semi-axes from length_x/2 and length_y/2
  • Transforms coordinates to local IB frame using rotation matrix
  • Uses the standard ellipse equation (x/a)² + (y/b)² ≤ 1 for interior check
  • Follows the same pattern as other 2D patch subroutines

The GPU parallel loop is properly configured with appropriate private variables and data movement.

✏️ Tip: You can disable this entire section by setting review_details to false in your review settings.

Comment on lines 84 to 90
# Patch: Cylinder Immersed Boundary
"patch_ib(1)%geometry": 6,
"patch_ib(1)%x_centroid": 1.5e-03,
"patch_ib(1)%y_centroid": 3.0e-03,
"patch_ib(1)%length_x": 0.4e-03,
"patch_ib(1)%length_y": 0.2e-03,
"patch_ib(1)%slip": "F",
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 | 🟡 Minor

Update the comment to reflect the actual geometry.

The comment on line 84 says "Cylinder Immersed Boundary" but this patch defines an ellipse (geometry=6). Update the comment to match the actual geometry being configured.

-            # Patch: Cylinder Immersed Boundary
+            # Patch: Ellipse Immersed Boundary
             "patch_ib(1)%geometry": 6,
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
# Patch: Cylinder Immersed Boundary
"patch_ib(1)%geometry": 6,
"patch_ib(1)%x_centroid": 1.5e-03,
"patch_ib(1)%y_centroid": 3.0e-03,
"patch_ib(1)%length_x": 0.4e-03,
"patch_ib(1)%length_y": 0.2e-03,
"patch_ib(1)%slip": "F",
# Patch: Ellipse Immersed Boundary
"patch_ib(1)%geometry": 6,
"patch_ib(1)%x_centroid": 1.5e-03,
"patch_ib(1)%y_centroid": 3.0e-03,
"patch_ib(1)%length_x": 0.4e-03,
"patch_ib(1)%length_y": 0.2e-03,
"patch_ib(1)%slip": "F",
🤖 Prompt for AI Agents
In `@examples/2D_ibm_ellipse/case.py` around lines 84 - 90, Update the misleading
comment "Cylinder Immersed Boundary" to describe the actual geometry being
configured (an ellipse) where the patch keys like "patch_ib(1)%geometry" are set
to 6; change the comment near the block that sets "patch_ib(1)%geometry",
"patch_ib(1)%x_centroid", "patch_ib(1)%y_centroid", "patch_ib(1)%length_x",
"patch_ib(1)%length_y", and "patch_ib(1)%slip" so it accurately states "Ellipse
Immersed Boundary" (or equivalent) to match geometry=6.

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.

6 issues found across 10 files

Prompt for AI agents (all 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:1121">
P1: After the MPI reductions, the centroid is divided by `num_cells_local` (the pre-reduction local count) instead of the globally reduced `num_cells`, which can produce NaNs or incorrect centers of mass. Divide by the global count returned by `s_mpi_allreduce_integer_sum` (and guard against zero) to compute a valid centroid.</violation>
</file>

<file name="src/pre_process/m_check_ib_patches.fpp">

<violation number="1" location="src/pre_process/m_check_ib_patches.fpp:65">
P2: Geometry 6 IB patches bypass all validation: the new branch just prints “Ellipse Patch” and never enforces required parameters (centroid, radii, dimensionality). This lets malformed ellipse IB definitions pass unchecked and can cause runtime failures when s_ib_ellipse later uses unset data. Replace the print with proper validation mirroring other geometries.</violation>
</file>

<file name="toolchain/mfc/run/case_dicts.py">

<violation number="1" location="toolchain/mfc/run/case_dicts.py:369">
P1: Assigning `ParamType.REAL.analytic()` (which returns a plain dict) directly into `SIMULATION` is inconsistent with other entries that use `ParamType` Enum members. Downstream code that builds `_properties = {k: v.value for k, v in ALL.items()}` will fail with an `AttributeError` because dicts don't have a `.value` attribute. Use `ParamType.REAL` instead to maintain consistency.</violation>

<violation number="2" location="toolchain/mfc/run/case_dicts.py:371">
P1: Same issue as with `vel`: assigning `ParamType.REAL.analytic()` (a plain dict) will cause `AttributeError` when building `_properties`. Use `ParamType.REAL` instead.</violation>
</file>

<file name="src/common/m_compute_levelset.fpp">

<violation number="1" location="src/common/m_compute_levelset.fpp:377">
P1: levelset and normal fields are never set for cells outside the ellipse, leaving undefined geometry data for most of the domain. Add an `else` branch that computes the signed distance/normal outside the ellipse instead of skipping the assignment.</violation>
</file>

<file name="toolchain/mfc/case.py">

<violation number="1" location="toolchain/mfc/case.py:227">
P1: Analytic IB codegen replaces `x`, `y`, `z` with grid-center arrays using undefined indices (`j`, `k`) inside the time-stepper loop, causing compilation failures or wrong values when users reference spatial variables.</violation>
</file>

Reply with feedback, questions, or to request a fix. Tag @cubic-dev-ai to re-run a review.

call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2))
call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3))
call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
center_of_mass = center_of_mass / real(num_cells_local, wp)
Copy link
Contributor

Choose a reason for hiding this comment

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

P1: After the MPI reductions, the centroid is divided by num_cells_local (the pre-reduction local count) instead of the globally reduced num_cells, which can produce NaNs or incorrect centers of mass. Divide by the global count returned by s_mpi_allreduce_integer_sum (and guard against zero) to compute a valid centroid.

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

<comment>After the MPI reductions, the centroid is divided by `num_cells_local` (the pre-reduction local count) instead of the globally reduced `num_cells`, which can produce NaNs or incorrect centers of mass. Divide by the global count returned by `s_mpi_allreduce_integer_sum` (and guard against zero) to compute a valid centroid.</comment>

<file context>
@@ -1080,6 +1081,60 @@ contains
+            call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2))
+            call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3))
+            call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
+            center_of_mass =  center_of_mass / real(num_cells_local, wp)
+
+            ! assign the centroid offset as a vector pointing from the true COM to the "centroid" in the input file and replace the current centroid
</file context>
Suggested change
center_of_mass = center_of_mass / real(num_cells_local, wp)
if (num_cells > 0) then
center_of_mass = center_of_mass / real(num_cells, wp)
else
center_of_mass = 0._wp
end if

SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic()
Copy link
Contributor

Choose a reason for hiding this comment

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

P1: Same issue as with vel: assigning ParamType.REAL.analytic() (a plain dict) will cause AttributeError when building _properties. Use ParamType.REAL instead.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At toolchain/mfc/run/case_dicts.py, line 371:

<comment>Same issue as with `vel`: assigning `ParamType.REAL.analytic()` (a plain dict) will cause `AttributeError` when building `_properties`. Use `ParamType.REAL` instead.</comment>

<file context>
@@ -368,9 +366,9 @@ def analytic(self):
+        SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
         SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
-        SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
+        SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic()
 
     for cmp_id, cmp in enumerate(["x", "y", "z"]):
</file context>
Suggested change
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL.analytic()
SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL


for dir_id in range(1, 4):
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
Copy link
Contributor

Choose a reason for hiding this comment

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

P1: Assigning ParamType.REAL.analytic() (which returns a plain dict) directly into SIMULATION is inconsistent with other entries that use ParamType Enum members. Downstream code that builds _properties = {k: v.value for k, v in ALL.items()} will fail with an AttributeError because dicts don't have a .value attribute. Use ParamType.REAL instead to maintain consistency.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At toolchain/mfc/run/case_dicts.py, line 369:

<comment>Assigning `ParamType.REAL.analytic()` (which returns a plain dict) directly into `SIMULATION` is inconsistent with other entries that use `ParamType` Enum members. Downstream code that builds `_properties = {k: v.value for k, v in ALL.items()}` will fail with an `AttributeError` because dicts don't have a `.value` attribute. Use `ParamType.REAL` instead to maintain consistency.</comment>

<file context>
@@ -368,9 +366,9 @@ def analytic(self):
 
     for dir_id in range(1, 4):
-        SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL
+        SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
         SIMULATION[f"patch_ib({ib_id})%angles({dir_id})"] = ParamType.REAL
-        SIMULATION[f"patch_ib({ib_id})%angular_vel({dir_id})"] = ParamType.REAL
</file context>
Suggested change
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL.analytic()
SIMULATION[f"patch_ib({ib_id})%vel({dir_id})"] = ParamType.REAL

xy_local = matmul(inverse_rotation, xy_local)

! we will get NaNs in the levelset if we compute this outside the ellipse
if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
Copy link
Contributor

Choose a reason for hiding this comment

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

P1: levelset and normal fields are never set for cells outside the ellipse, leaving undefined geometry data for most of the domain. Add an else branch that computes the signed distance/normal outside the ellipse instead of skipping the assignment.

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

<comment>levelset and normal fields are never set for cells outside the ellipse, leaving undefined geometry data for most of the domain. Add an `else` branch that computes the signed distance/normal outside the ellipse instead of skipping the assignment.</comment>

<file context>
@@ -336,6 +339,62 @@ contains
+                xy_local = matmul(inverse_rotation, xy_local)
+
+                ! we will get NaNs in the levelset if we compute this outside the ellipse
+                if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
+
+                    normal_vector = xy_local
</file context>


# function that defines how we will replace variable names with
# values from the case file
def rhs_replace(match):
Copy link
Contributor

Choose a reason for hiding this comment

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

P1: Analytic IB codegen replaces x, y, z with grid-center arrays using undefined indices (j, k) inside the time-stepper loop, causing compilation failures or wrong values when users reference spatial variables.

Prompt for AI agents
Check if this issue is valid — if so, understand the root cause and fix it. At toolchain/mfc/case.py, line 227:

<comment>Analytic IB codegen replaces `x`, `y`, `z` with grid-center arrays using undefined indices (`j`, `k`) inside the time-stepper loop, causing compilation failures or wrong values when users reference spatial variables.</comment>

<file context>
@@ -181,6 +196,72 @@ def rhs_replace(match):
+
+            # function that defines how we will replace variable names with
+            # values from the case file
+            def rhs_replace(match):
+                return {
+                    'x': 'x_cc(i)', 'y': 'y_cc(j)', 'z': 'z_cc(k)',
</file context>

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

🤖 Fix all issues with AI agents
In `@src/common/m_ib_patches.fpp`:
- Line 781: The inline comment "! Transferring the rectangle's centroid and
length information" is a copy-paste error; update it to refer to the ellipse
instead (e.g., "! Transferring the ellipse's centroid and axes/length
information") so it accurately describes the code handling ellipse centroid and
length/axes values; locate and change this exact comment string in
m_ib_patches.fpp (search for "Transferring the rectangle's centroid") and
replace "rectangle" with "ellipse" (optionally adjust "length" to "axes" for
clarity).
♻️ Duplicate comments (2)
src/common/m_ib_patches.fpp (1)

788-791: Fix copy-paste comment: "rectangle" should be "ellipse".

📝 Proposed fix
-        ! Checking whether the rectangle covers a particular cell in the
+        ! Checking whether the ellipse covers a particular cell in the
         ! domain and verifying whether the current patch has the permission
         ! to write to that cell. If both queries check out, the primitive
         ! variables of the current patch are assigned to this cell.
src/common/m_compute_levelset.fpp (1)

376-391: Level set only computed inside ellipse - exterior cells remain uninitialized.

The condition at line 377 restricts computation to points inside the ellipse. Cells outside will retain their initial values (likely zero from the caller's initialization), which may cause incorrect ghost point interpolation for cells near but outside the ellipse boundary.

Compare to s_circle_levelset (lines 49-62) which computes dist - radius for ALL grid points, providing a valid signed distance everywhere.

🐛 Proposed fix: Compute signed distance for all points

For points outside the ellipse, you could use an algebraic distance approximation:

-                if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
+                ! Algebraic distance approximation valid everywhere
+                levelset%sf(i, j, 0, ib_patch_id) = sqrt(ellipse_coeffs(1) * ellipse_coeffs(2)) * &
+                    (sqrt((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2) - 1._wp)
 
+                ! Normal vector computation only valid inside/near surface
+                if ((xy_local(1)/ellipse_coeffs(1))**2 + (xy_local(2)/ellipse_coeffs(2))**2 <= 1._wp) then
                     normal_vector = xy_local
                     normal_vector(2) = normal_vector(2)*(ellipse_coeffs(1)/ellipse_coeffs(2))**2._wp
                     normal_vector = normal_vector/sqrt(dot_product(normal_vector, normal_vector))
                     levelset_norm%sf(i, j, 0, ib_patch_id, :) = matmul(rotation, normal_vector)
-
-                    quadratic_coeffs(1) = ...
-                    ...
-                    levelset%sf(i, j, 0, ib_patch_id) = ...
                 end if

Alternatively, compute the exact closest-point-on-ellipse distance for exterior points, though this requires iterative root finding.

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

1115-1121: Previous critical bug is fixed, but consider guarding against zero cell count.

The division now correctly uses the globally reduced num_cells instead of num_cells_local. However, if an IB patch has no cells (e.g., misconfigured geometry), num_cells could be zero, causing a division-by-zero crash.

🛡️ Optional: Add zero-guard for robustness
             call s_mpi_allreduce_sum(center_of_mass(1), center_of_mass(1))
             call s_mpi_allreduce_sum(center_of_mass(2), center_of_mass(2))
             call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3))
             call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
+            if (num_cells == 0) then
+                call s_mpi_abort("s_compute_centroid_offset: No cells found for IB patch")
+            end if
             center_of_mass = center_of_mass/real(num_cells, wp)
src/common/m_compute_levelset.fpp (1)

356-357: Remove unused variables k and idx.

These variables are declared but never used in the subroutine.

♻️ Proposed fix
-        integer :: i, j, k !< Loop index variables
-        integer :: idx !< Shortest path direction indicator
+        integer :: i, j !< Loop index variables
📜 Review details

Configuration used: defaults

Review profile: CHILL

Plan: Pro

📥 Commits

Reviewing files that changed from the base of the PR and between 7bd28cd and 541e727.

📒 Files selected for processing (3)
  • src/common/m_compute_levelset.fpp
  • src/common/m_ib_patches.fpp
  • src/simulation/m_ibm.fpp
🧰 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/common/m_ib_patches.fpp
  • src/simulation/m_ibm.fpp
  • src/common/m_compute_levelset.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/common/m_ib_patches.fpp
  • src/simulation/m_ibm.fpp
  • src/common/m_compute_levelset.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_ibm.fpp
🧠 Learnings (4)
📚 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_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_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_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/common/m_compute_levelset.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). (8)
  • GitHub Check: Github (ubuntu, no-mpi, single, no-debug, false)
  • GitHub Check: Github (ubuntu, mpi, debug, false)
  • GitHub Check: Github (ubuntu, mpi, no-debug, false)
  • GitHub Check: Github (ubuntu, mpi, no-debug, true)
  • GitHub Check: Github (ubuntu, mpi, debug, true)
  • GitHub Check: Code Cleanliness Check
  • GitHub Check: Coverage Test on CodeCov
  • GitHub Check: Build & Publish
🔇 Additional comments (8)
src/simulation/m_ibm.fpp (2)

1084-1136: LGTM! Well-structured centroid offset computation.

The subroutine correctly:

  • Targets only geometry types requiring offset computation (airfoils, STL models)
  • Uses MPI reduction for distributed computation
  • Transforms the offset back to local IB coordinates for consistent usage
  • Initializes offset to zero for inapplicable geometries

96-105: LGTM! Correct integration of centroid offset computation.

The call ordering is correct: rotation matrix is computed first (required for the offset rotation at line 1131), and the centroid offset is computed before the GPU data transfer.

src/common/m_compute_levelset.fpp (2)

95-97: LGTM! Correct coordinate frame for centroid offset subtraction.

The centroid_offset is stored in local IB coordinates (transformed at m_ibm.fpp:1131), and here it's correctly subtracted after xy_local has been rotated into the same local frame.


191-193: LGTM! Consistent with 2D airfoil implementation.

src/common/m_ib_patches.fpp (4)

280-282: LGTM! Centroid offset correctly applied in local coordinates.

The centroid_offset is pre-rotated into local IB coordinates during initialization (m_ibm.fpp:1131). Here it's correctly subtracted after xy_local is transformed to the same local frame. Both vectors share the same coordinate basis.


435-437: LGTM! Consistent centroid offset handling for 3D airfoil.


770-809: LGTM! Clean ellipse IB marker implementation.

The subroutine correctly:

  • Extracts ellipse parameters from patch data
  • Transforms coordinates to local IB frame
  • Uses standard ellipse membership test (x/a)² + (y/b)² ≤ 1

122-124: LGTM! Correct dispatch for ellipse geometry.

✏️ Tip: You can disable this entire section by setting review_details to false in your review settings.

@codeant-ai
Copy link
Contributor

codeant-ai bot commented Jan 16, 2026

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

@codeant-ai codeant-ai bot added size:L This PR changes 100-499 lines, ignoring generated files and removed size:L This PR changes 100-499 lines, ignoring generated files labels Jan 16, 2026
Comment on lines 96 to +97
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
xy_local = xy_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Coordinate-frame mixup: centroid_offset is subtracted after applying inverse_rotation, which treats the offset as if it were expressed in the rotated/local frame; if centroid_offset is defined in the global frame (as center is), this produces incorrect local coordinates and wrong nearest-point logic. Subtract the centroid_offset before rotation so the subtraction happens in the same (global) coordinate frame as the initial subtraction from center. [logic error]

Severity Level: Critical 🚨
- ❌ Airfoil 2D levelset geometry computed incorrectly.
- ⚠️ Nearest-point selection for airfoil grids unreliable.
- ⚠️ Normal vectors for IB forcing may be incorrect.
- ⚠️ Affects pitching airfoil test cases using centroid_offset.
Suggested change
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
xy_local = xy_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
xy_local = xy_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordinate
Steps of Reproduction ✅
1. Build and run the pitching airfoil test described in the PR description (the provided
pitching-airfoil configuration used to validate the change). The levelset computation for
the airfoil is invoked from module m_compute_levelset, subroutine s_airfoil_levelset in
src/common/m_compute_levelset.fpp.

2. Execution enters s_airfoil_levelset and iterates the 2D grid (do i = 0, m; do j = 0,
n). See code computing local coordinates at src/common/m_compute_levelset.fpp lines 95-97
where xy_local is set, rotated, then centroid_offset is subtracted.

3. If the airfoil patch has a non-zero patch_ib(ib_patch_id)%centroid_offset (set up by
the IB initialization for that patch), the subtraction at line 97 occurs after the
inverse_rotation at line 96, i.e., in the already-rotated/local frame. If centroid_offset
is provided in the global/frame-consistent coordinates (the usual convention for patch_ib
fields), this produces an incorrect xy_local.

4. Observe incorrect nearest-point selection and normals: levelset%sf and levelset_norm%sf
entries computed in s_airfoil_levelset (later in the same loop) will be shifted or have
incorrect normals. Reproduce by comparing levelset output or visualized geometry between
runs with centroid_offset==0 and centroid_offset!=0; the incorrect subtraction order
produces a displaced/rotated airfoil geometry in the test case.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/common/m_compute_levelset.fpp
**Line:** 96:97
**Comment:**
	*Logic Error: Coordinate-frame mixup: `centroid_offset` is subtracted after applying `inverse_rotation`, which treats the offset as if it were expressed in the rotated/local frame; if `centroid_offset` is defined in the global frame (as `center` is), this produces incorrect local coordinates and wrong nearest-point logic. Subtract the `centroid_offset` before rotation so the subtraction happens in the same (global) coordinate frame as the initial subtraction from `center`.

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.

moving_immersed_boundary_flag = .true.
end if
call s_update_ib_rotation_matrix(i)
call s_compute_centroid_offset(i)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Performance / correctness: s_compute_centroid_offset is called for every IB inside the initial setup loop before IB markers are populated; this causes an expensive full-domain scan per IB and also runs before ib_markers%sf is initialized/populated (leading to zero-cell counts). Remove the call from this early loop and compute centroid offsets later after IB patches/markers are applied and buffers populated. [performance]

Severity Level: Critical 🚨
- ❌ Startup initialization spends excessive CPU time.
- ⚠️ Centroid offsets computed from empty ib_markers.
- ❌ Ghost-point velocity uses wrong centroids.
- ⚠️ Pitching airfoil initial motion affected.
Suggested change
call s_compute_centroid_offset(i)
Steps of Reproduction ✅
1. Start the code path that calls the IBM setup routine: the public subroutine s_ibm_setup
in src/simulation/m_ibm.fpp is executed (see s_ibm_setup loop at lines 98-108).

2. Inside s_ibm_setup (lines 98-108) the code iterates do i = 1, num_ibs and
unconditionally calls s_compute_centroid_offset(i) at line 106 before any IB marker
population steps.

3. s_compute_centroid_offset (defined at lines 1155-1207) scans ib_markers%sf over the
full domain to count cells (it depends on ib_markers contents). At this point ib_markers
was only allocated earlier (s_initialize_ibm_module) but not populated by
s_apply_ib_patches / s_populate_ib_buffers yet, so entries are effectively empty/zero.

4. Observable effects when reproducing: you will see s_compute_centroid_offset perform a
full-domain nested loop (i=0..m, j=0..n, k=0..p) for every IB (costly) and compute zero
cell counts / zero centroid results because ib_markers%sf contains no IB assignments yet.
This reproduces wasted work and incorrect centroid results on startup.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_ibm.fpp
**Line:** 106:106
**Comment:**
	*Performance: Performance / correctness: s_compute_centroid_offset is called for every IB inside the initial setup loop before IB markers are populated; this causes an expensive full-domain scan per IB and also runs before `ib_markers%sf` is initialized/populated (leading to zero-cell counts). Remove the call from this early loop and compute centroid offsets later after IB patches/markers are applied and buffers populated.

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.

call s_mpi_allreduce_sum(center_of_mass(3), center_of_mass(3))
call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
center_of_mass = center_of_mass/real(num_cells, wp)
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggestion: Division-by-zero and incorrect averaging: the code reduces the summed coordinates and the integer cell count and unconditionally divides by num_cells. If no cells are found for the IB across all ranks this will divide by zero. Reduce the total cell-volume (or count) and check the global total before dividing; if zero, set a safe default and return. [possible bug]

Severity Level: Critical 🚨
- ❌ Simulation startup can crash (division by zero).
- ❌ Pitching airfoil initialisation fails.
- ⚠️ Any IB with geometry types 4,5,11,12 affected.
- ⚠️ MPI job may produce NaNs across ranks.
Suggested change
! also reduce the accumulated cell-volume across ranks and guard against zero total volume
call s_mpi_allreduce_sum(total_volume_local, total_volume)
if (total_volume <= 0._wp) then
! No volume found for this IB across all ranks: set zero offset and exit safely
patch_ib(ib_marker)%centroid_offset(:) = [0._wp, 0._wp, 0._wp]
return
end if
center_of_mass = center_of_mass/total_volume
Steps of Reproduction ✅
1. Trigger the startup path that runs s_ibm_setup (src/simulation/m_ibm.fpp lines 98-108)
so that s_compute_centroid_offset is invoked (calls at line 106).

2. Inside s_compute_centroid_offset (subroutine at lines 1155-1207) the routine tallies
num_cells_local on each rank and then calls s_mpi_allreduce_integer_sum to produce
num_cells (see lines 1188-1192).

3. If no grid cell was assigned to this IB on any rank (num_cells == 0 after allreduce)
the code executes center_of_mass = center_of_mass/real(num_cells, wp) at line 1192,
producing a division-by-zero (or invalid NaN) at runtime.

4. Reproduce by running a minimal case where the IB geometry guarded by the if-condition
(geometry == 4,5,11,12) is present but ib_markers are empty (this occurs during current
setup ordering), observe crash/NaN at startup originating from s_compute_centroid_offset
lines 1188-1192.
Prompt for AI Agent 🤖
This is a comment left during a code review.

**Path:** src/simulation/m_ibm.fpp
**Line:** 1193:1193
**Comment:**
	*Possible Bug: Division-by-zero and incorrect averaging: the code reduces the summed coordinates and the integer cell count and unconditionally divides by `num_cells`. If no cells are found for the IB across all ranks this will divide by zero. Reduce the total cell-volume (or count) and check the global total before dividing; if zero, set a safe default and return.

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.

@codeant-ai
Copy link
Contributor

codeant-ai bot commented Jan 16, 2026

CodeAnt AI Incremental review completed.

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.

Airfoil Centroid in Wrong Location Parametrized 1-Way MIBM Movement

1 participant