From cae2acb43835bbf8c8ea5ef504dde1b3e83e24e3 Mon Sep 17 00:00:00 2001 From: ruru25254649 Date: Thu, 22 Oct 2026 06:38:49 +0900 Subject: [PATCH 1/3] Fix the immersed-boundary moment of inertia Three errors in how a moving IB's moment of inertia is set: - s_ibm_setup computed it before the IB markers were filled. Geometries without a closed form (airfoil, STL, cuboid, cylinder, ...) sum over the body's marked cells, found none, and stored 0/0 = NaN. In 2D that is the only time the moment is computed; in 3D the time stepper reads it before recomputing it. A two-way-coupled airfoil, or a cuboid started with a nonzero angular velocity, went NaN. 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 (#1994). - The rectangle used m*(lx^2 + ly^2)/6. length_x and length_y are full side lengths (f_is_inside_cuboid tests |x| <= length_x/2), so the moment about the centroid is m*(lx^2 + ly^2)/12, the same convention the ellipse branch already uses. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WPhkY8p1EHhGY7jScizJmC --- src/simulation/m_ibm.fpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 614d0dc4b5..8f1491f259 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 real(wp), dimension(3) :: position, closest_point_along_axis, vector_to_axis, normal_axis integer :: i, j, k, count, ib_marker @@ -1575,19 +1581,13 @@ 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 if (p == 0) then @@ -1628,7 +1628,7 @@ contains end do ! write the final moment assuming the points are all uniform density - moment = moment*patch%mass/(count*cell_volume) + moment = moment*patch%mass/count end if end subroutine s_compute_moment_of_inertia From 38693260caff121d229f737ccbbd292299626557 Mon Sep 17 00:00:00 2001 From: ruru25254649 Date: Sun, 25 Oct 2026 16:36:59 +0900 Subject: [PATCH 2/3] Give a rank holding no cell of the body the dummy moment, not 0/0 Same guard and dummy value as #1992 (the non-rotating case already returns 1), so the two changes agree on this block. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WPhkY8p1EHhGY7jScizJmC --- src/simulation/m_ibm.fpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 8f1491f259..9daf187724 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1627,8 +1627,13 @@ contains end do end do - ! write the final moment assuming the points are all uniform density - moment = moment*patch%mass/count + ! 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 From f93805d9abc545276d15116572a4a15824d88355 Mon Sep 17 00:00:00 2001 From: ruru25254649 Date: Sun, 25 Oct 2026 17:51:14 +0900 Subject: [PATCH 3/3] Treat only a zero axis as non-rotating in the moment of inertia In 3D the time stepper passes the angular momentum as the axis. Comparing its norm with sgm_eps = 1e-16 is a dimensional threshold. A two-way-coupled cuboid with I = 6e-14 kg m^2 spinning at 1e-3 rad/s has L = 6e-17, so the first RK stage gets the dummy moment 1. The next stage multiplies the blended angular velocity by that dummy and divides by the real moment again: omega_z goes from 1e-3 to 8.3e9 rad/s in one step and the run stops on ICFL. Only the exact zero vector now returns the dummy value. The axis is scaled by its largest component before normalising, so a small one does not underflow. The same case keeps omega_z = 1.000e-3 for 20 steps. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01WPhkY8p1EHhGY7jScizJmC --- src/simulation/m_ibm.fpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 9daf187724..1dd32de4ff 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1573,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 + 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 @@ -1590,14 +1590,17 @@ contains count = 0 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