From d21eccba86ad23467dd068a15f4dc8c86bae275a Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 29 Sep 2026 23:12:36 -0400 Subject: [PATCH 1/2] Let prescribed kinematics drive airfoil and STL immersed boundaries These geometries move their working centroid to the centre of mass and keep the difference in centroid_offset. Account for it everywhere it matters: - s_compute_centroid_offset is collective; call it for every global patch on every rank (#1895) - s_model_levelset subtracts centroid_offset, matching the marker test (#1896) - prescribed kinematics rotate kin_offset - centroid_offset about the hinge - checkpoints record the offsets and restarts restore them (#1902) The validator no longer forbids kin_model on geometries 4, 5, 11 and 12. Co-Authored-By: Claude --- docs/documentation/case.md | 2 + src/simulation/m_compute_levelset.fpp | 2 + src/simulation/m_data_output.fpp | 34 ++++++- src/simulation/m_ibm.fpp | 128 +++++++++++++++++--------- toolchain/mfc/case_validator.py | 7 -- 5 files changed, 124 insertions(+), 49 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 12bec5a38c..390794bf00 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -426,6 +426,8 @@ Additional details on this specification can be found in [NACA airfoil](https:// - `kin_model = 1` prescribes hinged flapping kinematics at run time (no analytic expressions, so the binary is shared across parameter values): roll \f$\phi\f$ about the lab \f$x\f$ axis through `kin_hinge` and pitch \f$\theta\f$ about the body spanwise (\f$y\f$) axis through the hinge, composed as \f$R = R_x(\phi) R_y(\theta)\f$. With \f$\tau = t - t_0\f$ and amplitude envelope \f$A(\tau)\f$ (0 before onset, raised cosine over `kin_ramp`, then 1): \f$\phi = A \phi_0 \sin(2\pi f \tau)\f$, \f$\theta = \theta_m + A \theta_0 \sin(2\pi f \tau + \psi)\f$. The centroid follows \f$x_c = x_h + R\,\mathbf{r}_\mathrm{off}\f$ and the ghost-cell velocities use the lab-frame angular velocity \f$\dot\phi \mathbf{e}_x + \dot\theta R_x(\phi)\mathbf{e}_y\f$. Set the initial `x[y,z]_centroid` and `angles` consistently with \f$t = 0\f$ so pre-process marks the body in the right place. - `kin_model = 2` is the smoothed linear pitch-ramp-and-hold of the AIAA low-Reynolds-number canonical cases (Eldredge et al. 2009, Ol et al. 2010) about the hinge, with no roll: \f$\theta(t) = \theta_m + \frac{\theta_0}{2}\left[1 + \frac{1}{a t_p}\log\frac{\cosh(a\tau)}{\cosh(a(\tau - t_p))}\right]\f$, \f$\tau = t - t_0\f$, \f$t_p = \theta_0/\Omega\f$, so the angle rises from `kin_theta_mean` by `kin_theta0` at nominal rate `kin_pitch_rate` starting at `kin_t0`, smoothed by `kin_smooth`. The same hinge, offset and centroid conventions as `kin_model = 1` apply. + +- Airfoils and STL models (`geometry` 4, 5, 11, 12) can be driven by `kin_model` too. For these the solver moves the working centroid to the centre of mass of the marked cells and keeps the difference as an internal offset \f$\mathbf{d}\f$; `kin_offset` still refers to the centroid the case file names, and the kinematics rotate `kin_offset` \f$- \mathbf{d}\f$ about the hinge. Name a centroid that lies inside the body: the rank that owns a patch is chosen by its centroid, so an STL placed with `model_translate` while its centroid stays at the origin is handed to ranks that do not hold it. Available variables: `x` (`x_cc(i)`), `y` (`y_cc(j)`), `z` (`z_cc(k)`), `t` (current simulation time), and `r` (the IB patch radius). The same intrinsic functions and `pi` constant apply; bare `e` is not available. diff --git a/src/simulation/m_compute_levelset.fpp b/src/simulation/m_compute_levelset.fpp index 61fe0e6182..84ad68398e 100644 --- a/src/simulation/m_compute_levelset.fpp +++ b/src/simulation/m_compute_levelset.fpp @@ -568,6 +568,8 @@ contains xyz_local(3) = z_cc(k) - center(3) end if xyz_local = matmul(inverse_rotation, xyz_local) + ! match the marker test in s_apply_ib_patches, which subtracts the offset too + xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset ! 3D models if (p > 0) then diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index a1f3425804..3a45c7d0cf 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1408,7 +1408,7 @@ contains end subroutine s_close_ib_force_history - !> @brief Writes IB state records to restart_data/ib_state.dat. Must be called only on rank 0. + !> @brief Writes IB state records to restart_data/ib_state.dat. Called on every rank. impure subroutine s_write_ib_state_file(time_step) integer, intent(in) :: time_step @@ -1420,9 +1420,41 @@ contains else call s_write_serial_ib_state(time_step) end if + call s_write_centroid_offsets(time_step) end subroutine s_write_ib_state_file + !> Write each global patch's centroid_offset to restart_data/ib_offset_.dat, so a restart continues about the same centre + !! of mass instead of re-measuring it from the body voxelised at the restart attitude. Collective. + impure subroutine s_write_centroid_offsets(step) + + integer, intent(in) :: step + character(len=path_len + 2*name_len) :: file_loc + real(wp), dimension(:,:), allocatable :: off_loc, off_glb + integer :: gid, i, k, file_unit + + allocate (off_loc(num_gbl_ibs, 3), off_glb(num_gbl_ibs, 3)) + off_loc = -huge(1._wp) ! ranks not holding a patch lose the max-reduction + do i = 1, num_ibs + off_loc(patch_ib(i)%gbl_patch_id,:) = patch_ib(i)%centroid_offset + end do + do gid = 1, num_gbl_ibs + do k = 1, 3 + call s_mpi_allreduce_max(off_loc(gid, k), off_glb(gid, k)) + end do + end do + if (proc_rank == 0) then + write (file_loc, '(A,I0,A)') trim(case_dir) // '/restart_data/ib_offset_', step, '.dat' + open (newunit=file_unit, file=trim(file_loc), status='replace', action='write') + do gid = 1, num_gbl_ibs + if (off_glb(gid, 1) > -huge(1._wp)) write (file_unit, '(I0,3(1X,ES24.16))') gid, off_glb(gid,:) + end do + close (file_unit) + end if + deallocate (off_loc, off_glb) + + end subroutine s_write_centroid_offsets + !> Write flow probe data at the current time step impure subroutine s_write_probe_files(t_step, q_cons_vf, accel_mag) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 483c065712..00711a1a85 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -79,7 +79,7 @@ contains !> Initializes the values of various IBM variables, such as ghost points and image points. impure subroutine s_ibm_setup() - integer :: i, j, k + integer :: i, j, k, gid real(wp) :: t_init !< initial time for prescribed kinematics integer(kind=8) :: max_num_gps @@ -124,8 +124,13 @@ contains $:GPU_UPDATE(device='[ib_markers%sf, corrected_gps%sf]') call s_apply_ib_patches(ib_markers) $:GPU_UPDATE(host='[ib_markers%sf]') + ! Loop over global ids: the offset is a collective reduction, and ranks hold different patches + do gid = 1, num_gbl_ibs + call s_get_neighborhood_idx(gid, i) + call s_compute_centroid_offset(gid, i) + end do + call s_restore_centroid_offsets(t_init) do i = 1, num_ibs - if (patch_ib(i)%moving_ibm /= 0) call s_compute_centroid_offset(i) ! offsets are computed after IB markers are generated $:GPU_UPDATE(device='[patch_ib(i)]') end do @@ -1137,8 +1142,9 @@ contains thetad = damp*patch_ib(i)%kin_theta0*sin(arg_th) + amp*patch_ib(i)%kin_theta0*omega*cos(arg_th) end if - ! centroid relative to the hinge: Rx(phi) Ry(theta) offset - r = patch_ib(i)%kin_offset + ! centroid relative to the hinge: Rx(phi) Ry(theta) offset. kin_offset points to the case-file centroid; a + ! centroid moved to the centre of mass (geometries 4, 5, 11, 12) sits centroid_offset short of it + r = patch_ib(i)%kin_offset - patch_ib(i)%centroid_offset c(1) = cos(theta)*r(1) + sin(theta)*r(3) c(2) = r(2) c(3) = -sin(theta)*r(1) + cos(theta)*r(3) @@ -1311,27 +1317,38 @@ contains !> 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_local, decoded_gbl_id - integer(kind=8) :: num_cells + !> Move a moving airfoil/STL patch's centroid to its centre of mass, keeping the difference in centroid_offset. Collective: + !! every rank calls it for every global id, contributing zeros for patches it does not hold. + subroutine s_compute_centroid_offset(gid, ib_marker) + + integer, intent(in) :: gid !< global patch id + integer, intent(in) :: ib_marker !< local index on this rank; <= 0 if not held + integer :: i, j, k, num_cells_local, decoded_gbl_id, geom + integer(kind=8) :: num_cells, needs_loc, needs_glb 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 + needs_loc = 0_8 + if (ib_marker > 0) then + geom = patch_ib(ib_marker)%geometry + if (patch_ib(ib_marker)%moving_ibm /= 0 .and. (geom == 4 .or. geom == 5 .or. geom == 11 .or. geom == 12)) & + & needs_loc = 1_8 + end if + call s_mpi_allreduce_integer_sum(needs_loc, needs_glb) + if (needs_glb == 0_8) then + if (ib_marker > 0) patch_ib(ib_marker)%centroid_offset(:) = 0._wp + return + end if + center_of_mass_local = [0._wp, 0._wp, 0._wp] + num_cells_local = 0 + if (ib_marker > 0) then ! 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) /= 0) then call s_decode_patch_periodicity(ib_markers%sf(i, j, k), decoded_gbl_id) - if (decoded_gbl_id == patch_ib(ib_marker)%gbl_patch_id) then + if (decoded_gbl_id == gid) 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) @@ -1340,36 +1357,65 @@ contains end do end do end do + end if - ! reduce the mass contribution over all MPI ranks and compute COM - call s_mpi_allreduce_integer_sum(int(num_cells_local, 8), num_cells) - if (num_cells /= 0) then - 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)) - center_of_mass = center_of_mass/real(num_cells, wp) - else - patch_ib(ib_marker)%centroid_offset = [0._wp, 0._wp, 0._wp] - return - end if - - ! 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] + ! reduce the mass contribution over all MPI ranks and compute COM + call s_mpi_allreduce_integer_sum(int(num_cells_local, 8), num_cells) + 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)) + if (ib_marker <= 0) return + if (num_cells == 0) then + 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) + + ! 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) end subroutine s_compute_centroid_offset + !> On restart, replace the re-measured centroid offsets with the ones the run was using (restart_data/ib_offset_.dat, + !! written with each checkpoint), and re-place kinematics-driven bodies about them. No file: the measured offsets stand. + impure subroutine s_restore_centroid_offsets(t_init) + + real(wp), intent(in) :: t_init + character(len=path_len + 2*name_len) :: file_loc + logical :: file_exist + integer :: gid, i, ios, file_unit, step + real(wp), dimension(3) :: off + + step = t_step_start + if (cfl_dt) step = n_start + if (step == 0) return + write (file_loc, '(A,I0,A)') trim(case_dir) // '/restart_data/ib_offset_', step, '.dat' + inquire (file=trim(file_loc), exist=file_exist) + if (.not. file_exist) return + open (newunit=file_unit, file=trim(file_loc), status='old', action='read', iostat=ios) + if (ios /= 0) return + do + read (file_unit, *, iostat=ios) gid, off + if (ios /= 0) exit + call s_get_neighborhood_idx(gid, i) + if (i > 0) then + patch_ib(i)%centroid_offset = off + if (patch_ib(i)%moving_ibm /= 0 .and. patch_ib(i)%kin_model > 0) call s_prescribed_kinematics(i, t_init) + end if + end do + close (file_unit) + + end subroutine s_restore_centroid_offsets + !> Computes the moment of inertia for an immersed boundary subroutine s_compute_moment_of_inertia(patch, axis, moment) diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index e2f5718a0c..001e425913 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -803,13 +803,6 @@ def check_ibm(self): self.prohibit(kin_model not in (0, 1, 2), f"patch_ib({i})%kin_model must be 0, 1 or 2") self.prohibit(kin_model > 0 and self.get(f"patch_ib({i})%moving_ibm", 0) != 1, f"patch_ib({i})%kin_model requires moving_ibm = 1") self.prohibit(kin_model > 0 and p <= 0, f"patch_ib({i})%kin_model requires a 3D case (p > 0)") - # Geometries 4, 5, 11 and 12 have their centroid replaced by the marked-cell centre of mass, with the - # difference kept in centroid_offset and re-applied when the patch is drawn. Prescribed kinematics write - # the centroid outright every stage, so the body would render centroid_offset away from the hinge. - self.prohibit( - kin_model > 0 and self.get(f"patch_ib({i})%geometry", 0) in (4, 5, 11, 12), - f"patch_ib({i})%kin_model is not supported for geometries 4, 5, 11 and 12, whose centroid is offset to the centre of mass", - ) self.prohibit(kin_model == 1 and (self.get(f"patch_ib({i})%kin_freq", 0) or 0) <= 0, f"patch_ib({i})%kin_freq must be > 0 when kin_model = 1") self.prohibit(kin_model == 2 and (self.get(f"patch_ib({i})%kin_pitch_rate", 0) or 0) <= 0, f"patch_ib({i})%kin_pitch_rate must be > 0 when kin_model = 2") self.prohibit(kin_model == 2 and (self.get(f"patch_ib({i})%kin_smooth", 0) or 0) <= 0, f"patch_ib({i})%kin_smooth must be > 0 when kin_model = 2") From 42f38f1aebbcce0ffe86c66b381c23923f645b9b Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 30 Sep 2026 09:21:15 -0400 Subject: [PATCH 2/2] Skip centroid-offset work when no patch needs an offset Only moving airfoils and STLs carry a centroid_offset, but setup ran one collective per global patch and each checkpoint three, so particle clouds (thousands of patches) paid for work they never use. One reduction now decides whether any patch needs it; otherwise nothing runs and no ib_offset file is written. The checkpoint writer also reduces in a single call, with only each patch's owner contributing. Co-Authored-By: Claude --- src/simulation/m_data_output.fpp | 27 ++++++++++++++----------- src/simulation/m_ibm.fpp | 34 ++++++++++++++++++++++---------- 2 files changed, 39 insertions(+), 22 deletions(-) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 3a45c7d0cf..63e8ac9615 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -1431,23 +1431,26 @@ contains integer, intent(in) :: step character(len=path_len + 2*name_len) :: file_loc real(wp), dimension(:,:), allocatable :: off_loc, off_glb - integer :: gid, i, k, file_unit - - allocate (off_loc(num_gbl_ibs, 3), off_glb(num_gbl_ibs, 3)) - off_loc = -huge(1._wp) ! ranks not holding a patch lose the max-reduction - do i = 1, num_ibs - off_loc(patch_ib(i)%gbl_patch_id,:) = patch_ib(i)%centroid_offset - end do - do gid = 1, num_gbl_ibs - do k = 1, 3 - call s_mpi_allreduce_max(off_loc(gid, k), off_glb(gid, k)) - end do + integer :: gid, i, ib_idx, n_own, file_unit + + if (.not. centroid_offsets_active) return + ! column 4 flags a patch that has an offset; only the owning rank contributes, so the sum is the value + allocate (off_loc(num_gbl_ibs, 4), off_glb(num_gbl_ibs, 4)) + off_loc = 0._wp + n_own = num_local_ibs + if (num_procs == 1) n_own = num_ibs + do i = 1, n_own + ib_idx = i + if (num_procs > 1) ib_idx = local_ib_patch_ids(i) + if (.not. f_needs_centroid_offset(patch_ib(ib_idx))) cycle + off_loc(patch_ib(ib_idx)%gbl_patch_id,:) = [patch_ib(ib_idx)%centroid_offset, 1._wp] end do + call s_mpi_allreduce_vectors_sum(off_loc, off_glb, num_gbl_ibs, 4) if (proc_rank == 0) then write (file_loc, '(A,I0,A)') trim(case_dir) // '/restart_data/ib_offset_', step, '.dat' open (newunit=file_unit, file=trim(file_loc), status='replace', action='write') do gid = 1, num_gbl_ibs - if (off_glb(gid, 1) > -huge(1._wp)) write (file_unit, '(I0,3(1X,ES24.16))') gid, off_glb(gid,:) + if (off_glb(gid, 4) > 0.5_wp) write (file_unit, '(I0,3(1X,ES24.16))') gid, off_glb(gid,1:3) end do close (file_unit) end if diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 00711a1a85..870fbdfaa8 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -48,6 +48,7 @@ module m_ibm $:GPU_DECLARE(create='[num_gps]') #endif logical :: moving_immersed_boundary_flag + logical :: centroid_offsets_active = .false. !< some patch is a moving airfoil or STL, which carries a centroid_offset ! IB MPI buffers integer, allocatable :: send_ids(:), recv_ids(:) @@ -80,6 +81,7 @@ contains impure subroutine s_ibm_setup() integer :: i, j, k, gid + integer(kind=8) :: n_need real(wp) :: t_init !< initial time for prescribed kinematics integer(kind=8) :: max_num_gps @@ -124,12 +126,17 @@ contains $:GPU_UPDATE(device='[ib_markers%sf, corrected_gps%sf]') call s_apply_ib_patches(ib_markers) $:GPU_UPDATE(host='[ib_markers%sf]') - ! Loop over global ids: the offset is a collective reduction, and ranks hold different patches - do gid = 1, num_gbl_ibs - call s_get_neighborhood_idx(gid, i) - call s_compute_centroid_offset(gid, i) - end do - call s_restore_centroid_offsets(t_init) + ! One reduction decides whether any patch needs an offset, so cases without one (particle clouds have thousands of + ! patches) skip the per-patch collectives. Then loop over global ids: ranks hold different patches. + call s_mpi_allreduce_integer_sum(int(count([(f_needs_centroid_offset(patch_ib(i)), i=1, num_ibs)]), 8), n_need) + centroid_offsets_active = n_need > 0_8 + if (centroid_offsets_active) then + do gid = 1, num_gbl_ibs + call s_get_neighborhood_idx(gid, i) + call s_compute_centroid_offset(gid, i) + end do + call s_restore_centroid_offsets(t_init) + end if do i = 1, num_ibs $:GPU_UPDATE(device='[patch_ib(i)]') end do @@ -1323,15 +1330,13 @@ contains integer, intent(in) :: gid !< global patch id integer, intent(in) :: ib_marker !< local index on this rank; <= 0 if not held - integer :: i, j, k, num_cells_local, decoded_gbl_id, geom + integer :: i, j, k, num_cells_local, decoded_gbl_id integer(kind=8) :: num_cells, needs_loc, needs_glb real(wp), dimension(1:3) :: center_of_mass, center_of_mass_local needs_loc = 0_8 if (ib_marker > 0) then - geom = patch_ib(ib_marker)%geometry - if (patch_ib(ib_marker)%moving_ibm /= 0 .and. (geom == 4 .or. geom == 5 .or. geom == 11 .or. geom == 12)) & - & needs_loc = 1_8 + if (f_needs_centroid_offset(patch_ib(ib_marker))) needs_loc = 1_8 end if call s_mpi_allreduce_integer_sum(needs_loc, needs_glb) if (needs_glb == 0_8) then @@ -1385,6 +1390,15 @@ contains end subroutine s_compute_centroid_offset + !> A moving airfoil or STL (geometries 4, 5, 11, 12) moves its centroid to the centre of mass and keeps a centroid_offset + pure logical function f_needs_centroid_offset(patch) + + type(ib_patch_parameters), intent(in) :: patch + + f_needs_centroid_offset = patch%moving_ibm /= 0 .and. any(patch%geometry == [4, 5, 11, 12]) + + end function f_needs_centroid_offset + !> On restart, replace the re-measured centroid offsets with the ones the run was using (restart_data/ib_offset_.dat, !! written with each checkpoint), and re-place kinematics-driven bodies about them. No file: the measured offsets stand. impure subroutine s_restore_centroid_offsets(t_init)