Skip to content

Analytic mibm velocities and airfoil centroid #1111

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
danieljvickers wants to merge 18 commits into
base: master
Choose a base branch
from
+238 −62
Open

Analytic mibm velocities and airfoil centroid #1111

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
89a0ac5
Added ellipse marker patch
danieljvickers Dec 12, 2025
ef446f5
Added ellipse IB patch
danieljvickers Dec 12, 2025
2f9b1be
Added an ellipse immersed boundary patch
danieljvickers Dec 12, 2025
d54bc5b
initial implementation of analytic moving IB
danieljvickers Jan 14, 2026
697b5a0
Works in simple example
danieljvickers Jan 14, 2026
6ed8350
Added beginning of centroid offset for calculations
danieljvickers Jan 15, 2026
7bd28cd
Works with pitching airfoil
danieljvickers Jan 15, 2026
f40422a
Spelling and formatting
danieljvickers Jan 15, 2026
541e727
Hey, the AI finally caught a bug
danieljvickers Jan 15, 2026
cb7b845
Added empty analytic function that can be found in documentation
danieljvickers Jan 15, 2026
54c9421
Merge branch 'master' into analytic-mibm-velocities-and-airfoil-centroid
danieljvickers Jan 16, 2026
37f2730
Fixed multi-rank issue on frontier
Jan 16, 2026
8b00a82
fpp files must end in a empty line
danieljvickers Jan 17, 2026
6817228
Merge branch 'analytic-mibm-velocities-and-airfoil-centroid' of githu…
danieljvickers Jan 17, 2026
f6c961e
formatting
danieljvickers Jan 17, 2026
a815445
Refactored the IB code into it's own subroutine
danieljvickers Jan 18, 2026
e94637d
Updated the force calculations to work with multiple fluids
danieljvickers Jan 18, 2026
675d974
logic corrections to dynamic viscosity calculation
danieljvickers Jan 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions src/common/include/case.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,3 +6,8 @@
#:def analytical()

#:enddef

! For moving immersed boundaries in simulation
#:def mib_analytical()

#:enddef
2 changes: 2 additions & 0 deletions src/common/m_compute_levelset.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,7 @@ contains
do j = 0, n
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 coordinate
xy_local = xy_local - patch_ib(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset

if (xy_local(2) >= 0._wp) then
! finds the location on the airfoil grid with the minimum distance (closest)
Expand Down Expand Up @@ -189,6 +190,7 @@ contains

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(ib_patch_id)%centroid_offset ! airfoils are a patch that require a centroid offset

if (xyz_local(2) >= center(2)) then
do k = 1, Np
Expand Down
1 change: 1 addition & 0 deletions src/common/m_derived_types.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -296,6 +296,7 @@ module m_derived_types
!! is specified through its x-, y- and z-coordinates, respectively.
real(wp) :: step_x_centroid, step_y_centroid, step_z_centroid !<
!! Centroid locations of intermediate steps in the time_stepper module
real(wp), dimension(1:3) :: centroid_offset ! offset of center of mass from computed cell center for odd-shaped IBs

real(wp), dimension(1:3) :: angles
real(wp), dimension(1:3) :: step_angles
Expand Down
2 changes: 2 additions & 0 deletions src/common/m_ib_patches.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -279,6 +279,7 @@ contains
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

if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then
xa = xy_local(1)/ca_in
Expand Down Expand Up @@ -433,6 +434,7 @@ contains
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

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


if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then

Expand Down
1 change: 1 addition & 0 deletions src/pre_process/m_global_parameters.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -572,6 +572,7 @@ contains
patch_ib(i)%angular_vel(:) = 0._wp
patch_ib(i)%mass = dflt_real
patch_ib(i)%moment = dflt_real
patch_ib(i)%centroid_offset(:) = 0._wp

! sets values of a rotation matrix which can be used when calculating rotations
patch_ib(i)%rotation_matrix = 0._wp
Expand Down
99 changes: 81 additions & 18 deletions src/simulation/m_ibm.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,13 +95,15 @@ contains
integer :: i, j, k
integer :: max_num_gps, max_num_inner_gps

! 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
call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel)
moving_immersed_boundary_flag = .true.
end if
call s_update_ib_rotation_matrix(i)
call s_compute_centroid_offset(i)
Copy link
Contributor

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

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.

