diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 219bddb1f..c993d49ed 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 d6f1d1b8b..b0a485207 100644 --- a/src/simulation/m_compute_levelset.fpp +++ b/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 diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 7902ff5ab..bd50f071d 100644 --- a/src/simulation/m_data_output.fpp +++ b/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_.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) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index d31c4addb..dd4111a2f 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(:) @@ -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_.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 304850b47..6371da5ec 100644 --- a/toolchain/mfc/case_validator.py +++ b/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")