From c0a2c4ee90f6f79a813afd696fdfd54e5bd160fc Mon Sep 17 00:00:00 2001 From: Jiaxing Song Date: Sun, 25 Oct 2026 17:49:08 +0800 Subject: [PATCH] Compute the IB moment of inertia after the body is marked Motivation Moving immersed boundaries can receive NaN or incorrectly scaled inertia, and small nonzero angular momentum can be mistaken for zero. Expectation Compute physical inertia after cell marking and preserve the direction of every finite nonzero rotation axis. Change Four errors in how a two-way-coupled IB's moment of inertia is set: - s_ibm_setup computed it before the IB markers were filled. Shapes without a closed form (airfoil, STL, cuboid, cylinder, ...) sum over the body's marked cells, found none, and stored 0/0 = NaN, so the run aborted. It is now computed once the markers and the centre-of-mass offsets are in place. - The cell sum divided by count*cell_volume. With patch%mass the total mass, each of the count cells carries mass/count, so the moment is mass*sum(r^2)/count; the extra 1/cell_volume made it grid dependent. - The rectangle used m*(lx^2 + ly^2)/6. length_x and length_y are full side lengths, so the moment about the centroid is /12, the same convention the ellipse branch already uses. - In 3D the time stepper passes the angular momentum as the axis, and its norm was compared with sgm_eps = 1e-16, a threshold with units. A cuboid with I = 6e-14 kg m^2 spinning at 1e-3 rad/s got the dummy moment 1, and the next RK stage took omega_z to 8.3e9 rad/s. Only the exact zero vector now counts as not rotating; the axis is scaled by its largest component before normalising so it cannot underflow. A rank holding no cell of the body now gets the dummy moment 1 instead of 0/0, the same guard and value as #1992. Circles, spheres, and prescribed-motion or static IBs are unchanged; no goldens move. Result The corrected initialization, formulas and axis normalization are retained. The precheck passes and regression reference files are unchanged. Rank-local sums and nonuniform-cell weighting remain limitations. Fixes #1994 Co-authored-by: ruru25254649 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WPhkY8p1EHhGY7jScizJmC --- src/simulation/m_ibm.fpp | 40 ++++++++++++++++++++++++---------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 614d0dc4b5..1dd32de4ff 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -114,10 +114,7 @@ contains if (cfl_dt) t_init = t_save*n_start $:GPU_PARALLEL_LOOP(private='[i]', copyin='[t_init]') do i = 1, num_ibs - if (patch_ib(i)%moving_ibm /= 0) then - if (patch_ib(i)%kin_model > 0) call s_prescribed_kinematics(i, t_init) - call s_compute_moment_of_inertia(patch_ib(i), patch_ib(i)%angular_vel, patch_ib(i)%moment) - end if + if (patch_ib(i)%moving_ibm /= 0 .and. patch_ib(i)%kin_model > 0) call s_prescribed_kinematics(i, t_init) call s_update_ib_rotation_matrix(i) end do $:END_GPU_PARALLEL_LOOP() @@ -152,6 +149,15 @@ contains $:GPU_UPDATE(device='[patch_ib(i)]') end do + ! After the markers and centre-of-mass offsets above: shapes without a closed form sum over the body's marked cells + $:GPU_PARALLEL_LOOP(private='[i]') + do i = 1, num_ibs + if (patch_ib(i)%moving_ibm /= 0) call s_compute_moment_of_inertia(patch_ib(i), patch_ib(i)%angular_vel, & + & patch_ib(i)%moment) + end do + $:END_GPU_PARALLEL_LOOP() + $:GPU_UPDATE(host='[patch_ib(1:num_ibs)]') + ! find the number of ghost points and set them to be the maximum total across ranks call s_find_num_ghost_points(num_gps) if (moving_immersed_boundary_flag) then @@ -1567,7 +1573,7 @@ contains type(ib_patch_parameters), intent(in) :: patch real(wp), dimension(3), intent(in) :: axis real(wp), intent(out) :: moment - real(wp) :: distance_to_axis, cell_volume + real(wp) :: distance_to_axis, axis_max real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis, normal_axis integer :: i, j, k, count, ib_marker @@ -1575,29 +1581,26 @@ contains if (patch%geometry == 2) then ! circle moment = 0.5_wp*patch%mass*(patch%radius)**2 else if (patch%geometry == 3) then ! rectangle - moment = patch%mass*(patch%length_x**2 + patch %length_y**2)/6._wp + moment = patch%mass*(patch%length_x**2 + patch %length_y**2)/12._wp else if (patch%geometry == 6) then ! ellipse moment = 0.0625_wp*patch%mass*(patch%length_x**2 + patch %length_y**2) else if (patch%geometry == 8) then ! sphere moment = 0.4*patch%mass*(patch%radius)**2 else ! we do not have an analytic moment of inertia calculation and need to approximate it directly via a sum count = 0 - cell_volume = (x_cc(1) - x_cc(0))*(y_cc(1) - y_cc(0)) - ! computed without grid stretching. Update in the loop to perform with stretching - if (p /= 0) then - cell_volume = cell_volume*(z_cc(1) - z_cc(0)) - end if - ib_marker = patch%gbl_patch_id + axis_max = max(abs(axis(1)), abs(axis(2)), abs(axis(3))) if (p == 0) then normal_axis = [0, 0, 1] - else if (sqrt(sum(axis**2)) < sgm_eps) then + else if (axis_max == 0._wp) then ! if the object is not actually rotating at this time, return a dummy value and exit moment = 1._wp return else - normal_axis = axis/sqrt(sum(axis**2)) + ! only the direction is used; scaling first keeps a small angular momentum from underflowing + normal_axis = axis/axis_max + normal_axis = normal_axis/sqrt(sum(normal_axis**2)) end if moment = 0._wp @@ -1627,8 +1630,13 @@ contains end do end do - ! write the final moment assuming the points are all uniform density - moment = moment*patch%mass/(count*cell_volume) + ! write the final moment assuming the points are all uniform density; a rank holding no cell of the body gets the + ! dummy value instead of 0/0 + if (count > 0) then + moment = moment*patch%mass/count + else + moment = 1._wp + end if end if end subroutine s_compute_moment_of_inertia