end do
$:GPU_ENTER_DATA(copyin='[patch_ib]')

Expand Down Expand Up @@ -1012,24 +1014,24 @@ contains
! compute the surface integrals of the IB via a volume integraion method described in
! "A coupled IBM/Euler-Lagrange framework for simulating shock-induced particle size segregation"
! by Archana Sridhar and Jesse Capecelatro
subroutine s_compute_ib_forces(q_prim_vf, dynamic_viscosity)
subroutine s_compute_ib_forces(q_prim_vf, fluid_pp)

! real(wp), dimension(idwbuff(1)%beg:idwbuff(1)%end, &
! idwbuff(2)%beg:idwbuff(2)%end, &
! idwbuff(3)%beg:idwbuff(3)%end), intent(in) :: pressure
type(scalar_field), dimension(1:sys_size), intent(in) :: q_prim_vf
real(wp), intent(in) :: dynamic_viscosity
type(physical_parameters), dimension(1:num_fluids), intent(in) :: fluid_pp

integer :: gp_id, i, j, k, l, q, ib_idx
integer :: gp_id, i, j, k, l, q, ib_idx, fluid_idx
real(wp), dimension(num_ibs, 3) :: forces, torques
real(wp), dimension(1:3, 1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2 ! viscous stress tensor with temp vectors to hold divergence calculations
real(wp), dimension(1:3) :: local_force_contribution, radial_vector, local_torque_contribution, vel
real(wp) :: cell_volume, dx, dy, dz
real(wp) :: cell_volume, dx, dy, dz, dynamic_viscosity

forces = 0._wp
torques = 0._wp

$:GPU_PARALLEL_LOOP(private='[ib_idx,radial_vector,local_force_contribution,cell_volume,local_torque_contribution, viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2, dx, dy, dz]', copy='[forces,torques]', copyin='[ib_markers,patch_ib,dynamic_viscosity]', collapse=3)
$:GPU_PARALLEL_LOOP(private='[ib_idx,fluid_idx, radial_vector,local_force_contribution,cell_volume,local_torque_contribution, dynamic_viscosity, viscous_stress_div, viscous_stress_div_1, viscous_stress_div_2, viscous_cross_1, viscous_cross_2, dx, dy, dz]', copy='[forces,torques]', copyin='[ib_markers,patch_ib]', collapse=3)
do i = 0, m
do j = 0, n
do k = 0, p
Expand All @@ -1044,26 +1046,33 @@ contains
dx = x_cc(i + 1) - x_cc(i)
dy = y_cc(j + 1) - y_cc(j)

! Get the pressure contribution to force via a finite difference to compute the 2D components of the gradient of the pressure and cell volume
local_force_contribution(1) = -1._wp*(q_prim_vf(E_idx)%sf(i + 1, j, k) - q_prim_vf(E_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
local_force_contribution(2) = -1._wp*(q_prim_vf(E_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx)%sf(i, j - 1, k))/(2._wp*dy)
cell_volume = abs(dx*dy)
! add the 3D component of the pressure gradient, if we are working in 3 dimensions
if (num_dims == 3) then
dz = z_cc(k + 1) - z_cc(k)
local_force_contribution(3) = -1._wp*(q_prim_vf(E_idx)%sf(i, j, k + 1) - q_prim_vf(E_idx)%sf(i, j, k - 1))/(2._wp*dz)
cell_volume = abs(cell_volume*dz)
else
local_force_contribution(3) = 0._wp
end if
local_force_contribution(:) = 0._wp
do fluid_idx = 0, num_fluids - 1
! Get the pressure contribution to force via a finite difference to compute the 2D components of the gradient of the pressure and cell volume
local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i + 1, j, k) - q_prim_vf(E_idx + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
local_force_contribution(2) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 18, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

P1: Pressure-gradient accumulation overwrites Y/Z components with the X component, producing incorrect force/torque directions. Each component should accumulate from its own value.

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

<comment>Pressure-gradient accumulation overwrites Y/Z components with the X component, producing incorrect force/torque directions. Each component should accumulate from its own value.</comment>

<file context>
@@ -1046,26 +1046,33 @@ contains
+                        do fluid_idx = 0, num_fluids - 1
+                            ! Get the pressure contribution to force via a finite difference to compute the 2D components of the gradient of the pressure and cell volume
+                            local_force_contribution(1) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i + 1, j, k) - q_prim_vf(E_idx + fluid_idx)%sf(i - 1, j, k))/(2._wp*dx) ! force is the negative pressure gradient
+                            local_force_contribution(2) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
+                            cell_volume = abs(dx*dy)
+                            ! add the 3D component of the pressure gradient, if we are working in 3 dimensions
</file context>
Suggested change
local_force_contribution(2) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
Fix with Cubic

cell_volume = abs(dx*dy)
! add the 3D component of the pressure gradient, if we are working in 3 dimensions
if (num_dims == 3) then
dz = z_cc(k + 1) - z_cc(k)
local_force_contribution(3) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j, k + 1) - q_prim_vf(E_idx + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
cell_volume = abs(cell_volume*dz)
Comment on lines +1049 to +1059
Copy link
Contributor

@coderabbitai coderabbitai bot Jan 18, 2026

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Fix pressure-gradient accumulation for y/z components.
Line 1053 and Line 1058 use local_force_contribution(1) on the RHS, which overwrites components 2/3 and mixes axes. This yields incorrect force vectors.

🐛 Proposed fix 
-                            local_force_contribution(2) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
+                            local_force_contribution(2) = local_force_contribution(2) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j + 1, k) - q_prim_vf(E_idx + fluid_idx)%sf(i, j - 1, k))/(2._wp*dy)
@@
-                                local_force_contribution(3) = local_force_contribution(1) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j, k + 1) - q_prim_vf(E_idx + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
+                                local_force_contribution(3) = local_force_contribution(3) - (q_prim_vf(E_idx + fluid_idx)%sf(i, j, k + 1) - q_prim_vf(E_idx + fluid_idx)%sf(i, j, k - 1))/(2._wp*dz)
🤖 Prompt for AI Agents 
In `@src/simulation/m_ibm.fpp` around lines 1049 - 1059, The y and z components of
the pressure-force accumulator are accidentally using
local_force_contribution(1) on the right-hand side, which mixes axes; in the
block that computes local_force_contribution(2) and, inside the num_dims==3
branch, local_force_contribution(3), change the RHS to reference the same
component being updated (use local_force_contribution(2) for the y update and
local_force_contribution(3) for the z update) so each axis subtracts its own
finite-difference pressure gradient (retain use of q_prim_vf(E_idx +
fluid_idx)%sf(...) and the existing dx, dy, dz and cell_volume logic).

end if
end do

! Update the force values atomically to prevent race conditions
call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)

