Skip to content

IBM force calculation simplification (Dr Bala recommendations) #1234

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
mrvandenboom wants to merge 4 commits into
base: master
Choose a base branch
from
Open

IBM force calculation simplification (Dr Bala recommendations) #1234

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
33ee6b6
Force calculation improvements
mrvandenboom Feb 1, 2026
615f99d
Merge branch 'master' into bala_force_changes
sbryngelson Feb 22, 2026
448d067
Merge branch 'master' into bala_force_changes
sbryngelson Feb 22, 2026
25ce36d
Merge branch 'master' into bala_force_changes
sbryngelson Feb 23, 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
42 changes: 11 additions & 31 deletions src/simulation/m_ibm.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1100,9 +1100,6 @@ contains
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 (viscous) then
Expand All @@ -1113,46 +1110,29 @@ contains
dynamic_viscosity = dynamic_viscosity + (q_prim_vf(fluid_idx + advxb - 1)%sf(i, j, k)*dynamic_viscosities(fluid_idx))
end do

! get the linear force component first
! get the linear force components 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)
viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dx) ! get the x derivative of the viscous stress tensor
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, 1:3) ! add te x components of the derivative to the force
do l = 1, 3
! take the cross products for the torque component
call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
end do

viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dx) ! get the x derivative of the cross product
local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(1, 1:3) ! apply the cross product derivative to the torque
viscous_stress_div(1, 1:3) = (viscous_stress_div_2(1, 1:3) - viscous_stress_div_1(1, 1:3))/(2._wp*dx) ! get x derivative of the first-row of viscous stress tensor
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(1, 1:3) ! add the x components of the divergence to the force

call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j - 1, k)
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j + 1, k)
viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dy)
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, 1:3)
do l = 1, 3
call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
end do

viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dy)
local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(2, 1:3)
viscous_stress_div(2, 1:3) = (viscous_stress_div_2(2, 1:3) - viscous_stress_div_1(2, 1:3))/(2._wp*dy) ! get y derivative of the second-row of viscous stress tensor
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(2, 1:3) ! add the y components of the divergence to the force

if (num_dims == 3) then
call s_compute_viscous_stress_tensor(viscous_stress_div_1, q_prim_vf, dynamic_viscosity, i, j, k - 1)
call s_compute_viscous_stress_tensor(viscous_stress_div_2, q_prim_vf, dynamic_viscosity, i, j, k + 1)
viscous_stress_div = (viscous_stress_div_2 - viscous_stress_div_1)/(2._wp*dz)
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, 1:3)
do l = 1, 3
call s_cross_product(radial_vector, viscous_stress_div_1(l, 1:3), viscous_cross_1(l, 1:3))
call s_cross_product(radial_vector, viscous_stress_div_2(l, 1:3), viscous_cross_2(l, 1:3))
end do
viscous_stress_div = (viscous_cross_2 - viscous_cross_1)/(2._wp*dz)
local_torque_contribution(1:3) = local_torque_contribution(1:3) + viscous_stress_div(3, 1:3)
viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the second-row of viscous stress tensor
Copy link

Copilot AI Feb 22, 2026

Choose a reason for hiding this comment

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

Inline comment is incorrect: this statement computes the z-derivative of the third row (index 3) of the viscous stress tensor, not the second row. Please update the comment to avoid confusion for future maintenance.

Suggested change
viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the second-row of viscous stress tensor
viscous_stress_div(3, 1:3) = (viscous_stress_div_2(3, 1:3) - viscous_stress_div_1(3, 1:3))/(2._wp*dz) ! get z derivative of the third-row of viscous stress tensor

Copilot uses AI. Check for mistakes.
local_force_contribution(1:3) = local_force_contribution(1:3) + viscous_stress_div(3, 1:3) ! add the z components of the divergence to the force
end if

Copy link
Contributor

@coderabbitai coderabbitai bot Feb 22, 2026

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Trailing whitespace — this is causing the CI formatting failure.

Line 1130 appears to contain trailing whitespace. The pipeline log confirms: formatting diff found in m_ibm.fpp. Remove the trailing spaces to pass the ./mfc.sh format check.

🔧 Fix trailing whitespace 
-                            
+
📝 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
🤖 Prompt for AI Agents 
Verify each finding against the current code and only fix it if needed.

In `@src/simulation/m_ibm.fpp` at line 1130, Remove the trailing whitespace on the
blank line flagged by CI (a line that contains only spaces) in the m_ibm.fpp
source so the formatter passes; locate the empty/blank line in the file, delete
the trailing spaces (convert it to a truly empty line), run ./mfc.sh format to
verify no formatting diffs remain, and commit the change.

end if

call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution)

! Update the force values atomically to prevent race conditions
Copy link

Copilot AI Feb 22, 2026

Choose a reason for hiding this comment

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

This comment says only the force is updated atomically, but the loop performs atomic updates for both forces and torques. Consider updating the comment to reflect both arrays to prevent misleading documentation.

Suggested change
! Update the force values atomically to prevent race conditions
! Update the force and torque values atomically to prevent race conditions

Copilot uses AI. Check for mistakes.
do l = 1, 3
$:GPU_ATOMIC(atomic='update')
forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)
Expand Down
Loading