Let prescribed kinematics drive airfoil and STL immersed boundaries

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.
 

src/simulation/m_compute_levelset.fpp

@@ -572,6 +572,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

src/simulation/m_data_output.fpp

@@ -1414,7 +1414,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
@@ -1426,9 +1426,44 @@ 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_<step>.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, 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, 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
+        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)
 

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(:)
@@ -79,7 +80,8 @@ 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
+        integer(kind=8) :: n_need
         real(wp)        :: t_init  !< initial time for prescribed kinematics
         real(wp)        :: max_num_gps_rank
         integer(kind=8) :: max_num_gps
@@ -126,8 +128,18 @@ contains
         call s_apply_ib_patches(ib_markers)
         $:GPU_UPDATE(host='[ib_markers%sf]')
         call s_check_every_patch_marked()
+        ! 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
-            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 +1149,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)
@@ -1318,27 +1331,36 @@ 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)
+    !> 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)      :: 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
-        integer(kind=8)          :: num_cells
+        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
+            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
+            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)
@@ -1347,36 +1369,74 @@ 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
 
+    !> 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_<step>.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)
 

toolchain/mfc/case_validator.py

@@ -808,13 +808,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")