! get the viscous stress and add its contribution if that is considered
! TODO :: This is really bad code
! if (.false.) then
if (viscous) then
! compute the volume-weighted local dynamic viscosity
dynamic_viscosity = 0._wp
do fluid_idx = 1, num_fluids
! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
if (fluid_pp(fluid_idx)%Re(1) /= 0._wp) dynamic_viscosity = dynamic_viscosity + q_prim_vf(fluid_idx + advxb - 1)%sf(j, k, l)*(1._wp/fluid_pp(fluid_idx)%Re(1))
Copy link
Contributor

@cubic-dev-ai cubic-dev-ai bot Jan 18, 2026 •
edited
Loading

Choose a reason for hiding this comment

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

P2: The dynamic viscosity computation indexes q_prim_vf with (j, k, l) even though l is not initialized in this loop. This can read the wrong cell (or an undefined index), producing incorrect forces/torques. Use the current cell indices (i, j, k).

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

<comment>The dynamic viscosity computation indexes q_prim_vf with (j, k, l) even though l is not initialized in this loop. This can read the wrong cell (or an undefined index), producing incorrect forces/torques. Use the current cell indices (i, j, k).</comment>

<file context>
@@ -1070,7 +1070,7 @@ contains
                             do fluid_idx = 1, num_fluids
                                 ! local dynamic viscosity is the dynamic viscosity of the fluid times alpha of the fluid
