diff --git a/docs/documentation/contributing.md b/docs/documentation/contributing.md index dd3c0f33c..366ea1d84 100644 --- a/docs/documentation/contributing.md +++ b/docs/documentation/contributing.md @@ -218,7 +218,7 @@ Both human reviewers and AI code reviewers reference this section. ### Parameter Plumbing - **Derived-type parameters are not auto-broadcast.** `generated_bcast.fpp` covers namelist *scalars* only. Each derived type (`chem_params`, `lag_params`, `rburn`) needs a hand-written `_emit_` in `toolchain/mfc/params/generators/fortran_gen.py` plus its call site in that generator's simulation branch, and, if it is read on device, an explicit ``$:GPU_UPDATE(device='[name]')`` in both the target's `m_global_parameters.fpp` and `src/simulation/m_start_up.fpp` — `GPU_DECLARE` alone does not make it device-resident. Regrouping existing scalars into a derived type silently drops their broadcast, leaving every non-root rank holding the `dflt_real` sentinel. Single-rank golden files cannot catch this, so pair such a change with a `ppn=2` test and confirm it fails without the emitter. -- **A `patch_ib` member that immersed-boundary ghost-point code reads must also be set in `s_add_cloud_particle`** (`src/simulation/m_particle_cloud.fpp`). `particle_cloud_ibs` is allocated without default initialization, and `s_reduce_ib_patch_array` copies the whole struct into `patch_ib`, overwriting the defaults assigned in `s_assign_default_values_to_user_inputs`. Anything left unset reaches the solver as uninitialized memory, and only where the allocation is not already zero-filled. A platform-only NaN is the signature of this class: a garbage `v_blow` once failed an AMD lane with `ICFL is NaN` while every NVIDIA lane and all local runs passed. +- **A `patch_ib` member that immersed-boundary ghost-point code reads must also be set for particle-cloud IBs in `s_assign_particle_cloud_ib_defaults`** (`src/simulation/m_start_up.fpp`). Pre-process writes only position, kinematics and radius to the IB state file; simulation builds every other property there, writing into a reused `patch_ib` slot. Anything it leaves unset keeps whatever that slot held, which may be a namelist patch's value or uninitialized memory, and shows up only where that memory is not already zero-filled. A platform-only NaN is the signature of this class: a garbage `v_blow` once failed an AMD lane with `ICFL is NaN` while every NVIDIA lane and all local runs passed. - **Runtime checks go where they run.** Shared constraints belong in `src/common/m_checker_common.fpp`, simulation-only ones in `src/simulation/m_checker.fpp`, and pre- and post-process ones in their own `m_checker.fpp`. Those two `s_check_inputs` are currently empty; that is still the correct home for their checks, not `m_checker_common`. - **Analytic initial conditions are compiled into the binary** and their expressions are AST-validated at case load, so syntax errors and unknown variables surface immediately and by name. Each IC variable maps to an `eqn_idx` expression in `QPVF_IDX_VARS` (`toolchain/mfc/case.py`); adding a patch-settable conserved variable means updating that map and the Fortran `eqn_idx` builder together, because a mismatch is a silent wrong index. - **Under `--case-optimization` the baked-in constants are dropped from the namelist**, so changing one requires a rebuild rather than a case-file edit. diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp index 4a456d661..17a986d6a 100644 --- a/src/common/m_helper.fpp +++ b/src/common/m_helper.fpp @@ -10,6 +10,7 @@ module m_helper use m_derived_types use m_global_parameters + use m_constants, only: BC_PERIODIC use ieee_arithmetic !< For checking NaN implicit none @@ -18,7 +19,8 @@ module m_helper public :: s_comp_n_from_prim, s_comp_n_from_cons, s_initialize_bubbles_model, s_initialize_nonpoly, s_simpson, s_transcoeff, & & s_int_to_str, s_transform_vec, s_transform_triangle, s_transform_model, s_swap, f_cross, f_create_transform_matrix, & & f_create_bbox, s_print_2D_array, f_xor, f_logical_to_int, associated_legendre, real_ylm, double_factorial, factorial, & - & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul + & f_cut_on, f_cut_off, s_downsample_data, s_upsample_data, s_cross_product, f_unit_vector, s_prng, modmul, & + & f_local_rank_owns_location contains @@ -699,4 +701,45 @@ contains end subroutine s_upsample_data + !> @brief True if `location` falls within this rank's own subdomain (a strict partition - each location is owned by exactly one + !! rank, unlike an overlapping multi-rank ghost-stencil neighborhood). Used by both pre_process (to decide which + !! generated/namelist IBs a rank writes to its own restart_data/ib_state_0.dat chunk) and simulation (to decide which IBs a rank + !! owns for its own periodic IB-state writes) so the two stay consistent. glb_bounds_in is the global domain extent, used only + !! to project a location that falls just outside the domain (floating-point edge case) onto the domain so some rank still claims + !! it. + function f_local_rank_owns_location(location, glb_bounds_in) result(owns_location) + + $:GPU_ROUTINE(parallelism='[seq]') + + real(wp), dimension(3), intent(in) :: location + type(bounds_info), dimension(3), intent(in) :: glb_bounds_in + logical :: owns_location + real(wp), dimension(3) :: projected_location + + owns_location = .true. + +#ifdef MFC_MPI + if (num_procs > 1) then + projected_location(:) = location(:) + + ! catch the edge case where the location lies just outside the computational domain + #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] + if (num_dims >= ${ID}$) then + if (ib_bc_${X}$%beg /= BC_PERIODIC) then + ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it + if (location(${ID}$) < glb_bounds_in(${ID}$)%beg) then + projected_location(${ID}$) = glb_bounds_in(${ID}$)%beg + else if (glb_bounds_in(${ID}$)%end < location(${ID}$)) then + projected_location(${ID}$) = glb_bounds_in(${ID}$)%end - 1.0e-10_wp + end if + end if + owns_location = owns_location .and. ${X}$_cb(-1) <= projected_location(${ID}$) & + & .and. projected_location(${ID}$) < ${X}$_cb(${DIM}$) + end if + #:endfor + end if +#endif + + end function f_local_rank_owns_location + end module m_helper diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 4c3a6712d..05bed7f79 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/src/post_process/m_global_parameters.fpp @@ -93,6 +93,7 @@ module m_global_parameters !> @name Boundary conditions in the x-, y- and z-coordinate directions !> @{ type(int_bounds_info) :: bc_x, bc_y, bc_z + type(int_bounds_info) :: ib_bc_x, ib_bc_y, ib_bc_z !< bc_x/y/z before decomposition overwrites them with MPI neighbor ranks type(bc_xyz_info) :: bc !> @} @@ -364,7 +365,7 @@ contains ! Particle clouds expand into individual IB patches at simulation startup, so num_ibs as read ! from the case file counts only the namelist patches. Match the global count the simulation - ! arrives at (s_reduce_ib_patch_array) so the IB state records can be read back. + ! arrives at (s_read_ib_restart_data) so the IB state records can be read back. do i = 1, num_particle_clouds num_ibs = num_ibs + particle_cloud(i)%num_particles end do diff --git a/src/post_process/m_start_up.fpp b/src/post_process/m_start_up.fpp index 86a00c001..c97133fff 100644 --- a/src/post_process/m_start_up.fpp +++ b/src/post_process/m_start_up.fpp @@ -952,6 +952,12 @@ contains end if call s_mpi_bcast_user_inputs() + + ! Save original BCs before decomposition overwrites them with MPI neighbor ranks + ib_bc_x = bc_x + ib_bc_y = bc_y + ib_bc_z = bc_z + call s_initialize_parallel_io() output_offsets = (/offset_x, offset_y, offset_z/) call s_mpi_decompose_computational_domain(write_silo_ghost_offsets=format == format_silo, adjust_local_domains=.false., & diff --git a/src/pre_process/m_data_output.fpp b/src/pre_process/m_data_output.fpp index e6b4847a0..85e9300ed 100644 --- a/src/pre_process/m_data_output.fpp +++ b/src/pre_process/m_data_output.fpp @@ -28,7 +28,7 @@ module m_data_output private public :: s_write_serial_data_files, s_write_parallel_data_files, s_write_data_files, s_initialize_data_output_module, & - & s_finalize_data_output_module + & s_finalize_data_output_module, s_write_ib_state_0_file type(scalar_field), allocatable, dimension(:) :: q_cons_temp @@ -733,6 +733,75 @@ contains end subroutine s_initialize_data_output_module + !> @brief Writes the initial IB layout (namelist patch_ib entries, then generated particle-cloud beds) that simulation reads + !! back at startup (s_read_ib_restart_data, src/simulation/m_start_up.fpp), in the layouts simulation's own IB state writers + !! use. Under file_per_process each rank writes only the IBs it owns to restart_data/lustre_0/ib_state_0_.dat; otherwise + !! only rank 0 calls this and writes every IB, in global-id order, to restart_data/ib_state_0.dat. + impure subroutine s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + + type(bounds_info), dimension(3), intent(in) :: glb_bounds + type(ib_patch_parameters), dimension(:), intent(in) :: particle_cloud_ibs + integer, intent(in) :: num_particle_cloud_ibs + character(LEN=len_trim(case_dir) + 2*name_len) :: file_loc + integer :: i, ios, file_unit + logical, dimension(num_ibs) :: owned + real(wp), dimension(3) :: centroid + + if (file_per_process) then + do i = 1, num_ibs + centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, patch_ib(i)%z_centroid] + owned(i) = f_local_rank_owns_location(centroid, glb_bounds) + end do + + if (proc_rank == 0) call s_create_directory(trim(case_dir) // '/restart_data/lustre_0') + call s_mpi_barrier() + call s_delay_file_access(proc_rank) + write (file_loc, '(A,i7.7,A)') '/restart_data/lustre_0/ib_state_0_', proc_rank, '.dat' + else + owned = .true. + call s_create_directory(trim(case_dir) // '/restart_data') + file_loc = '/restart_data/ib_state_0.dat' + end if + file_loc = trim(case_dir) // trim(file_loc) + + open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='replace', iostat=ios) + if (ios /= 0) call s_mpi_abort('Cannot open IB state output file: ' // trim(file_loc)) + + if (file_per_process) write (file_unit) count(owned) + num_particle_cloud_ibs + do i = 1, num_ibs + if (owned(i)) call s_write_ib_state_record(patch_ib(i), i) + end do + do i = 1, num_particle_cloud_ibs + call s_write_ib_state_record(particle_cloud_ibs(i), particle_cloud_ibs(i)%gbl_patch_id) + end do + + close (file_unit) + + contains + + !> Writes one 20-field IB state record, prefixed by its global id under file_per_process. + subroutine s_write_ib_state_record(ib_patch, gbl_id) + + type(ib_patch_parameters), intent(in) :: ib_patch + integer, intent(in) :: gbl_id + real(wp), dimension(20) :: ib_buf + + ib_buf = 0._wp + ib_buf(8:10) = ib_patch%vel + ib_buf(11:13) = ib_patch%angular_vel + ib_buf(14:16) = ib_patch%angles + ib_buf(17) = ib_patch%x_centroid + ib_buf(18) = ib_patch%y_centroid + ib_buf(19) = ib_patch%z_centroid + ib_buf(20) = ib_patch%radius + + if (file_per_process) write (file_unit) gbl_id + write (file_unit) ib_buf + + end subroutine s_write_ib_state_record + + end subroutine s_write_ib_state_0_file + !> Resets s_write_data_files pointer impure subroutine s_finalize_data_output_module diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index 7d85c2ac9..551f1432e 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/src/pre_process/m_global_parameters.fpp @@ -52,7 +52,8 @@ module m_global_parameters ! Cell indices (InDices With BUFFer): includes buffer except in pre_process type(int_bounds_info) :: idwbuff(1:3) type(int_bounds_info) :: bc_x, bc_y, bc_z !< Boundary conditions in the x-, y- and z-coordinate directions - type(bc_xyz_info) :: bc !< Combined BC storage (used by the shared beta-buffer routines; pre-process-local) + type(int_bounds_info) :: ib_bc_x, ib_bc_y, ib_bc_z !< bc_x/y/z before decomposition overwrites them with MPI neighbor ranks + type(bc_xyz_info) :: bc !< Combined BC storage (used by the shared beta-buffer routines; pre-process-local) ! simplex_params: auto-generated in generated_decls.fpp ! shear_num/shear_indices/shear_BC_flip_*, bc: in m_global_parameters_common integer :: fd_order !< Finite-difference order for CoM/probe derivative approximations @@ -361,6 +362,27 @@ contains patch_ib(i)%rotation_matrix_inverse = patch_ib(i)%rotation_matrix end do + num_particle_clouds = 0 + do i = 1, num_particle_clouds_max + particle_cloud(i)%x_centroid = 0._wp + particle_cloud(i)%y_centroid = 0._wp + particle_cloud(i)%z_centroid = 0._wp + particle_cloud(i)%length_x = dflt_real + particle_cloud(i)%length_y = dflt_real + particle_cloud(i)%length_z = dflt_real + particle_cloud(i)%num_particles = 0 + particle_cloud(i)%radius = dflt_real + particle_cloud(i)%mass = dflt_real + particle_cloud(i)%min_spacing = 0._wp + particle_cloud(i)%shell_inner_radius = dflt_real + particle_cloud(i)%shell_outer_radius = dflt_real + particle_cloud(i)%moving_ibm = 0 + particle_cloud(i)%seed = 0 + particle_cloud(i)%cloud_geometry = 1 + particle_cloud(i)%packing_method = dflt_int + particle_cloud(i)%periodic = 0 + end do + do i = 1, num_ib_airfoils_max ib_airfoil(i)%c = dflt_real ib_airfoil(i)%p = dflt_real diff --git a/src/pre_process/m_mpi_proxy.fpp b/src/pre_process/m_mpi_proxy.fpp index 45c1a3b50..826af94eb 100644 --- a/src/pre_process/m_mpi_proxy.fpp +++ b/src/pre_process/m_mpi_proxy.fpp @@ -135,6 +135,20 @@ contains call MPI_BCAST(patch_ib(i)%slip, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) end do + ! manual: particle_cloud (runtime loop to num_particle_clouds; irregular member subset) + do i = 1, num_particle_clouds + #:for VAR in ['x_centroid', 'y_centroid', 'z_centroid', 'length_x', 'length_y', 'length_z', & + & 'radius', 'mass', 'min_spacing', 'shell_inner_radius', 'shell_outer_radius'] + call MPI_BCAST(particle_cloud(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) + #:endfor + call MPI_BCAST(particle_cloud(i)%num_particles, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%cloud_geometry, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%packing_method, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + call MPI_BCAST(particle_cloud(i)%periodic, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + end do + ! manual: ib_airfoil (kept manual alongside patch_ib) do i = 1, num_ib_airfoils_max #:for VAR in ['c', 'p', 't', 'm'] diff --git a/src/simulation/m_particle_cloud.fpp b/src/pre_process/m_particle_cloud.fpp similarity index 80% rename from src/simulation/m_particle_cloud.fpp rename to src/pre_process/m_particle_cloud.fpp index 6e0843352..fba6452ca 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/pre_process/m_particle_cloud.fpp @@ -1,18 +1,22 @@ !> !! @file m_particle_cloud.fpp !! @brief Generates particle beds: converts particle_cloud specifications into -!! individual sphere/circle particle_cloud_ibs entries before reduction. +!! individual sphere/circle particle_cloud_ibs entries, written to +!! restart_data/ib_state_0.dat for simulation to read at startup. #:include 'macros.fpp' !> @brief Generates particle beds by converting particle_cloud patch specifications into individual immersed boundary patches before -!! domain reduction. Each rank runs the same deterministic placement so no MPI broadcast of particle positions is needed. +!! writing them to the initial IB state file. Under file_per_process it runs on every rank: each rank computes the same +!! deterministic placement (so no MPI broadcast of particle positions is needed) and keeps only the particles +!! f_local_rank_owns_location says are its own, so every generated particle is written by exactly one rank. Otherwise only rank 0 +!! runs it and keeps every particle. module m_particle_cloud use m_global_parameters use m_constants use m_mpi_common - use m_collisions + use m_helper implicit none @@ -22,13 +26,14 @@ module m_particle_cloud contains - !> Generate all particle beds and fill particle_cloud_ibs. Called on all ranks before s_reduce_ib_patch_array. Each packing - !! method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / s_particle_cloud_rejection_pack) - !! and hands back only the entries that fall within this rank's IB neighborhood. Only the first num_particle_cloud_ibs of them - !! are actually written - callers must use that count, not size(particle_cloud_ibs), since the remainder of the array is left - !! uninitialized. - impure subroutine s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) + !> Generate all particle beds and fill particle_cloud_ibs, keeping only the particles this rank owns under file_per_process (see + !! module docs). Each packing method owns and allocates its own per-cloud working array (see s_particle_cloud_lattice / + !! s_particle_cloud_rejection_pack) and hands back only the entries this rank keeps. Only the first num_particle_cloud_ibs of + !! them are actually written - callers must use that count, not size(particle_cloud_ibs), since the remainder of the array is + !! left uninitialized. + impure subroutine s_generate_particle_clouds(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs integer, intent(out) :: num_particle_cloud_ibs type(ib_patch_parameters), allocatable :: cloud_ibs(:) @@ -58,9 +63,9 @@ contains ! packing is box-only - the hemisphere-shell + lattice combination is rejected in case_validator.py. select case (particle_cloud(cloud_idx)%packing_method) case (1) ! rejection (random) packing method - call s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + call s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) case (2) ! lattice packing method - call s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + call s_particle_cloud_lattice(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) case default call s_mpi_abort("Particle cloud packing method is not a known packing method of MFC. Exiting.") end select @@ -79,18 +84,19 @@ contains end subroutine s_generate_particle_clouds !> Rejection-samples particle centres into a box or hemisphere-shell region with a minimum centre-to-centre spacing. Rejection - !! sampling needs every placed particle tracked (regardless of which rank's neighborhood it falls in) to detect overlaps - !! deterministically, so cloud_ibs is allocated here to the cloud's full requested particle count and only pared down to this - !! rank's neighborhood afterwards, via s_reduce_particle_cloud_ibs. Only the per-candidate geometry sampling differs between box - !! and hemisphere shell; it is delegated to s_sample_cloud_candidate, and every other step (overlap rejection via the spatial - !! hash, acceptance, reduction) is geometry-independent. - subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + !! sampling needs every placed particle tracked (regardless of which rank owns it) to detect overlaps deterministically, so + !! cloud_ibs is allocated here to the cloud's full requested particle count and only pared down to this rank's own subdomain + !! afterwards, via s_reduce_particle_cloud_ibs. Only the per-candidate geometry sampling differs between box and hemisphere + !! shell; it is delegated to s_sample_cloud_candidate, and every other step (overlap rejection via the spatial hash, acceptance, + !! reduction) is geometry-independent. + subroutine s_particle_cloud_rejection_pack(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) integer, intent(in) :: cloud_idx integer, intent(inout) :: glbl_idx + type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs integer, intent(out) :: num_cloud_ibs - integer :: ib_idx, n_placed, geom, seed, alloc_stat + integer :: ib_idx, n_placed, seed, alloc_stat integer(8) :: n_attempts, max_attempts real(wp) :: min_dist, rx, ry, rz logical :: overlaps, reject, periodic_pack @@ -120,12 +126,6 @@ contains nz_bins = max(1, ceiling(length_z/min_dist)) if (num_dims < 3) nz_bins = 1 - if (num_dims < 3) then - geom = 2 ! circle for 2D - else - geom = 8 ! sphere for 3D - end if - max_attempts = int(particle_cloud(cloud_idx)%num_particles, 8)*1000_8 n_placed = 0 n_attempts = 0 @@ -164,7 +164,7 @@ contains hash_head(slot) = n_placed glbl_idx = glbl_idx + 1 - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, rx, ry, rz, cloud_ibs) + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, rx, ry, rz, cloud_ibs) end if end do @@ -174,7 +174,7 @@ contains deallocate (placed, hash_head, chain_next) - call s_reduce_particle_cloud_ibs(cloud_ibs, ib_idx) + if (file_per_process) call s_reduce_particle_cloud_ibs(cloud_ibs, glb_bounds, ib_idx) num_cloud_ibs = ib_idx end subroutine s_particle_cloud_rejection_pack @@ -243,16 +243,17 @@ contains !> Places particles on the optimally dense lattice for the cloud region: a triangular lattice in 2D, a face-centered cubic !! lattice in 3D. The lattice spacing is set by the particle density (num_particles over the region area/volume); if that !! spacing falls below the required centre-to-centre distance (2*radius + min_spacing), the region is too dense and the run is - !! aborted. No two lattice sites can overlap, so unlike rejection packing each site's IB neighborhood membership - !! (get_neighbor_bounds() must already have run) is checked as it is generated and only in-neighborhood sites are stored; - !! cloud_ibs is therefore allocated to the neighborhood-sized cap rather than the cloud's full particle count. - subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, cloud_ibs, num_cloud_ibs) + !! aborted. No two lattice sites can overlap, so unlike rejection packing each site's local ownership + !! (f_local_rank_owns_location) is checked as it is generated and only this rank's sites are stored; cloud_ibs is therefore + !! allocated to a worst-case cap rather than the cloud's full particle count. + subroutine s_particle_cloud_lattice(cloud_idx, glbl_idx, glb_bounds, cloud_ibs, num_cloud_ibs) integer, intent(in) :: cloud_idx integer, intent(inout) :: glbl_idx + type(bounds_info), dimension(3), intent(in) :: glb_bounds type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: cloud_ibs integer, intent(out) :: num_cloud_ibs - integer :: ib_idx, n_placed, n_target, geom + integer :: ib_idx, n_placed, n_target integer :: row, col, ncx, ncy, ix, jy, kz, b real(wp) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist real(wp) :: spacing, row_dy, cell, x0, px, py @@ -274,11 +275,9 @@ contains n_placed = 0 if (num_dims < 3) then - geom = 2 ! circle for 2D ! Triangular lattice: area per particle = (sqrt(3)/2)*spacing**2. spacing = sqrt(2._wp*(xmax - xmin)*(ymax - ymin)/(sqrt(3._wp)*real(n_target, wp))) else - geom = 8 ! sphere for 3D ! Face-centered cubic lattice: volume per particle = spacing**3/sqrt(2). spacing = (sqrt(2._wp)*(xmax - xmin)*(ymax - ymin)*(zmax - zmin)/real(n_target, wp))**(1._wp/3._wp) end if @@ -301,9 +300,8 @@ contains do while (px <= xmax .and. n_placed < n_target) glbl_idx = glbl_idx + 1 centroid = [px, py, particle_cloud(cloud_idx)%z_centroid] - if (f_neighborhood_ranks_own_location(centroid)) then - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), centroid(3), & - & cloud_ibs) + if (.not. file_per_process .or. f_local_rank_owns_location(centroid, glb_bounds)) then + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, centroid(1), centroid(2), centroid(3), cloud_ibs) end if n_placed = n_placed + 1 col = col + 1 @@ -328,9 +326,9 @@ contains centroid = [xmin + real(ix, wp)*cell + bx_off(b), ymin + real(jy, wp)*cell + by_off(b), & & zmin + real(kz, wp)*cell + bz_off(b)] glbl_idx = glbl_idx + 1 - if (f_neighborhood_ranks_own_location(centroid)) then - call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, centroid(1), centroid(2), & - & centroid(3), cloud_ibs) + if (.not. file_per_process .or. f_local_rank_owns_location(centroid, glb_bounds)) then + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, centroid(1), centroid(2), centroid(3), & + & cloud_ibs) end if n_placed = n_placed + 1 end do @@ -344,60 +342,52 @@ contains end subroutine s_particle_cloud_lattice - !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx. The caller decides whether - !! this particle belongs in the array (neighborhood membership, for lattice packing, or unconditionally for rejection packing - - !! see s_particle_cloud_lattice / s_particle_cloud_rejection_pack) and supplies its already-assigned, absolute global patch id - !! via glbl_idx - s_reduce_ib_patch_array copies gbl_patch_id as-is. Shared by all packing methods so the per-particle - !! ib_patch_parameters setup stays in one place. - subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, px, py, pz, particle_cloud_ibs) + !> Compacts cloud_ibs(1:num_cloud_ibs) in place, discarding entries this rank does not own (f_local_rank_owns_location) and + !! updating num_cloud_ibs to the retained count. Used by rejection packing, which cannot filter as it places particles (see + !! s_particle_cloud_rejection_pack), to pare its full, unfiltered placement down to this rank's own subdomain. + subroutine s_reduce_particle_cloud_ibs(cloud_ibs, glb_bounds, num_cloud_ibs) + + type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs + type(bounds_info), dimension(3), intent(in) :: glb_bounds + integer, intent(inout) :: num_cloud_ibs + integer :: i, write_idx + real(wp), dimension(3) :: centroid + + write_idx = 0 + do i = 1, num_cloud_ibs + centroid = [cloud_ibs(i)%x_centroid, cloud_ibs(i)%y_centroid, 0._wp] + if (num_dims == 3) centroid(3) = cloud_ibs(i)%z_centroid + if (f_local_rank_owns_location(centroid, glb_bounds)) then + write_idx = write_idx + 1 + if (write_idx /= i) cloud_ibs(write_idx) = cloud_ibs(i) + end if + end do + num_cloud_ibs = write_idx + + end subroutine s_reduce_particle_cloud_ibs + + !> Writes a single placed particle into particle_cloud_ibs at the next free slot, advancing ib_idx, tagged with its + !! already-assigned, absolute global patch id via glbl_idx. Only the fields s_write_ib_state_0_file writes are set; simulation + !! fills every other property in s_assign_particle_cloud_ib_defaults (src/simulation/m_start_up.fpp). + subroutine s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, px, py, pz, particle_cloud_ibs) - integer, intent(in) :: cloud_idx, glbl_idx, geom + integer, intent(in) :: cloud_idx, glbl_idx integer, intent(inout) :: ib_idx real(wp), intent(in) :: px, py, pz type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs ib_idx = ib_idx + 1 @:PROHIBIT(ib_idx > size(particle_cloud_ibs), & - & "Too many particle-cloud IBs in one rank's neighborhood. Modify case file or increase num_ib_patches_max_namelist.") + & "Too many particle-cloud IBs on one rank. Modify case file or increase num_ib_patches_max_namelist.") particle_cloud_ibs(ib_idx)%gbl_patch_id = glbl_idx - particle_cloud_ibs(ib_idx)%geometry = geom particle_cloud_ibs(ib_idx)%x_centroid = px particle_cloud_ibs(ib_idx)%y_centroid = py particle_cloud_ibs(ib_idx)%z_centroid = pz - particle_cloud_ibs(ib_idx)%step_x_centroid = px - particle_cloud_ibs(ib_idx)%step_y_centroid = py - particle_cloud_ibs(ib_idx)%step_z_centroid = pz - particle_cloud_ibs(ib_idx)%angles(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angles(:) = 0._wp particle_cloud_ibs(ib_idx)%vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_vel(:) = 0._wp particle_cloud_ibs(ib_idx)%angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%step_angular_vel(:) = 0._wp - particle_cloud_ibs(ib_idx)%force(:) = 0._wp - particle_cloud_ibs(ib_idx)%torque(:) = 0._wp - particle_cloud_ibs(ib_idx)%centroid_offset(:) = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix = 0._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(1, 1) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(2, 2) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix(3, 3) = 1._wp - particle_cloud_ibs(ib_idx)%rotation_matrix_inverse = particle_cloud_ibs(ib_idx)%rotation_matrix + particle_cloud_ibs(ib_idx)%angles(:) = 0._wp particle_cloud_ibs(ib_idx)%radius = particle_cloud(cloud_idx)%radius - particle_cloud_ibs(ib_idx)%mass = particle_cloud(cloud_idx)%mass - particle_cloud_ibs(ib_idx)%moment = dflt_real - particle_cloud_ibs(ib_idx)%moving_ibm = particle_cloud(cloud_idx)%moving_ibm - particle_cloud_ibs(ib_idx)%slip = .false. - - ! Particles are inert surfaces. These must be set explicitly: particle_cloud_ibs is - ! allocated (not default-initialized) and s_reduce_ib_patch_array copies the whole - ! struct into patch_ib, overwriting the defaults from - ! s_assign_default_values_to_user_inputs -- so anything left unset here reaches the - ! solver as uninitialized memory (a nonzero v_blow injects a garbage wall-normal - ! velocity and NaNs the field). - particle_cloud_ibs(ib_idx)%v_blow = 0._wp - particle_cloud_ibs(ib_idx)%inj_species = 0 - particle_cloud_ibs(ib_idx)%burn_rate_exp = 0._wp - particle_cloud_ibs(ib_idx)%burn_rate_pref = 0._wp end subroutine s_add_cloud_particle @@ -500,29 +490,6 @@ contains end subroutine s_check_cloud_particle_overlap - !> Compacts cloud_ibs(1:num_ibs) in place, discarding entries outside this rank's IB neighborhood (get_neighbor_bounds() must - !! already have run) and updating num_ibs to the retained count. Used by rejection packing, which cannot filter as it places - !! particles (see s_particle_cloud_rejection_pack), to pare its full, unfiltered placement down to this rank's neighborhood. - subroutine s_reduce_particle_cloud_ibs(cloud_ibs, num_cloud_ibs) - - type(ib_patch_parameters), intent(inout), dimension(:) :: cloud_ibs - integer, intent(inout) :: num_cloud_ibs - integer :: i, write_idx - real(wp), dimension(3) :: centroid - - write_idx = 0 - do i = 1, num_cloud_ibs - centroid = [cloud_ibs(i)%x_centroid, cloud_ibs(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = cloud_ibs(i)%z_centroid - if (f_neighborhood_ranks_own_location(centroid)) then - write_idx = write_idx + 1 - if (write_idx /= i) cloud_ibs(write_idx) = cloud_ibs(i) - end if - end do - num_cloud_ibs = write_idx - - end subroutine s_reduce_particle_cloud_ibs - !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1). function f_xorshift(seed) result(rval) diff --git a/src/pre_process/m_start_up.fpp b/src/pre_process/m_start_up.fpp index 04e399737..fd21f95d9 100644 --- a/src/pre_process/m_start_up.fpp +++ b/src/pre_process/m_start_up.fpp @@ -28,6 +28,7 @@ module m_start_up use m_check_patches use m_check_ib_patches + use m_particle_cloud use m_helper use m_checker_common use m_checker @@ -39,7 +40,8 @@ module m_start_up private public :: s_read_input_file, s_check_input_file, s_read_grid_data_files, s_read_ic_data_files, s_read_serial_grid_data_files, & & s_read_serial_ic_data_files, s_read_parallel_grid_data_files, s_read_parallel_ic_data_files, s_check_grid_data_files, & - & s_initialize_modules, s_initialize_mpi_domain, s_finalize_modules, s_apply_initial_condition, s_save_data, s_read_grid + & s_initialize_modules, s_initialize_mpi_domain, s_finalize_modules, s_apply_initial_condition, s_save_data, s_read_grid, & + & s_write_ib_state_0 abstract interface @@ -137,6 +139,27 @@ contains end subroutine s_check_input_file + !> @brief Generates the particle-cloud beds (if any) and writes the initial IB state file that simulation reads back at startup + !! (src/simulation/m_start_up.fpp:s_read_ib_restart_data). Must run after the domain is decomposed (s_initialize_mpi_domain) and + !! the grid is populated (s_read_grid). Under file_per_process every rank computes the same deterministic placement and keeps + !! only the IBs f_local_rank_owns_location says are its own; otherwise rank 0 alone generates and writes every IB. + impure subroutine s_write_ib_state_0() + + type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) + integer :: num_particle_cloud_ibs + type(bounds_info), dimension(3) :: glb_bounds + + if (.not. ib) return + if (.not. file_per_process .and. proc_rank /= 0) return + + glb_bounds = (/x_domain_glb, y_domain_glb, z_domain_glb/) + + call s_generate_particle_clouds(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + call s_write_ib_state_0_file(glb_bounds, particle_cloud_ibs, num_particle_cloud_ibs) + deallocate (particle_cloud_ibs) + + end subroutine s_write_ib_state_0 + !> The goal of this subroutine is to read in any preexisting grid data as well as based on the imported grid, complete the !! necessary global computational domain parameters. impure subroutine s_read_serial_grid_data_files @@ -623,6 +646,12 @@ contains ! Broadcasting the user inputs to all of the processors and performing the parallel computational domain decomposition. ! Neither procedure has to be carried out if pre-process is in fact not truly executed in parallel. call s_mpi_bcast_user_inputs() + + ! Save original BCs before decomposition overwrites them with MPI neighbor ranks + ib_bc_x = bc_x + ib_bc_y = bc_y + ib_bc_z = bc_z + call s_initialize_parallel_io() ! Save the global domain bounds before decomposition overwrites x/y/z_domain with each processor's local sub-domain bounds diff --git a/src/pre_process/p_main.f90 b/src/pre_process/p_main.f90 index e16ff6857..08cdc7576 100644 --- a/src/pre_process/p_main.f90 +++ b/src/pre_process/p_main.f90 @@ -24,6 +24,8 @@ program p_main call s_read_grid() + call s_write_ib_state_0() + allocate (proc_time(0:num_procs - 1)) call s_apply_initial_condition(start, finish) diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index 54ad9465d..c8148fcd1 100644 --- a/src/simulation/m_collisions.fpp +++ b/src/simulation/m_collisions.fpp @@ -21,7 +21,7 @@ module m_collisions implicit none private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module, & - & f_local_rank_owns_location, f_neighborhood_ranks_own_location, ib_gbl_idx_lookup, collisions_active + & f_neighborhood_ranks_own_location, ib_gbl_idx_lookup, collisions_active ! overlap distances for computing collisions integer, allocatable, dimension(:,:) :: collision_lookup real(wp), allocatable, dimension(:,:) :: wall_overlap_distances @@ -128,7 +128,7 @@ contains overlap_distance = patch_ib(pid1)%radius + patch_ib(pid2)%radius - norm2(normal_vector) if (overlap_distance > 0._wp) then ! if the two patches are close enough to collide normal_vector = normal_vector/norm2(normal_vector) - if (f_local_rank_owns_location(centroid_1)) then + if (f_local_rank_owns_location(centroid_1, glb_bounds)) then ! compute constants of the collision effective_mass = 1.0_wp/((1.0_wp/patch_ib(pid1)%mass) + (1._wp/(patch_ib(pid2)%mass))) k = spring_stiffness*effective_mass @@ -207,7 +207,7 @@ contains ! ensure the local rank owns that collision before proceeding collision_location = [patch_ib(patch_id)%x_centroid, patch_ib(patch_id)%y_centroid, 0._wp] if (num_dims == 3) collision_location(3) = patch_ib(patch_id)%z_centroid - if (f_local_rank_owns_location(collision_location)) then + if (f_local_rank_owns_location(collision_location, glb_bounds)) then k = spring_stiffness*patch_ib(patch_id)%mass eta = damping_parameter*patch_ib(patch_id)%mass @@ -386,41 +386,6 @@ contains end subroutine s_detect_wall_collisions - !> @brief function checks if this local MPI processor owns this specific collision - function f_local_rank_owns_location(location) result(owns_collision) - - $:GPU_ROUTINE(parallelism='[seq]') - - real(wp), dimension(3), intent(in) :: location - logical :: owns_collision - real(wp), dimension(3) :: projected_location - - owns_collision = .true. - -#ifdef MFC_MPI - if (num_procs > 1) then - projected_location(:) = location(:) - - ! catch the edge case where th collision lies just outside the computational domain - #:for X, ID, DIM in [('x', 1, 'm'), ('y', 2, 'n'), ('z', 3, 'p')] - if (num_dims >= ${ID}$) then - if (ib_bc_${X}$%beg /= BC_PERIODIC) then - ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it - if (location(${ID}$) < glb_bounds(${ID}$)%beg) then - projected_location(${ID}$) = glb_bounds(${ID}$)%beg - else if (glb_bounds(${ID}$)%end < location(${ID}$)) then - projected_location(${ID}$) = glb_bounds(${ID}$)%end - 1.0e-10_wp - end if - end if - owns_collision = owns_collision .and. ${X}$_cb(-1) <= projected_location(${ID}$) & - & .and. projected_location(${ID}$) < ${X}$_cb(${DIM}$) - end if - #:endfor - end if -#endif - - end function f_local_rank_owns_location - !> @brief function checks if this local MPI processor owns this specific collision function f_neighborhood_ranks_own_location(location) result(owns_collision) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 66690f25e..ce54d7550 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1445,7 +1445,7 @@ contains end if ! check if in local domain - if (f_local_rank_owns_location(centroid)) then + if (f_local_rank_owns_location(centroid, glb_bounds)) then local_output_idx = local_output_idx + 1 local_ib_patch_ids(local_output_idx) = output_idx end if diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 4a9a6d231..58b7d55a9 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -43,7 +43,6 @@ module m_start_up use m_ibm use m_ib_patches use m_model - use m_particle_cloud use m_collisions use m_compile_specific use m_checker_common @@ -913,33 +912,19 @@ contains if (model_eqns == model_eqns_6eq) call s_initialize_internal_energy_equations(q_cons_ts(1)%vf) if (ib) then - block - type(ib_patch_parameters), allocatable :: particle_cloud_ibs(:) - integer :: num_particle_cloud_ibs - - call s_instantiate_STL_models() - call s_initialize_ib_airfoils() - call s_get_neighbor_bounds() - - if (cfl_dt .and. n_start > 0) then - call s_read_ib_restart_data(n_start) - allocate (particle_cloud_ibs(0)) - num_particle_cloud_ibs = 0 - else if (t_step_start > 0) then - call s_read_ib_restart_data(t_step_start) - allocate (particle_cloud_ibs(0)) - num_particle_cloud_ibs = 0 - else - call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) - end if - call s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) - deallocate (particle_cloud_ibs) - end block - call s_ibm_setup() - if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) then - call s_write_ib_data_file(0) - call s_write_ib_state_file(0) + call s_instantiate_STL_models() + call s_initialize_ib_airfoils() + call s_get_neighbor_bounds() + if (cfl_dt .and. n_start > 0) then + call s_read_ib_restart_data(n_start) + else if (t_step_start > 0) then + call s_read_ib_restart_data(t_step_start) + else + call s_read_ib_restart_data(0) end if + call s_build_ib_neighborhood() + call s_ibm_setup() + if (t_step_start == 0 .or. (cfl_dt .and. n_start == 0)) call s_write_ib_data_file(0) end if if (bodyForces .or. synthetic_turbulence) call s_initialize_body_forces_module() if (acoustic_source) call s_precalculate_acoustic_spatial_sources() @@ -1160,200 +1145,190 @@ contains end subroutine s_finalize_modules - !> @brief Reads IB kinematic state from restart_data/ib_state.dat on restart. Rank 0 reads the last num_ibs records and - !! broadcasts to all ranks. Overwrites patch_ib vel, angular_vel, angles, and centroid. - impure subroutine s_read_ib_restart_data(t_step) + !> @brief Fills the properties of a generated particle-cloud IB that the IB state file does not carry (geometry, mass, + !! moving_ibm, inert surface, identity rotation matrix, zeroed step state). This is the only place they are set - pre_process + !! (s_add_cloud_particle) generates only position, kinematics and radius. cloud_ib_idx is the global patch id minus the number + !! of namelist patches; pre_process numbers particles cloud by cloud. + subroutine s_assign_particle_cloud_ib_defaults(cloud_ib_idx, ib_patch) - integer, intent(in) :: t_step - character(len=path_len + 2*name_len) :: file_loc - integer :: i, ios, file_unit, ierr - integer :: r, nlocal, gbl_id - integer, parameter :: NFIELDS_PER_IB = 20 - real(wp) :: ib_buf(NFIELDS_PER_IB) - logical :: file_exist - character(len=10) :: t_step_string + integer, intent(in) :: cloud_ib_idx + type(ib_patch_parameters), intent(inout) :: ib_patch + integer :: cloud_idx, idx_in_cloud - if (file_per_process) then - call s_int_to_str(t_step, t_step_string) + idx_in_cloud = cloud_ib_idx + do cloud_idx = 1, num_particle_clouds - 1 + if (idx_in_cloud <= particle_cloud(cloud_idx)%num_particles) exit + idx_in_cloud = idx_in_cloud - particle_cloud(cloud_idx)%num_particles + end do - do r = 0, num_procs - 1 - write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', r, '.dat' - file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc) + ib_patch%geometry = merge(2, 8, num_dims < 3) + ib_patch%step_x_centroid = 0._wp + ib_patch%step_y_centroid = 0._wp + ib_patch%step_z_centroid = 0._wp + ib_patch%step_angles(:) = 0._wp + ib_patch%step_vel(:) = 0._wp + ib_patch%step_angular_vel(:) = 0._wp + ib_patch%force(:) = 0._wp + ib_patch%torque(:) = 0._wp + ib_patch%centroid_offset(:) = 0._wp + ib_patch%rotation_matrix = 0._wp + ib_patch%rotation_matrix(1, 1) = 1._wp + ib_patch%rotation_matrix(2, 2) = 1._wp + ib_patch%rotation_matrix(3, 3) = 1._wp + ib_patch%rotation_matrix_inverse = ib_patch%rotation_matrix + ib_patch%mass = particle_cloud(cloud_idx)%mass + ib_patch%moment = dflt_real + ib_patch%moving_ibm = particle_cloud(cloud_idx)%moving_ibm + ib_patch%slip = .false. + ib_patch%v_blow = 0._wp + ib_patch%inj_species = 0 + ib_patch%burn_rate_exp = 0._wp + ib_patch%burn_rate_pref = 0._wp + + end subroutine s_assign_particle_cloud_ib_defaults + + !> @brief Loads the IBs this rank owns from the IB state file for t_step into patch_ib(1:num_ibs), all of which are local. Under + !! file_per_process the rank reads only its own restart_data/lustre_/ib_state__.dat, which holds exactly + !! its IBs; otherwise every rank reads every record of restart_data/ib_state_.dat and keeps the ones + !! f_local_rank_owns_location assigns it. Records carry kinematics, position and radius; every other property comes from the + !! namelist patch (global id <= num_ibs) or the particle cloud the IB was generated from. Written by pre_process at t_step = 0 + !! (src/pre_process/m_data_output.fpp:s_write_ib_state_0_file) and by s_write_ib_state_file on later steps. + impure subroutine s_read_ib_restart_data(t_step) - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (.not. file_exist) call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc)) - - open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios) - if (ios /= 0) call s_mpi_abort('Error opening IB state restart file: ' // trim(file_loc)) - - read (file_unit, iostat=ios) nlocal - if (ios /= 0) call s_mpi_abort('Error reading IB state file header: ' // trim(file_loc)) - - do i = 1, nlocal - read (file_unit, iostat=ios) gbl_id - if (ios /= 0) call s_mpi_abort('Error reading IB patch ID: ' // trim(file_loc)) - read (file_unit, iostat=ios) ib_buf - if (ios /= 0) call s_mpi_abort('Error reading IB state data: ' // trim(file_loc)) - - patch_ib(gbl_id)%vel = ib_buf(8:10) - patch_ib(gbl_id)%angular_vel = ib_buf(11:13) - patch_ib(gbl_id)%angles = ib_buf(14:16) - patch_ib(gbl_id)%x_centroid = ib_buf(17) - patch_ib(gbl_id)%y_centroid = ib_buf(18) - patch_ib(gbl_id)%z_centroid = ib_buf(19) - end do + integer, intent(in) :: t_step + type(ib_patch_parameters), allocatable, dimension(:) :: namelist_ibs + character(len=path_len + 2*name_len) :: file_loc + integer :: i, ios, file_unit, gbl_id, n_records + integer, parameter :: NFIELDS_PER_IB = 20 + real(wp) :: ib_buf(NFIELDS_PER_IB) + character(len=10) :: t_step_string - close (file_unit) - end do + moving_immersed_boundary_flag = any(patch_ib(1:num_ibs)%moving_ibm /= 0) & + & .or. any(particle_cloud(1:num_particle_clouds)%moving_ibm /= 0) + + allocate (namelist_ibs(num_ibs)) + namelist_ibs(:) = patch_ib(1:num_ibs) + num_gbl_ibs = num_ibs + sum(particle_cloud(1:num_particle_clouds)%num_particles) + + if (file_per_process) then + call s_int_to_str(t_step, t_step_string) + write (file_loc, '(A,I0,A,i7.7,A)') 'ib_state_', t_step, '_', proc_rank, '.dat' + file_loc = trim(case_dir) // '/restart_data/lustre_' // trim(t_step_string) // '/' // trim(file_loc) else write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat' file_loc = trim(case_dir) // trim(file_loc) + end if - if (proc_rank == 0) then - inquire (FILE=trim(file_loc), EXIST=file_exist) - if (.not. file_exist) then - call s_mpi_abort('Cannot open IB state file for restart: ' // trim(file_loc)) - end if + open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', action='read', iostat=ios) + if (ios /= 0) call s_mpi_abort('Error opening IB state file: ' // trim(file_loc)) - open (newunit=file_unit, file=trim(file_loc), form='unformatted', access='stream', status='old', iostat=ios) - if (ios /= 0) call s_mpi_abort('Error opening IB state restart file: ' // trim(file_loc)) + n_records = num_gbl_ibs + if (file_per_process) then + read (file_unit, iostat=ios) n_records + if (ios /= 0) call s_mpi_abort('Error reading IB state file header: ' // trim(file_loc)) + end if - do i = 1, num_ibs - read (file_unit, iostat=ios) ib_buf - if (ios /= 0) call s_mpi_abort('Error reading IB state restart file') + num_ibs = 0 + do i = 1, n_records + gbl_id = i + if (file_per_process) read (file_unit, iostat=ios) gbl_id + if (ios == 0) read (file_unit, iostat=ios) ib_buf + if (ios /= 0) call s_mpi_abort('Error reading IB state file: ' // trim(file_loc)) - patch_ib(i)%vel = ib_buf(8:10) - patch_ib(i)%angular_vel = ib_buf(11:13) - patch_ib(i)%angles = ib_buf(14:16) - patch_ib(i)%x_centroid = ib_buf(17) - patch_ib(i)%y_centroid = ib_buf(18) - patch_ib(i)%z_centroid = ib_buf(19) - end do + if (.not. file_per_process) then + if (.not. f_local_rank_owns_location(ib_buf(17:19), glb_bounds)) cycle + end if - close (file_unit) + num_ibs = num_ibs + 1 + @:PROHIBIT(num_ibs > num_local_ibs_max, & + & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") + if (gbl_id <= size(namelist_ibs)) then + patch_ib(num_ibs) = namelist_ibs(gbl_id) + else + call s_assign_particle_cloud_ib_defaults(gbl_id - size(namelist_ibs), patch_ib(num_ibs)) end if + patch_ib(num_ibs)%gbl_patch_id = gbl_id + patch_ib(num_ibs)%vel = ib_buf(8:10) + patch_ib(num_ibs)%angular_vel = ib_buf(11:13) + patch_ib(num_ibs)%angles = ib_buf(14:16) + patch_ib(num_ibs)%x_centroid = ib_buf(17) + patch_ib(num_ibs)%y_centroid = ib_buf(18) + patch_ib(num_ibs)%z_centroid = ib_buf(19) + patch_ib(num_ibs)%radius = ib_buf(20) + local_ib_patch_ids(num_ibs) = num_ibs + end do -#ifdef MFC_MPI - do i = 1, num_ibs - call MPI_BCAST(patch_ib(i)%vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%angular_vel, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%angles, 3, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%x_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%y_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - call MPI_BCAST(patch_ib(i)%z_centroid, 1, mpi_p, 0, MPI_COMM_WORLD, ierr) - end do -#endif - end if + close (file_unit) + deallocate (namelist_ibs) + + num_local_ibs = num_ibs end subroutine s_read_ib_restart_data - !> @brief Merges patch_ib (namelist patches, fixed at num_ib_patches_max_namelist) with particle_cloud_ibs (already filtered by - !! s_generate_particle_clouds to this rank's IB neighborhood, each entry already tagged with its final, absolute gbl_patch_id) - !! and reduces to only the patches in or near the local computational domain. patch_ib is never reallocated; the local subset is - !! written in-place from the front. particle_cloud_ibs is owned by the caller and freed there after this returns. - !! num_particle_cloud_ibs is the number of entries s_generate_particle_clouds actually wrote into particle_cloud_ibs - it may be - !! allocated to a larger worst-case capacity, so size() of it must never be used as the valid-entry count. - subroutine s_reduce_ib_patch_array(particle_cloud_ibs, num_particle_cloud_ibs) - - type(ib_patch_parameters), intent(in), dimension(:) :: particle_cloud_ibs - integer, intent(in) :: num_particle_cloud_ibs - real(wp), dimension(3) :: centroid - integer :: i - integer :: num_namelist_ibs, num_bed_ibs - - num_namelist_ibs = num_ibs - num_bed_ibs = 0 - do i = 1, num_particle_clouds - num_bed_ibs = num_bed_ibs + particle_cloud(i)%num_particles - end do + !> @brief Completes this rank's IB neighborhood once s_read_ib_restart_data has loaded only the IBs it owns: every rank sends + !! its own IBs to, and receives the owned IBs of, each distinct rank in ib_neighbor_ranks, appending them to patch_ib after its + !! own. + subroutine s_build_ib_neighborhood() - ! Check for moving IBs across both namelist and particle cloud patches. - moving_immersed_boundary_flag = .false. - do i = 1, num_namelist_ibs - if (patch_ib(i)%moving_ibm /= 0) then - moving_immersed_boundary_flag = .true. - exit - end if - end do - if (.not. moving_immersed_boundary_flag) then - do i = 1, num_particle_clouds - if (particle_cloud(i)%moving_ibm /= 0) then - moving_immersed_boundary_flag = .true. - exit - end if - end do - end if +#ifdef MFC_MPI + integer, allocatable, dimension(:) :: nbr_ranks, recv_counts, requests + type(ib_patch_parameters), allocatable, dimension(:,:) :: recv_ibs + integer :: i, n_nbrs, nreqs, patch_bytes, ierr +#endif call s_compute_ib_neighbor_ranks() #ifdef MFC_MPI - if (num_procs == 1) then - ! single-rank: all patches are local; append particle bed entries directly into patch_ib. - do i = 1, num_particle_cloud_ibs - patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) + if (num_procs > 1) then + ! A rank can fill several table slots (periodicity, few ranks) or be its own neighbor; exchange once per distinct rank + allocate (nbr_ranks(size(ib_neighbor_ranks))) + nbr_ranks = reshape(ib_neighbor_ranks, [size(ib_neighbor_ranks)]) + n_nbrs = 0 + do i = 1, size(nbr_ranks) + if (nbr_ranks(i) < 0 .or. nbr_ranks(i) == proc_rank) cycle + if (any(nbr_ranks(1:n_nbrs) == nbr_ranks(i))) cycle + n_nbrs = n_nbrs + 1 + nbr_ranks(n_nbrs) = nbr_ranks(i) end do - num_gbl_ibs = num_namelist_ibs + num_particle_cloud_ibs - @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & - & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - num_ibs = num_gbl_ibs - num_local_ibs = num_gbl_ibs - do i = 1, num_gbl_ibs - local_ib_patch_ids(i) = i + + allocate (recv_counts(n_nbrs), requests(2*n_nbrs)) + do i = 1, n_nbrs + call MPI_IRECV(recv_counts(i), 1, MPI_INTEGER, nbr_ranks(i), 500, MPI_COMM_WORLD, requests(2*i - 1), ierr) + call MPI_ISEND(num_local_ibs, 1, MPI_INTEGER, nbr_ranks(i), 500, MPI_COMM_WORLD, requests(2*i), ierr) end do - else - ! multi-rank: compact namelist patches in-place (write_idx <= read_idx, no aliasing), then append local particle beds. - num_ibs = 0 - num_local_ibs = 0 - num_gbl_ibs = num_namelist_ibs + num_bed_ibs - do i = 1, num_namelist_ibs - centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid - if (f_neighborhood_ranks_own_location(centroid)) then - num_ibs = num_ibs + 1 - patch_ib(num_ibs) = patch_ib(i) - patch_ib(num_ibs)%gbl_patch_id = i - if (f_local_rank_owns_location(centroid)) then - num_local_ibs = num_local_ibs + 1 - @:PROHIBIT(num_local_ibs > num_local_ibs_max, & - & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") - local_ib_patch_ids(num_local_ibs) = num_ibs - end if + call MPI_WAITALL(2*n_nbrs, requests, MPI_STATUSES_IGNORE, ierr) + + patch_bytes = storage_size(patch_ib(1))/8 + allocate (recv_ibs(max(1, maxval(recv_counts)), n_nbrs)) + nreqs = 0 + do i = 1, n_nbrs + if (recv_counts(i) > 0) then + nreqs = nreqs + 1 + call MPI_IRECV(recv_ibs(:,i), recv_counts(i)*patch_bytes, MPI_BYTE, nbr_ranks(i), 501, MPI_COMM_WORLD, & + & requests(nreqs), ierr) end if - end do - ! particle_cloud_ibs entries already passed the neighborhood check at generation time so no need to recheck it here. - do i = 1, num_particle_cloud_ibs - centroid = [particle_cloud_ibs(i)%x_centroid, particle_cloud_ibs(i)%y_centroid, 0._wp] - if (num_dims == 3) centroid(3) = particle_cloud_ibs(i)%z_centroid - num_ibs = num_ibs + 1 - @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, & - & "Local IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - patch_ib(num_ibs) = particle_cloud_ibs(i) - if (f_local_rank_owns_location(centroid)) then - num_local_ibs = num_local_ibs + 1 - @:PROHIBIT(num_local_ibs > num_local_ibs_max, & - & "Too many IBs on a single processor rank. Modify case file or increase limit of num_local_ibs_max to resolve.") - local_ib_patch_ids(num_local_ibs) = num_ibs + if (num_local_ibs > 0) then + nreqs = nreqs + 1 + call MPI_ISEND(patch_ib, num_local_ibs*patch_bytes, MPI_BYTE, nbr_ranks(i), 501, MPI_COMM_WORLD, & + & requests(nreqs), ierr) end if end do + call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr) + + do i = 1, n_nbrs + @:PROHIBIT(num_ibs + recv_counts(i) > num_ib_patches_max_namelist, & + & "IB neighborhood exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") + patch_ib(num_ibs + 1:num_ibs + recv_counts(i)) = recv_ibs(1:recv_counts(i),i) + num_ibs = num_ibs + recv_counts(i) + end do + + deallocate (nbr_ranks, recv_counts, requests, recv_ibs) end if -#else - ! no-MPI: all patches are local; append particle bed entries directly into patch_ib. - do i = 1, num_particle_cloud_ibs - patch_ib(num_namelist_ibs + i) = particle_cloud_ibs(i) - end do - num_gbl_ibs = num_namelist_ibs + num_particle_cloud_ibs - @:PROHIBIT(num_gbl_ibs > num_ib_patches_max_namelist, & - & "Total IB count exceeds patch_ib capacity. Increase num_ib_patches_max_namelist.") - num_ibs = num_gbl_ibs - num_local_ibs = num_gbl_ibs - do i = 1, num_gbl_ibs - local_ib_patch_ids(i) = i - end do #endif @:ALLOCATE(ib_gbl_idx_lookup(1:num_gbl_ibs)) - end subroutine s_reduce_ib_patch_array + end subroutine s_build_ib_neighborhood !> Build ib_neighbor_ranks(-1:1,-1:1,-1:1): MPI ranks of all neighbor domains. Uses two rounds of MPI_SENDRECV cascades - face !! neighbors are known from bc_*, edge neighbors are obtained in round 1, and (3D) corner neighbors in round 2. diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 1da790a90..c048f8029 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -1872,6 +1872,8 @@ def check_parallel_io_pre_process(self): m = self.get("m", 0) n = self.get("n", 0) + self.prohibit(file_per_process and not parallel_io, "file_per_process requires parallel_io = T") + if down_sample: self.prohibit(not parallel_io, "down sample requires parallel_io = T") self.prohibit(not igr, "down sample requires igr = T") diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 67559d36e..a0f080952 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -1359,6 +1359,9 @@ def _decl(targets: set, *names: str) -> None: "avg_state", "alt_soundspeed", "mixture_err", +) +_nv( + _ALL, "num_particle_clouds", "particle_cloud", ) diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index fff3b4552..fe1dc39d6 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -1309,6 +1309,38 @@ def alter_ib(dimInfo, six_eqn_model=False, viscous=False): ) ) + # Restart roundtrip regression: particle-cloud beds must survive a restart, not just namelist patch_ib patches - + # pre_process now generates them once and simulation reads that layout back on every start (fresh or restart). + cases.append( + define_case_d( + stack, + "IBM -> Particle Cloud -> Box -> Restart", + { + "ib": "T", + "num_ibs": 0, + "num_particle_clouds": 1, + "fd_order": 2, + "n": 49, + "particle_cloud(1)%cloud_geometry": 1, + "particle_cloud(1)%packing_method": 1, + "particle_cloud(1)%x_centroid": 0.5, + "particle_cloud(1)%y_centroid": 0.5, + "particle_cloud(1)%length_x": 0.6, + "particle_cloud(1)%length_y": 0.6, + "particle_cloud(1)%num_particles": 4, + "particle_cloud(1)%radius": 0.02, + "particle_cloud(1)%mass": 1.0, + "particle_cloud(1)%min_spacing": 0.005, + "particle_cloud(1)%moving_ibm": 0, + "particle_cloud(1)%seed": 12345, + "patch_icpp(1)%vel(1)": 0.001, + "patch_icpp(2)%vel(1)": 0.001, + "patch_icpp(3)%vel(1)": 0.001, + }, + restart_check=True, + ) + ) + if len(dimInfo[0]) == 3 and not viscous: cases.append( define_case_d(