From 33ee6b601890b081a8dded58095d5d590e5f5aef Mon Sep 17 00:00:00 2001 From: Mike VandenBoom Date: Sat, 31 Jan 2026 21:58:16 -0500 Subject: [PATCH 1/3] Force calculation improvements (cherry picked from commit 2a034dc397321beb956bc30c102e42edbd44c2dc) (cherry picked from commit 1d8b1092c0c5f47636320370920612cfc1fadfba) --- src/simulation/m_ibm.fpp | 42 +++++++++++----------------------------- 1 file changed, 11 insertions(+), 31 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index bfe09ec5ac..8850034889 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -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 @@ -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 + 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 + end if + call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution) + + ! Update the force values atomically to prevent race conditions do l = 1, 3 $:GPU_ATOMIC(atomic='update') forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume) From c0d288c13fe9a42c97c0eb1b3a3c48ca7ebdb74a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 24 Feb 2026 15:25:20 -0500 Subject: [PATCH 2/3] fix format --- src/simulation/m_ibm.fpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 8850034889..a50073d3fa 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1127,7 +1127,7 @@ contains 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 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 - + end if call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution) From 7c5baa935cf350dfe40194f49413cc23e8b1618f Mon Sep 17 00:00:00 2001 From: Mike VandenBoom Date: Thu, 12 Mar 2026 17:43:31 -0400 Subject: [PATCH 3/3] Clean up edits recommended by Claude --- src/simulation/m_ibm.fpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index f8b1869ae6..87b7a67eff 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -997,7 +997,7 @@ contains 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, 1:3) :: viscous_stress_div, viscous_stress_div_1, viscous_stress_div_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, dynamic_viscosity #:if not MFC_CASE_OPTIMIZATION and USING_AMD @@ -1021,7 +1021,7 @@ contains end do end if - $: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,dynamic_viscosities]', 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, dx, dy, dz]', copy='[forces,torques]', copyin='[ib_markers,patch_ib,dynamic_viscosities]', collapse=3) do i = 0, m do j = 0, n do k = 0, p @@ -1051,7 +1051,6 @@ contains end do ! get the viscous stress and add its contribution if that is considered - ! TODO :: This is really bad code if (viscous) then ! compute the volume-weighted local dynamic viscosity dynamic_viscosity = 0._wp @@ -1074,7 +1073,7 @@ contains 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(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 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 @@ -1082,7 +1081,7 @@ contains call s_cross_product(radial_vector, local_force_contribution, local_torque_contribution) - ! Update the force values atomically to prevent race conditions + ! Update the force and torque values atomically to prevent race conditions do l = 1, 3 $:GPU_ATOMIC(atomic='update') forces(ib_idx, l) = forces(ib_idx, l) + (local_force_contribution(l)*cell_volume)