-                                if (fluid_pp%Re(fluid_idx) /= 0._wp) dynamic_viscosity = dynamic_viscosity + q_cons_vf(fluid_idx + advxb - 1)%sf(j, k, l)*(1._wp/fluid_pp(fluid_idx)%Re(1))
+                                if (fluid_pp(fluid_idx)%Re(1) /= 0._wp) dynamic_viscosity = dynamic_viscosity + q_prim_vf(fluid_idx + advxb - 1)%sf(j, k, l)*(1._wp/fluid_pp(fluid_idx)%Re(1))
                             end do
 
</file context>
Suggested change
if (fluid_pp(fluid_idx)%Re(1) /= 0._wp) dynamic_viscosity = dynamic_viscosity + q_prim_vf(fluid_idx + advxb - 1)%sf(j, k, l)*(1._wp/fluid_pp(fluid_idx)%Re(1))
if (fluid_pp(fluid_idx)%Re(1) /= 0._wp) dynamic_viscosity = dynamic_viscosity + q_prim_vf(fluid_idx + advxb - 1)%sf(i, j, k)*(1._wp/fluid_pp(fluid_idx)%Re(1))
Fix with Cubic

end do

! get the linear force component first
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i - 1, j, k)
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i + 1, j, k)
Expand Down Expand Up @@ -1150,6 +1159,60 @@ contains

end subroutine s_finalize_ibm_module

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

! 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_local = [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_local = center_of_mass_local + [x_cc(i), y_cc(j), 0._wp]
if (num_dims == 3) center_of_mass_local(3) = center_of_mass_local(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_local(1), center_of_mass(1))
call s_mpi_allreduce_sum(center_of_mass_local(2), center_of_mass(2))
call s_mpi_allreduce_sum(center_of_mass_local(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

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

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.

Comment on lines +1193 to +1200
Copy link
Contributor

@coderabbitai coderabbitai bot Jan 16, 2026

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Add guard for division by zero when num_cells == 0.

The previous bug dividing by num_cells_local instead of num_cells was correctly fixed. However, if no cells are found for this IB across all MPI ranks (e.g., IB outside computational domain), num_cells will be 0, causing a floating-point exception or NaN propagation.

🛡️ Proposed guard 
             call s_mpi_allreduce_integer_sum(num_cells_local, num_cells)
+            if (num_cells == 0) then
+                ! No cells found for this IB - set zero offset and return
+                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)
🤖 Prompt for AI Agents 
In `@src/simulation/m_ibm.fpp` around lines 1186 - 1193, After the MPI reductions
(s_mpi_allreduce_sum on center_of_mass_local and s_mpi_allreduce_integer_sum on
num_cells) add a guard to avoid dividing by zero: check if num_cells == 0 and in
that case set center_of_mass = 0.0_wp (or another defined safe default) and skip
the division; otherwise perform center_of_mass = center_of_mass/real(num_cells,
wp). This change touches the variables center_of_mass, center_of_mass_local,
num_cells_local, num_cells and the reduction sequence (s_mpi_allreduce_sum /
s_mpi_allreduce_integer_sum) so place the conditional immediately after the
reductions and before the division.

! 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

end subroutine s_compute_centroid_offset

subroutine s_compute_moment_of_inertia(ib_marker, axis)

real(wp), dimension(3), intent(in) :: axis !< the axis about which we compute the moment. Only required in 3D.
Expand Down
90 changes: 53 additions & 37 deletions src/simulation/m_time_steppers.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
!! @brief Contains module m_time_steppers

#:include 'macros.fpp'
#:include 'case.fpp'

!> @brief The following module features a variety of time-stepping schemes.
!! Currently, it includes the following Runge-Kutta (RK) algorithms:
Expand Down Expand Up @@ -611,44 +612,10 @@ contains
if (ib) then
! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms
if (moving_immersed_boundary_flag) then
do i = 1, num_ibs
if (s == 1) then
patch_ib(i)%step_vel = patch_ib(i)%vel
patch_ib(i)%step_angular_vel = patch_ib(i)%angular_vel
patch_ib(i)%step_angles = patch_ib(i)%angles
patch_ib(i)%step_x_centroid = patch_ib(i)%x_centroid
patch_ib(i)%step_y_centroid = patch_ib(i)%y_centroid
patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid
end if

if (patch_ib(i)%moving_ibm > 0) then
patch_ib(i)%vel = (rk_coef(s, 1)*patch_ib(i)%step_vel + rk_coef(s, 2)*patch_ib(i)%vel)/rk_coef(s, 4)
patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, 2)*patch_ib(i)%angular_vel)/rk_coef(s, 4)

if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque
! compute the force and torque on the IB from the fluid
call s_compute_ib_forces(q_prim_vf, 1._wp/fluid_pp(1)%Re(1))

! update the velocity from the force value
patch_ib(i)%vel = patch_ib(i)%vel + rk_coef(s, 3)*dt*(patch_ib(i)%force/patch_ib(i)%mass)/rk_coef(s, 4)

! update the angular velocity with the torque value
patch_ib(i)%angular_vel = (patch_ib(i)%angular_vel*patch_ib(i)%moment) + (rk_coef(s, 3)*dt*patch_ib(i)%torque/rk_coef(s, 4)) ! add the torque to the angular momentum
call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) ! update the moment of inertia to be based on the direction of the angular momentum
patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i)%moment ! convert back to angular velocity with the new moment of inertia
end if

! Update the angle of the IB
patch_ib(i)%angles = (rk_coef(s, 1)*patch_ib(i)%step_angles + rk_coef(s, 2)*patch_ib(i)%angles + rk_coef(s, 3)*patch_ib(i)%angular_vel*dt)/rk_coef(s, 4)

! Update the position of the IB
patch_ib(i)%x_centroid = (rk_coef(s, 1)*patch_ib(i)%step_x_centroid + rk_coef(s, 2)*patch_ib(i)%x_centroid + rk_coef(s, 3)*patch_ib(i)%vel(1)*dt)/rk_coef(s, 4)
patch_ib(i)%y_centroid = (rk_coef(s, 1)*patch_ib(i)%step_y_centroid + rk_coef(s, 2)*patch_ib(i)%y_centroid + rk_coef(s, 3)*patch_ib(i)%vel(2)*dt)/rk_coef(s, 4)
patch_ib(i)%z_centroid = (rk_coef(s, 1)*patch_ib(i)%step_z_centroid + rk_coef(s, 2)*patch_ib(i)%z_centroid + rk_coef(s, 3)*patch_ib(i)%vel(3)*dt)/rk_coef(s, 4)
end if
end do
call s_update_mib(num_ibs, levelset, levelset_norm)
call s_propagate_immersed_boundaries(s)
end if

! update the ghost fluid properties point values based on IB state
if (qbmm .and. .not. polytropic) then
call s_ibm_correct_state(q_cons_ts(1)%vf, q_prim_vf, pb_ts(1)%sf, mv_ts(1)%sf)
else
Expand Down Expand Up @@ -802,6 +769,55 @@ contains

end subroutine s_apply_bodyforces

subroutine s_propagate_immersed_boundaries(s)

integer, intent(in) :: s
integer :: i

do i = 1, num_ibs
if (s == 1) then
patch_ib(i)%step_vel = patch_ib(i)%vel
patch_ib(i)%step_angular_vel = patch_ib(i)%angular_vel
patch_ib(i)%step_angles = patch_ib(i)%angles
patch_ib(i)%step_x_centroid = patch_ib(i)%x_centroid
patch_ib(i)%step_y_centroid = patch_ib(i)%y_centroid
patch_ib(i)%step_z_centroid = patch_ib(i)%z_centroid
end if

if (patch_ib(i)%moving_ibm > 0) then
patch_ib(i)%vel = (rk_coef(s, 1)*patch_ib(i)%step_vel + rk_coef(s, 2)*patch_ib(i)%vel)/rk_coef(s, 4)
patch_ib(i)%angular_vel = (rk_coef(s, 1)*patch_ib(i)%step_angular_vel + rk_coef(s, 2)*patch_ib(i)%angular_vel)/rk_coef(s, 4)

if (patch_ib(i)%moving_ibm == 1) then
! plug in analytic velocities for 1-way coupling, if it exists
@:mib_analytical()
else if (patch_ib(i)%moving_ibm == 2) then ! if we are using two-way coupling, apply force and torque
! compute the force and torque on the IB from the fluid
call s_compute_ib_forces(q_prim_vf, fluid_pp)

! update the velocity from the force value
patch_ib(i)%vel = patch_ib(i)%vel + rk_coef(s, 3)*dt*(patch_ib(i)%force/patch_ib(i)%mass)/rk_coef(s, 4)

! update the angular velocity with the torque value
patch_ib(i)%angular_vel = (patch_ib(i)%angular_vel*patch_ib(i)%moment) + (rk_coef(s, 3)*dt*patch_ib(i)%torque/rk_coef(s, 4)) ! add the torque to the angular momentum
call s_compute_moment_of_inertia(i, patch_ib(i)%angular_vel) ! update the moment of inertia to be based on the direction of the angular momentum
patch_ib(i)%angular_vel = patch_ib(i)%angular_vel/patch_ib(i)%moment ! convert back to angular velocity with the new moment of inertia
end if

! Update the angle of the IB
patch_ib(i)%angles = (rk_coef(s, 1)*patch_ib(i)%step_angles + rk_coef(s, 2)*patch_ib(i)%angles + rk_coef(s, 3)*patch_ib(i)%angular_vel*dt)/rk_coef(s, 4)

! Update the position of the IB
patch_ib(i)%x_centroid = (rk_coef(s, 1)*patch_ib(i)%step_x_centroid + rk_coef(s, 2)*patch_ib(i)%x_centroid + rk_coef(s, 3)*patch_ib(i)%vel(1)*dt)/rk_coef(s, 4)
patch_ib(i)%y_centroid = (rk_coef(s, 1)*patch_ib(i)%step_y_centroid + rk_coef(s, 2)*patch_ib(i)%y_centroid + rk_coef(s, 3)*patch_ib(i)%vel(2)*dt)/rk_coef(s, 4)
patch_ib(i)%z_centroid = (rk_coef(s, 1)*patch_ib(i)%step_z_centroid + rk_coef(s, 2)*patch_ib(i)%z_centroid + rk_coef(s, 3)*patch_ib(i)%vel(3)*dt)/rk_coef(s, 4)
end if
end do

call s_update_mib(num_ibs, levelset, levelset_norm)

end subroutine s_propagate_immersed_boundaries

!> This subroutine saves the temporary q_prim_vf vector
!! into the q_prim_ts vector that is then used in p_main
!! @param t_step current time-step
Expand Down
Loading