diff --git a/benchmarks/ibm/core b/benchmarks/ibm/core new file mode 100644 index 0000000000..e69de29bb2 diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 219bddb1f1..1cbb991a84 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -108,6 +108,7 @@ is equivalent to `"riemann_solver": 2`. Defined names appear in each parameter's | Parameter | Type | Description | | ---: | :----: | :--- | | `run_time_info` | Logical | Output run-time information | +| `phase_timing_wrt` | Logical | Append per-phase wall times to `phase_time_data.dat` | | `rdma_mpi` | Logical | (GPUs) Enable RDMA for MPI communication. | | `case_dir` | String | Case directory path | | `old_grid` | Logical | Use grid from previous simulation | @@ -116,6 +117,7 @@ is equivalent to `"riemann_solver": 2`. Defined names appear in each parameter's | `n_start_old` | Integer | Starting index from previous simulation | - `run_time_info` generates a text file that includes run-time information including the CFL number(s) at each time-step. +- `phase_timing_wrt` times every NVTX range of the time march on each rank and, at the end of the run, appends one row per range to `phase_time_data.dat` with the min, mean and max over ranks of its inclusive wall time per step. Nothing is communicated until the end of the run. Plot the rows with `misc/plot_phase_timings.py`. - `rdma_mpi` optimizes data transfers between GPUs using Remote Direct Memory Access (RDMA). The underlying MPI implementation and communication infrastructure must support this feature, detecting GPU pointers and performing RDMA accordingly. diff --git a/misc/plot_phase_timings.py b/misc/plot_phase_timings.py new file mode 100644 index 0000000000..3d2f59e0b4 --- /dev/null +++ b/misc/plot_phase_timings.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 +"""Plot the max-over-ranks wall time per step of each phase against rank count. + +Reads one or more phase_time_data.dat files written by simulation with phase_timing_wrt = T +(a directory is read as /phase_time_data.dat). Rows from all files are pooled; when +several runs share a rank count, the last one read wins. +""" + +import argparse +import os +from collections import defaultdict + +import matplotlib.pyplot as plt + + +def read_rows(paths): + data = defaultdict(dict) # phase -> {ranks: (mean, max)} + for path in paths: + if os.path.isdir(path): + path = os.path.join(path, "phase_time_data.dat") + with open(path) as f: + next(f) + for line in f: + ranks, _, _, _, t_mean, t_max, phase = line.split() + data[phase][int(ranks)] = (float(t_mean), float(t_max)) + return data + + +def main(): + parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) + parser.add_argument("paths", nargs="+", help="phase_time_data.dat files or the case directories holding them") + parser.add_argument("--phases", nargs="+", help="phases to plot (default: --top largest)") + parser.add_argument("--top", type=int, default=10, help="plot the N phases with the largest max time at the largest rank count") + parser.add_argument("--mean", action="store_true", help="also plot the mean over ranks (dashed)") + parser.add_argument("--relative", action="store_true", help="divide each phase by its time at the smallest rank count") + parser.add_argument("-o", "--output", default="phase_timings.png", help="output image") + args = parser.parse_args() + + data = read_rows(args.paths) + phases = args.phases or sorted(data, key=lambda ph: data[ph][max(data[ph])][1], reverse=True)[: args.top] + + fig, ax = plt.subplots(figsize=(9, 6)) + for phase in phases: + ranks = sorted(data[phase]) + t_mean, t_max = zip(*(data[phase][r] for r in ranks)) + scale = t_max[0] if args.relative else 1.0 + (line,) = ax.plot(ranks, [t / scale for t in t_max], "o-", label=phase) + if args.mean: + ax.plot(ranks, [t / scale for t in t_mean], "--", color=line.get_color()) + + ax.set_xscale("log", base=2) + ax.set_xlabel("Ranks") + ax.set_ylabel("max time / time at fewest ranks" if args.relative else "max over ranks [s/step]") + ax.grid(True, which="both", alpha=0.3) + ax.legend(fontsize="small") + fig.tight_layout() + fig.savefig(args.output, dpi=150) + print(f"Wrote {args.output}") + + +if __name__ == "__main__": + main() diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 09d598b19b..b3454705e4 100644 --- a/src/common/m_constants.fpp +++ b/src/common/m_constants.fpp @@ -28,8 +28,8 @@ module m_constants integer, parameter :: num_stl_models_max = 10 !> Maximum number of immersed boundary patches (legacy, not used for patch_ib sizing) !> Fixed capacity of patch_ib (namelist patches + local particle bed subset after reduction) - integer, parameter :: num_local_ibs_max = 8000 !< Maximum number of immersed boundary patches (patch_ib) - integer, parameter :: num_ib_patches_max_namelist = 216000 + integer, parameter :: num_local_ibs_max = 30000 !< Maximum number of immersed boundary patches (patch_ib) + integer, parameter :: num_ib_patches_max_namelist = 810000 integer, parameter :: num_particle_clouds_max = 10 !< Maximum number of particle bed patch specifications integer, parameter :: num_bc_patches_max = 10 !< Maximum number of boundary condition patches integer, parameter :: max_2d_fourier_modes = 10 !< Max Fourier mode index for 2D modal patch (geometry 13) diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 295ca1d9de..97ea3a1446 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -340,8 +340,9 @@ module m_derived_types end type ib_stl_parameters type ib_patch_parameters - integer :: geometry !< Type of geometry for the patch + integer :: geometry !< Type of geometry for the patch integer :: gbl_patch_id + integer :: owner_rank !< MPI rank whose subdomain holds the centroid; the only rank that sums this IB's force real(wp) :: x_centroid, y_centroid, z_centroid !< Geometric center coordinates of the patch !> Centroid locations of intermediate steps in the time_stepper module diff --git a/src/common/m_helper.fpp b/src/common/m_helper.fpp index 17a986d6a9..eec1fe7df7 100644 --- a/src/common/m_helper.fpp +++ b/src/common/m_helper.fpp @@ -20,7 +20,7 @@ module m_helper & 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_local_rank_owns_location + & f_local_rank_owns_location, s_sort_int_key_value contains @@ -366,6 +366,29 @@ contains end subroutine s_swap + !> Sort the key-value pair by the key + pure subroutine s_sort_int_key_value(keys, vals, n) + + integer, dimension(:), intent(inout) :: keys, vals + integer, intent(in) :: n + integer :: i, j, key, val + + do i = 2, n + key = keys(i); val = vals(i) + + j = i + do while (j > 1) + if (keys(j - 1) <= key) exit + j = j - 1 + end do + + keys(j + 1:i) = keys(j:i - 1) + vals(j + 1:i) = vals(j:i - 1) + keys(j) = key; vals(j) = val + end do + + end subroutine s_sort_int_key_value + !> Create a transformation matrix. function f_create_transform_matrix(param, center) result(out_matrix) diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp index 15bec952af..8d8445edc0 100644 --- a/src/common/m_mpi_common.fpp +++ b/src/common/m_mpi_common.fpp @@ -283,6 +283,63 @@ contains end subroutine mpi_bcast_time_step_values + !> Reduce the per-rank NVTX range timers onto rank 0 over the union of range names seen by any rank + impure subroutine s_mpi_reduce_nvtx_timers(names, num_names, t_min, t_sum, t_max, calls_max) + + character(len=nvtx_name_len), intent(out) :: names(nvtx_max_timers) + integer, intent(out) :: num_names + real(kind=8), dimension(nvtx_max_timers), intent(out) :: t_min, t_sum, t_max + integer(kind=8), dimension(nvtx_max_timers), intent(out) :: calls_max + real(kind=8) :: sec(nvtx_num_timers), t_loc(nvtx_max_timers) + integer(kind=8) :: calls_loc(nvtx_max_timers) + integer :: i, j, owner + +#ifdef MFC_MPI + integer :: ierr !< Generic flag used to identify and report MPI errors +#endif + + ! Each round, the lowest rank holding an unlisted name broadcasts it + num_names = 0 + do + owner = num_procs + do i = 1, nvtx_num_timers + if (f_nvtx_find(names(1:num_names), nvtx_timer_names(i)) == 0) then + owner = proc_rank + exit + end if + end do +#ifdef MFC_MPI + call MPI_ALLREDUCE(MPI_IN_PLACE, owner, 1, MPI_INTEGER, MPI_MIN, MPI_COMM_WORLD, ierr) +#endif + if (owner == num_procs .or. num_names == nvtx_max_timers) exit + num_names = num_names + 1 + if (proc_rank == owner) names(num_names) = nvtx_timer_names(i) +#ifdef MFC_MPI + call MPI_BCAST(names(num_names), nvtx_name_len, MPI_CHARACTER, owner, MPI_COMM_WORLD, ierr) +#endif + end do + + sec = f_nvtx_timer_seconds() + t_loc = 0._8 + calls_loc = 0_8 + do i = 1, num_names + j = f_nvtx_find(nvtx_timer_names(1:nvtx_num_timers), names(i)) + if (j == 0) cycle + t_loc(i) = sec(j) + calls_loc(i) = nvtx_timer_calls(j) + end do + +#ifdef MFC_MPI + call MPI_REDUCE(t_loc, t_min, num_names, MPI_DOUBLE_PRECISION, MPI_MIN, 0, MPI_COMM_WORLD, ierr) + call MPI_REDUCE(t_loc, t_sum, num_names, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, ierr) + call MPI_REDUCE(t_loc, t_max, num_names, MPI_DOUBLE_PRECISION, MPI_MAX, 0, MPI_COMM_WORLD, ierr) + call MPI_REDUCE(calls_loc, calls_max, num_names, MPI_INTEGER8, MPI_MAX, 0, MPI_COMM_WORLD, ierr) +#else + t_min = t_loc; t_sum = t_loc; t_max = t_loc; calls_max = calls_loc +#endif + + end subroutine s_mpi_reduce_nvtx_timers + !> Print a case file error with the prohibited condition and message, then abort execution. impure subroutine s_prohibit_abort(condition, message) diff --git a/src/common/m_nvtx.f90 b/src/common/m_nvtx.f90 index 66893c10d0..eec9c111a9 100644 --- a/src/common/m_nvtx.f90 +++ b/src/common/m_nvtx.f90 @@ -14,6 +14,15 @@ module m_nvtx character(len=256), private :: tempName + !> Per-rank wall-clock accumulation of every named range, enabled by phase_timing_wrt + integer, parameter :: nvtx_name_len = 64, nvtx_max_timers = 256 + logical :: nvtx_timing = .false. + integer :: nvtx_num_timers = 0 + character(len=nvtx_name_len) :: nvtx_timer_names(nvtx_max_timers) + integer(c_int64_t) :: nvtx_timer_ticks(nvtx_max_timers) = 0, nvtx_timer_calls(nvtx_max_timers) = 0 + integer, private :: depth = 0, stack_id(64) + integer(c_int64_t), private :: stack_tick(64) + type, bind(C) :: nvtxEventAttributes integer(c_int16_t) :: version = 1 integer(c_int16_t) :: size = 48 ! @@ -64,6 +73,8 @@ subroutine nvtxStartRange(name, id) integer, intent(in), optional :: id type(nvtxEventAttributes) :: event + if (nvtx_timing) call s_push_timer(name) + #if defined(MFC_GPU) && defined(__PGI) tempName = trim(name) // c_null_char @@ -85,6 +96,65 @@ subroutine nvtxEndRange call nvtxRangePop #endif + if (nvtx_timing .and. depth > 0) call s_pop_timer() + end subroutine nvtxEndRange + !> Index of name in list, or 0 if absent + pure integer function f_nvtx_find(list, name) result(idx) + + character(len=*), intent(in) :: list(:), name + + do idx = 1, size(list) + if (list(idx) == name) return + end do + idx = 0 + + end function f_nvtx_find + + !> Open a timed range, registering its name on first use (id 0 = untracked) + subroutine s_push_timer(name) + + character(len=*), intent(in) :: name + integer :: id + + id = f_nvtx_find(nvtx_timer_names(1:nvtx_num_timers), name) + if (id == 0 .and. nvtx_num_timers < nvtx_max_timers) then + nvtx_num_timers = nvtx_num_timers + 1 + nvtx_timer_names(nvtx_num_timers) = name + id = nvtx_num_timers + end if + depth = depth + 1 + stack_id(depth) = id + call system_clock(stack_tick(depth)) + + end subroutine s_push_timer + + !> Close the innermost timed range and accumulate its elapsed ticks + subroutine s_pop_timer() + + integer(c_int64_t) :: tick + integer :: id + + call system_clock(tick) + id = stack_id(depth) + if (id > 0) then + nvtx_timer_ticks(id) = nvtx_timer_ticks(id) + tick - stack_tick(depth) + nvtx_timer_calls(id) = nvtx_timer_calls(id) + 1 + end if + depth = depth - 1 + + end subroutine s_pop_timer + + !> Accumulated seconds of every registered range + function f_nvtx_timer_seconds() result(sec) + + real(c_double) :: sec(nvtx_num_timers) + integer(c_int64_t) :: rate + + call system_clock(count_rate=rate) + sec = real(nvtx_timer_ticks(1:nvtx_num_timers), c_double)/real(rate, c_double) + + end function f_nvtx_timer_seconds + end module m_nvtx diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp index 1a7a24f589..f5e41fee51 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_neighborhood_ranks_own_location, ib_gbl_idx_lookup, collisions_active + & f_neighborhood_ranks_own_location, collisions_active ! overlap distances for computing collisions integer, allocatable, dimension(:,:) :: collision_lookup real(wp), allocatable, dimension(:,:) :: wall_overlap_distances @@ -29,9 +29,6 @@ module m_collisions $:GPU_DECLARE(create='[spring_stiffness, damping_parameter]') $:GPU_DECLARE(create='[collision_lookup, wall_overlap_distances]') - integer, dimension(:), allocatable :: ib_gbl_idx_lookup - $:GPU_DECLARE(create='[ib_gbl_idx_lookup]') - !> true when any IB-IB or IB-wall contact was detected on this rank since the last adaptive-dt computation logical :: collisions_active @@ -251,7 +248,7 @@ contains integer, intent(out) :: num_considered_collisions integer :: i, j, k, z_bound, ii, jj, kk integer, dimension(2) :: decoded_pairs - integer :: gp_idx, gp_patch_id, neighbor_patch_id + integer :: gp_idx, gp_patch_id, neighbor_patch_id, local_idx integer :: pair_idx, out_idx logical :: already_found @@ -306,8 +303,10 @@ contains ! get the decoded pairs for checking if they exist, using ii,jj,kk as dummy indices call s_decode_patch_periodicity(raw_pairs(pair_idx, 1), decoded_pairs(1), ii, jj, kk) call s_decode_patch_periodicity(raw_pairs(pair_idx, 2), decoded_pairs(2), ii, jj, kk) - decoded_pairs(1) = ib_gbl_idx_lookup(decoded_pairs(1)) - decoded_pairs(2) = ib_gbl_idx_lookup(decoded_pairs(2)) + call s_get_neighborhood_idx(decoded_pairs(1), local_idx) + decoded_pairs(1) = local_idx + call s_get_neighborhood_idx(decoded_pairs(2), local_idx) + decoded_pairs(2) = local_idx ! skip self-collisions (an IB cannot collide with its own periodic image) if (decoded_pairs(1) == decoded_pairs(2)) cycle diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index e72afb668c..12589c77ec 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -353,6 +353,7 @@ contains ! Logistics (sim-specific) run_time_info = .false. + phase_timing_wrt = .false. t_step_old = dflt_int ! Computational domain parameters (sim-specific) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 86c80c5e6a..35628a63ff 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -19,11 +19,20 @@ module m_ib_patches use m_helper_basic use m_helper use m_mpi_common + use m_constants implicit none private; public :: s_apply_ib_patches, s_update_ib_rotation_matrix, s_instantiate_STL_models, s_decode_patch_periodicity, & - & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities, s_get_ib_bound + & s_encode_patch_periodicity, s_initialize_ib_airfoils, s_get_periodicities, s_get_ib_bound, s_get_neighborhood_idx, & + & s_update_ib_lookup, s_compact_ib_lookup, s_merge_ib_lookup + + !> lookup arrays for converting global IB indices to local indices + integer, dimension(num_ib_patches_max_namelist) :: ib_lookup_keys, ib_lookup_vals + $:GPU_DECLARE(create='[ib_lookup_keys, ib_lookup_vals]') + + !> Holds each step's arrivals while they are sorted and merged in. Host only. + integer, dimension(num_ib_patches_max_namelist) :: ib_new_keys, ib_new_vals contains @@ -728,6 +737,117 @@ contains end subroutine s_decode_patch_periodicity + !> binary search to retrieve the local IB patch index using the global index + subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx, num_entries) + + $:GPU_ROUTINE(parallelism='[seq]') + + integer, intent(in) :: gbl_idx + integer, intent(out) :: neighborhood_idx + integer, intent(in), optional :: num_entries + integer :: lo, hi, mid + + neighborhood_idx = -1 + lo = 1 + hi = num_ibs + if (present(num_entries)) hi = num_entries + + do while (lo <= hi) + mid = lo + (hi - lo)/2 + if (ib_lookup_keys(mid) == gbl_idx) then + neighborhood_idx = ib_lookup_vals(mid) + return + else if (ib_lookup_keys(mid) < gbl_idx) then + lo = mid + 1 + else + hi = mid - 1 + end if + end do + + end subroutine s_get_neighborhood_idx + + !> Completely rebuilds the ib lookup map, used at startup + subroutine s_update_ib_lookup() + + integer :: i + + @:PROHIBIT(num_ibs > num_ib_patches_max_namelist, & + & "num_ibs exceeds the IB lookup capacity. Increase num_ib_patches_max_namelist.") + + do i = 1, num_ibs + ib_lookup_keys(i) = patch_ib(i)%gbl_patch_id + ib_lookup_vals(i) = i + end do + call s_sort_int_key_value(ib_lookup_keys, ib_lookup_vals, num_ibs) + + $:GPU_UPDATE(device='[ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs)]') + + end subroutine s_update_ib_lookup + + !> Drop the entries whose patches left the neighborhood and renumber the survivors onto the patch_ib slots they were compacted + !! into. Keys are never reordered, so the map stays sorted for free and only the values move: O(num_ibs_old) against re-sorting + !! the whole map. + subroutine s_compact_ib_lookup(old_to_new, num_ibs_old) + + integer, dimension(:), intent(in) :: old_to_new !< old patch_ib slot -> new slot, -1 if dropped + integer, intent(in) :: num_ibs_old + integer :: i, k + + k = 0 + do i = 1, num_ibs_old + if (old_to_new(ib_lookup_vals(i)) < 0) cycle + k = k + 1 + ib_lookup_keys(k) = ib_lookup_keys(i) + ib_lookup_vals(k) = old_to_new(ib_lookup_vals(i)) + end do + @:ASSERT(k == num_ibs, 'IB lookup and patch_ib disagree on the surviving patch count') + + $:GPU_UPDATE(device='[ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs)]') + + end subroutine s_compact_ib_lookup + + !> Fold patch_ib(num_ibs_pre+1:num_ibs) into the map: sort just the arrivals, then merge the two sorted runs downward from the + !! top. O(num_ibs) plus the sort of the few arrivals. They are copied out first because the runs share this array and merging in + !! place would overwrite entries still to be read. + subroutine s_merge_ib_lookup(num_ibs_pre) + + integer, intent(in) :: num_ibs_pre + integer :: i, j, k, r + logical :: take_old + + call nvtxStartRange("MERGE-IB-LOOKUP") + + r = num_ibs - num_ibs_pre + if (r <= 0) return + + do i = 1, r + ib_new_keys(i) = patch_ib(num_ibs_pre + i)%gbl_patch_id + ib_new_vals(i) = num_ibs_pre + i + end do + call s_sort_int_key_value(ib_new_keys, ib_new_vals, r) + + i = num_ibs_pre; j = r; k = num_ibs + do while (j >= 1) + ! Fortran does not short-circuit .and., so the exhausted-head test stands on its own + take_old = .false. + if (i >= 1) take_old = ib_lookup_keys(i) > ib_new_keys(j) + + if (take_old) then + ib_lookup_keys(k) = ib_lookup_keys(i); ib_lookup_vals(k) = ib_lookup_vals(i) + i = i - 1 + else + ib_lookup_keys(k) = ib_new_keys(j); ib_lookup_vals(k) = ib_new_vals(j) + j = j - 1 + end if + k = k - 1 + end do + + $:GPU_UPDATE(device='[ib_lookup_keys(1:num_ibs), ib_lookup_vals(1:num_ibs)]') + + call nvtxEndRange() + + end subroutine s_merge_ib_lookup + !> Determine the periodic wrapping bounds in each direction subroutine s_get_periodicities(xp_lower, xp_upper, yp_lower, yp_upper, zp_lower, zp_upper) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 3b156d6085..8bfdc203ba 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -49,10 +49,12 @@ module m_ibm #endif logical :: moving_immersed_boundary_flag - ! IB MPI buffers - integer, allocatable :: send_ids(:), recv_ids(:) - real(wp), allocatable :: send_ft(:,:), recv_ft(:,:) - real(wp), allocatable :: recv_forces_snap(:,:), recv_torques_snap(:,:) + ! IB force reduction (s_communicate_ib_forces): distinct neighborhood ranks in ascending order, and buffers allocated once + integer, parameter :: ib_rec_len = 7 !< gbl_patch_id, force(1:3), torque(1:3) + integer :: n_ib_nbrs = 0 + integer, allocatable :: ib_nbrs(:), ib_send_counts(:), ib_reqs(:), ib_stats(:,:) + real(wp), allocatable :: ib_send_bufs(:,:), ib_recv_bufs(:,:) + private :: ib_rec_len, n_ib_nbrs, ib_nbrs, ib_send_counts, ib_reqs, ib_stats, ib_send_bufs, ib_recv_bufs contains @@ -110,11 +112,7 @@ contains ! allocate some arrays for MPI communication, if required by this simulation #ifdef MFC_MPI - if (num_procs > 1) then - @:ALLOCATE(send_ids(size(patch_ib)), send_ft(6, size(patch_ib))) - @:ALLOCATE(recv_forces_snap(size(patch_ib), 3), recv_torques_snap(size(patch_ib), 3), recv_ids(size(patch_ib)), & - & recv_ft(6, size(patch_ib))) - end if + if (num_procs > 1) call s_setup_ib_force_comm() #endif call s_update_ib_lookup() @@ -1294,10 +1292,14 @@ contains $:END_GPU_PARALLEL_LOOP() end if + call nvtxStartRange("APPLY-COLLISION-FORCES") call s_apply_collision_forces(ghost_points, num_gps, ib_markers, forces, torques) + call nvtxEndRange ! reduce the forces across local neighborhood ranks + call nvtxStartRange("COMMUNICATE-IB-FORCES") call s_communicate_ib_forces(forces, torques) + call nvtxEndRange ! consider body forces after reducing to avoid double counting do i = 1, num_ibs @@ -1493,153 +1495,156 @@ contains end subroutine s_wrap_periodic_ibs - !> @brief Swaps ownership of IBs and passes ownership of IBs to neighbor processors - !> Reduces forces and torques across the local neighborhood without a global allreduce. Accumulation phase: 2 passes per - !! dimension receiving from the low-index (-X) neighbor. Pass 1: add received values; save what was received as recv_snap. Pass - !! 2: send current (post-pass-1) values; add received; subtract recv_snap to remove double-counting of the direct contribution - !! already added in pass 1. Back-propagation phase: 2 passes per dimension receiving from the high-index (+X) neighbor, each - !! overwriting local forces with the neighbor's accumulated total. + !> Builds the ascending list of distinct neighborhood ranks used by s_communicate_ib_forces and allocates its buffers once. The + !! ascending order is the fixed order in which an owner adds its neighbors' contributions. + impure subroutine s_setup_ib_force_comm() + +#ifdef MFC_MPI + integer :: i, j, r, max_recs + integer, allocatable, dimension(:) :: flat + + if (num_procs == 1 .or. allocated(ib_nbrs)) return + @:PROHIBIT(num_gbl_ibs >= 2**min(digits(0._wp), 30), "Too many IBs to encode their ids exactly in real(wp) force messages") + + ! a rank can fill several table slots (periodicity, few ranks) or be its own neighbor; keep each distinct rank once + flat = reshape(ib_neighbor_ranks, [size(ib_neighbor_ranks)]) + allocate (ib_nbrs(size(flat))) + n_ib_nbrs = 0 + do i = 1, size(flat) + r = flat(i) + if (r < 0 .or. r == proc_rank) cycle + if (any(ib_nbrs(1:n_ib_nbrs) == r)) cycle + ! insert in ascending order + j = n_ib_nbrs + do while (j > 0) + if (ib_nbrs(j) < r) exit + ib_nbrs(j + 1) = ib_nbrs(j) + j = j - 1 + end do + ib_nbrs(j + 1) = r + n_ib_nbrs = n_ib_nbrs + 1 + end do + + ! an owner never receives more records than it owns, and a broadcast never carries more than its sender owns + max_recs = min(size(patch_ib), num_local_ibs_max) + allocate (ib_send_counts(n_ib_nbrs), ib_reqs(2*n_ib_nbrs), ib_stats(MPI_STATUS_SIZE, 2*n_ib_nbrs)) + allocate (ib_send_bufs(ib_rec_len*max_recs, max(n_ib_nbrs, 1)), ib_recv_bufs(ib_rec_len*max_recs, max(n_ib_nbrs, 1))) +#endif + + end subroutine s_setup_ib_force_comm + + !> Position of rank in ib_nbrs + function f_ib_nbr_slot(rank) result(slot) + + integer, intent(in) :: rank + integer :: slot + + do slot = 1, n_ib_nbrs + if (ib_nbrs(slot) == rank) return + end do + call s_mpi_abort('IB owner is not a neighborhood rank') + + end function f_ib_nbr_slot + + !> Reduces forces and torques so every rank holding an IB ends with bit-identical values. Each IB has one owner, the rank whose + !! subdomain holds its centroid (patch_ib%owner_rank, stamped at handoff). Phase 1: each rank sends its nonzero partials to the + !! owner, which adds them to its own partial in ascending source-rank order once all have arrived. Phase 2: each owner sends its + !! totals to every neighborhood rank, which overwrite their partials. One nonblocking message per neighbor per phase. subroutine s_communicate_ib_forces(forces, torques) real(wp), dimension(num_ibs, 3), intent(inout) :: forces, torques #ifdef MFC_MPI - integer :: i, j, k, l, pack_pos, unpack_pos, buf_size, ierr - integer :: send_neighbor, recv_neighbor, recv_count, tag - character(len=1), allocatable :: ib_force_send_buf(:), ib_force_recv_buf(:) - - if (num_procs == 1) return - - buf_size = storage_size(0)/8 + (storage_size(0)/8 + 6*storage_size(0._wp)/8)*size(patch_ib) - allocate (ib_force_send_buf(buf_size), ib_force_recv_buf(buf_size)) - - ! Accumulation phase: propagate contributions toward the high-index corner. - #:for X, ID in [('x', 1), ('y', 2), ('z', 3)] - if (num_dims >= ${ID}$) then - send_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0) - recv_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0) - - recv_forces_snap = 0._wp - recv_torques_snap = 0._wp - $:GPU_UPDATE(device='[recv_forces_snap, recv_torques_snap]') - tag = 300 - - do k = 1, min(2*ib_neighborhood_radius, num_procs_${X}$ - 1) - ! send forces to +${X}$ neighbor; receive from -${X}$ neighbor. Add received values then - pack_pos = 0 - if (num_ibs > 0) then - $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') - do i = 1, num_ibs - send_ids(i) = patch_ib(i)%gbl_patch_id - do l = 1, 3 - send_ft(l, i) = forces(i, l) - send_ft(l + 3, i) = torques(i, l) - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - $:GPU_UPDATE(host='[send_ids, send_ft]') - call MPI_PACK(num_ibs, 1, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) - call MPI_PACK(send_ids, num_ibs, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) - call MPI_PACK(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) - call MPI_SENDRECV(ib_force_send_buf, pack_pos, MPI_PACKED, send_neighbor, tag, ib_force_recv_buf, buf_size, & - & MPI_PACKED, recv_neighbor, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - - if (recv_neighbor /= MPI_PROC_NULL) then - unpack_pos = 0 - call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) - call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, MPI_INTEGER, & - & MPI_COMM_WORLD, ierr) - call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, MPI_COMM_WORLD, ierr) - $:GPU_UPDATE(device='[recv_ids(1:recv_count), recv_ft(:, 1:recv_count)]') - if (num_ibs > 0) then - $:GPU_PARALLEL_LOOP(private='[i, j, l]', copy='[forces, torques]') - do i = 1, recv_count - call s_get_neighborhood_idx(recv_ids(i), j) - if (j > 0) then - ! add forces and subtract recv_snap prevent double-counting - do l = 1, 3 - forces(j, l) = forces(j, l) + recv_ft(l, i) - recv_forces_snap(j, l) - torques(j, l) = torques(j, l) + recv_ft(l + 3, i) - recv_torques_snap(j, l) - recv_forces_snap(j, l) = recv_ft(l, i) - recv_torques_snap(j, l) = recv_ft(l + 3, i) - end do - end if - end do - $:END_GPU_PARALLEL_LOOP() - end if - end if - tag = tag + 2 - end do - end if - #:endfor - - ! Send final sums back to neighbors in -X direction - #:for X, ID in [('x', 1), ('y', 2), ('z', 3)] - if (num_dims >= ${ID}$) then - send_neighbor = merge(bc_${X}$%beg, MPI_PROC_NULL, bc_${X}$%beg >= 0) - recv_neighbor = merge(bc_${X}$%end, MPI_PROC_NULL, bc_${X}$%end >= 0) - - do k = 1, min(2*ib_neighborhood_radius, num_procs_${X}$ - 1) - pack_pos = 0 - if (num_ibs > 0) then - $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') - do i = 1, num_ibs - send_ids(i) = patch_ib(i)%gbl_patch_id - do l = 1, 3 - send_ft(l, i) = forces(i, l) - send_ft(l + 3, i) = torques(i, l) - end do - end do - $:END_GPU_PARALLEL_LOOP() - end if - $:GPU_UPDATE(host='[send_ids, send_ft]') - call MPI_PACK(num_ibs, 1, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) - call MPI_PACK(send_ids, num_ibs, MPI_INTEGER, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) - call MPI_PACK(send_ft, 6*num_ibs, mpi_p, ib_force_send_buf, buf_size, pack_pos, MPI_COMM_WORLD, ierr) - call MPI_SENDRECV(ib_force_send_buf, pack_pos, MPI_PACKED, send_neighbor, tag, ib_force_recv_buf, buf_size, & - & MPI_PACKED, recv_neighbor, tag, MPI_COMM_WORLD, MPI_STATUS_IGNORE, ierr) - if (recv_neighbor /= MPI_PROC_NULL) then - unpack_pos = 0 - call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_count, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) - call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ids, recv_count, MPI_INTEGER, & - & MPI_COMM_WORLD, ierr) - call MPI_UNPACK(ib_force_recv_buf, buf_size, unpack_pos, recv_ft, 6*recv_count, mpi_p, MPI_COMM_WORLD, ierr) - $:GPU_UPDATE(device='[recv_ids(1:recv_count), recv_ft(:, 1:recv_count)]') - if (num_ibs > 0) then - $:GPU_PARALLEL_LOOP(private='[i, j, l]', copy='[forces, torques]') - do i = 1, recv_count - call s_get_neighborhood_idx(recv_ids(i), j) - if (j > 0) then - do l = 1, 3 - forces(j, l) = recv_ft(l, i) - torques(j, l) = recv_ft(l + 3, i) - end do - end if - end do - $:END_GPU_PARALLEL_LOOP() - end if - end if - tag = tag + 2 - end do - end if - #:endfor + integer :: i, j, s, r, slot, pos, nvals, ierr + + if (num_procs == 1 .or. n_ib_nbrs == 0) return + + ! Phase 1: send each partial to the owner of its IB + do s = 1, n_ib_nbrs + call MPI_IRECV(ib_recv_bufs(:,s), size(ib_recv_bufs, 1), mpi_p, ib_nbrs(s), 610, MPI_COMM_WORLD, ib_reqs(s), ierr) + end do + + ib_send_counts = 0 + do i = 1, num_ibs + if (patch_ib(i)%owner_rank == proc_rank) cycle + if (all(forces(i,:) == 0._wp) .and. all(torques(i,:) == 0._wp)) cycle + slot = f_ib_nbr_slot(patch_ib(i)%owner_rank) + pos = ib_rec_len*ib_send_counts(slot) + ib_send_bufs(pos + 1, slot) = real(patch_ib(i)%gbl_patch_id, wp) + ib_send_bufs(pos + 2:pos + 4,slot) = forces(i,1:3) + ib_send_bufs(pos + 5:pos + 7,slot) = torques(i,1:3) + ib_send_counts(slot) = ib_send_counts(slot) + 1 + end do + + ! every neighbor gets a message, possibly empty, so every posted receive completes + do s = 1, n_ib_nbrs + call MPI_ISEND(ib_send_bufs(:,s), ib_rec_len*ib_send_counts(s), mpi_p, ib_nbrs(s), 610, MPI_COMM_WORLD, & + & ib_reqs(n_ib_nbrs + s), ierr) + end do + call MPI_WAITALL(2*n_ib_nbrs, ib_reqs, ib_stats, ierr) + + ! owner sum: own partial first, then each neighbor in ascending rank order, whatever order the messages arrived in + do s = 1, n_ib_nbrs + call MPI_GET_COUNT(ib_stats(:,s), mpi_p, nvals, ierr) + do r = 0, nvals/ib_rec_len - 1 + pos = ib_rec_len*r + call s_get_neighborhood_idx(nint(ib_recv_bufs(pos + 1, s)), j) + @:ASSERT(j > 0, 'IB force contribution for an IB this rank does not hold') + @:ASSERT(patch_ib(j)%owner_rank == proc_rank, 'IB force contribution for an IB this rank does not own') + forces(j,1:3) = forces(j,1:3) + ib_recv_bufs(pos + 2:pos + 4,s) + torques(j,1:3) = torques(j,1:3) + ib_recv_bufs(pos + 5:pos + 7,s) + end do + end do + + ! Phase 2: each owner sends its totals to every neighbor + do s = 1, n_ib_nbrs + call MPI_IRECV(ib_recv_bufs(:,s), size(ib_recv_bufs, 1), mpi_p, ib_nbrs(s), 611, MPI_COMM_WORLD, ib_reqs(s), ierr) + end do + + do i = 1, num_local_ibs + j = local_ib_patch_ids(i) + pos = ib_rec_len*(i - 1) + ib_send_bufs(pos + 1, 1) = real(patch_ib(j)%gbl_patch_id, wp) + ib_send_bufs(pos + 2:pos + 4,1) = forces(j,1:3) + ib_send_bufs(pos + 5:pos + 7,1) = torques(j,1:3) + end do + + do s = 1, n_ib_nbrs + call MPI_ISEND(ib_send_bufs(:,1), ib_rec_len*num_local_ibs, mpi_p, ib_nbrs(s), 611, MPI_COMM_WORLD, & + & ib_reqs(n_ib_nbrs + s), ierr) + end do + call MPI_WAITALL(2*n_ib_nbrs, ib_reqs, ib_stats, ierr) + + ! holders take the owner's bits + do s = 1, n_ib_nbrs + call MPI_GET_COUNT(ib_stats(:,s), mpi_p, nvals, ierr) + do r = 0, nvals/ib_rec_len - 1 + pos = ib_rec_len*r + call s_get_neighborhood_idx(nint(ib_recv_bufs(pos + 1, s)), j) + if (j <= 0) cycle + @:ASSERT(patch_ib(j)%owner_rank == ib_nbrs(s), 'IB force total from a rank that does not own the IB') + forces(j,1:3) = ib_recv_bufs(pos + 2:pos + 4,s) + torques(j,1:3) = ib_recv_bufs(pos + 5:pos + 7,s) + end do + end do #endif end subroutine s_communicate_ib_forces + !> @brief Swaps ownership of IBs and passes ownership of IBs to neighbor processors subroutine s_handoff_ib_ownership() - integer :: i, j, k, output_idx, local_output_idx - integer :: old_num_local_ibs - integer :: new_count, recv_count - integer :: pack_pos, unpack_pos, buf_size, patch_bytes - integer :: send_neighbor, recv_neighbor, ierr - integer :: dx, dy, dz, tag, nbr_idx, nreqs - real(wp), dimension(3) :: centroid - logical :: is_new - type(ib_patch_parameters) :: tmp_patch - integer, dimension(num_local_ibs_max) :: local_ib_idx_old + integer :: i, j, k, output_idx, local_output_idx + integer :: old_num_local_ibs, num_ibs_old, num_ibs_pre + integer :: new_count, recv_count + integer :: pack_pos, unpack_pos, buf_size, patch_bytes + integer :: send_neighbor, recv_neighbor, ierr + integer :: dx, dy, dz, tag, nbr_idx, nreqs + real(wp), dimension(3) :: centroid + logical :: is_new + type(ib_patch_parameters) :: tmp_patch + integer, dimension(num_local_ibs_max) :: local_ib_idx_old + integer, dimension(num_ib_patches_max_namelist) :: old_to_new ! old patch_ib slot -> slot after compaction ! 26 neighbors max in 3D (8 in 2D); each gets its own recv buffer integer, parameter :: max_nbrs = 26 character(len=1), allocatable :: send_buf(:), recv_bufs(:,:) @@ -1647,6 +1652,8 @@ contains integer, dimension(max_nbrs) :: recv_neighbor_list #ifdef MFC_MPI + call nvtxStartRange("HANDOFF-IB-OWNERSHIP") + if (num_procs > 1) then ! save a copy of the local IB's global indices to cross-reference for later. local_ib_idx_old = 0 @@ -1660,15 +1667,18 @@ contains $:GPU_UPDATE(host='[patch_ib]') ! delete any particles that no longer need to be tracked and coalesce the array + num_ibs_old = num_ibs output_idx = 0 local_output_idx = 0 do i = 1, num_ibs + old_to_new(i) = -1 centroid = [patch_ib(i)%x_centroid, patch_ib(i)%y_centroid, 0._wp] if (num_dims == 3) centroid(3) = patch_ib(i)%z_centroid ! delete if not in neighborhood if (f_neighborhood_ranks_own_location(centroid)) then output_idx = output_idx + 1 + old_to_new(i) = output_idx if (i /= output_idx) then patch_ib(output_idx) = patch_ib(i) end if @@ -1679,14 +1689,14 @@ contains @:PROHIBIT(local_output_idx > 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(local_output_idx) = output_idx + patch_ib(output_idx)%owner_rank = proc_rank end if end if end do num_ibs = output_idx num_local_ibs = local_output_idx - ! num_ibs shrinks here, so refresh it with patch_ib: s_update_ib_lookup scatters over it on the device $:GPU_UPDATE(device='[patch_ib, num_ibs]') - call s_update_ib_lookup() + call s_compact_ib_lookup(old_to_new, num_ibs_old) ! Broadcast newly-owned patches to all neighborhood neighbors patch_bytes = storage_size(tmp_patch)/8 @@ -1759,42 +1769,36 @@ contains call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr) ! Unpack all received buffers - do nbr_idx = 1, ((2*ib_neighborhood_radius + 1)**num_dims) - 1 + num_ibs_pre = num_ibs + do nbr_idx = 1, merge(26, 8, num_dims == 3) if (recv_neighbor_list(nbr_idx) == MPI_PROC_NULL) cycle unpack_pos = 0 call MPI_UNPACK(recv_bufs(:,nbr_idx), buf_size, unpack_pos, recv_count, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) do i = 1, recv_count call MPI_UNPACK(recv_bufs(:,nbr_idx), buf_size, unpack_pos, tmp_patch, patch_bytes, MPI_BYTE, MPI_COMM_WORLD, & & ierr) - call s_get_neighborhood_idx(tmp_patch%gbl_patch_id, j) + call s_get_neighborhood_idx(tmp_patch%gbl_patch_id, j, num_ibs_pre) if (j < 0) then num_ibs = num_ibs + 1 @:ASSERT(num_ibs <= size(patch_ib), 'patch_ib overflow in neighborhood handoff') patch_ib(num_ibs) = tmp_patch + else + ! an IB we already hold changed owner: record the new owner, keep our own copy of its state + patch_ib(j)%owner_rank = tmp_patch%owner_rank end if end do end do deallocate (send_buf, recv_bufs) $:GPU_UPDATE(device='[patch_ib, num_ibs]') - call s_update_ib_lookup() + call s_merge_ib_lookup(num_ibs_pre) end if + + call nvtxEndRange() #endif end subroutine s_handoff_ib_ownership - subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx) - - $:GPU_ROUTINE(parallelism='[seq]') - - integer, intent(in) :: gbl_idx - integer, intent(out) :: neighborhood_idx - integer :: i - - neighborhood_idx = ib_gbl_idx_lookup(gbl_idx) - - end subroutine s_get_neighborhood_idx - !> Abort if any immersed boundary marked no cell anywhere in the domain. !! !! A rank is given a patch when the patch CENTROID falls in its share of the domain, but for an STL @@ -1836,31 +1840,12 @@ contains end subroutine s_check_every_patch_marked - subroutine s_update_ib_lookup() - - integer :: i - - ib_gbl_idx_lookup = -1 - $:GPU_UPDATE(device='[ib_gbl_idx_lookup]') - - $:GPU_PARALLEL_LOOP(private='[i]') - do i = 1, num_ibs - ib_gbl_idx_lookup(patch_ib(i)%gbl_patch_id) = i - end do - $:END_GPU_PARALLEL_LOOP() - - $:GPU_UPDATE(host='[ib_gbl_idx_lookup]') - - end subroutine s_update_ib_lookup - !> Finalize the IBM module impure subroutine s_finalize_ibm_module() integer :: i @:DEALLOCATE(ib_markers%sf) - @:DEALLOCATE(corrected_gps%sf) - @:DEALLOCATE(ib_gbl_idx_lookup) do i = 1, num_ib_airfoils_max if (allocated(ib_airfoil_grids(i)%upper)) then @:DEALLOCATE(ib_airfoil_grids(i)%upper) @@ -1875,10 +1860,7 @@ contains end if if (collision_model > 0) call s_finalize_collisions_module() #ifdef MFC_MPI - if (num_procs > 1) then - @:DEALLOCATE(send_ids, send_ft) - @:DEALLOCATE(recv_forces_snap, recv_torques_snap, recv_ids, recv_ft) - end if + if (allocated(ib_nbrs)) deallocate (ib_nbrs, ib_send_counts, ib_reqs, ib_stats, ib_send_bufs, ib_recv_bufs) #endif end subroutine s_finalize_ibm_module diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index df311bca2b..d24ff9d668 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -58,7 +58,8 @@ module m_start_up private; public :: s_read_input_file, s_check_input_file, s_read_data_files, s_read_serial_data_files, & & s_read_parallel_data_files, s_initialize_internal_energy_equations, s_initialize_modules, s_initialize_gpu_vars, & - & s_initialize_mpi_domain, s_finalize_modules, s_perform_time_step, s_save_data, s_save_performance_metrics + & s_initialize_mpi_domain, s_finalize_modules, s_perform_time_step, s_save_data, s_save_performance_metrics, & + & s_save_phase_timings type(scalar_field), allocatable, dimension(:) :: q_cons_temp real(wp) :: dt_init @@ -722,6 +723,36 @@ contains end subroutine s_save_performance_metrics + !> Append the min/mean/max over ranks of each NVTX range's wall time per step to phase_time_data.dat + impure subroutine s_save_phase_timings(num_steps) + + integer, intent(in) :: num_steps + character(len=nvtx_name_len) :: names(nvtx_max_timers) + real(kind=8), dimension(nvtx_max_timers) :: t_min, t_sum, t_max + integer(kind=8), dimension(nvtx_max_timers) :: calls_max + integer :: i, num_names + logical :: file_exists + + nvtx_timing = .false. + call s_mpi_reduce_nvtx_timers(names, num_names, t_min, t_sum, t_max, calls_max) + if (proc_rank /= 0) return + + inquire (FILE='phase_time_data.dat', EXIST=file_exists) + if (file_exists) then + open (1, file='phase_time_data.dat', position='append', status='old') + else + open (1, file='phase_time_data.dat', status='new') + write (1, '(A10, A10, A14, 3(A16), 2X, A)') "Ranks", "Steps", "Calls/step", "min_s/step", "mean_s/step", & + & "max_s/step", "Phase" + end if + do i = 1, num_names + write (1, '(I10, I10, F14.3, 3(ES16.6), 2X, A)') num_procs, num_steps, real(calls_max(i), 8)/max(num_steps, 1), & + & [t_min(i), t_sum(i)/num_procs, t_max(i)]/max(num_steps, 1), trim(names(i)) + end do + close (1) + + end subroutine s_save_phase_timings + !> Save conservative variable data to disk at the current time step impure subroutine s_save_data(t_step, start, finish, io_time_avg, nt) @@ -1307,6 +1338,9 @@ contains nbr_ranks(n_nbrs) = nbr_ranks(i) end do + ! the IBs loaded from restart are this rank's own; neighbors learn their owner from the sender below + patch_ib(1:num_local_ibs)%owner_rank = proc_rank + 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) @@ -1335,6 +1369,7 @@ contains @: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) + patch_ib(num_ibs + 1:num_ibs + recv_counts(i))%owner_rank = nbr_ranks(i) num_ibs = num_ibs + recv_counts(i) end do @@ -1342,8 +1377,6 @@ contains end if #endif - @:ALLOCATE(ib_gbl_idx_lookup(1:num_gbl_ibs)) - 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 diff --git a/src/simulation/p_main.fpp b/src/simulation/p_main.fpp index 0111543fee..2e9937b67d 100644 --- a/src/simulation/p_main.fpp +++ b/src/simulation/p_main.fpp @@ -59,6 +59,7 @@ program p_main call nvtxEndRange ! INIT + nvtx_timing = phase_timing_wrt call nvtxStartRange("SIMULATION-TIME-MARCH") ! Time-stepping Loop do @@ -93,6 +94,8 @@ program p_main call nvtxEndRange ! Simulation + if (phase_timing_wrt) call s_save_phase_timings(t_step - t_step_start) + deallocate (proc_time, io_proc_time) call nvtxStartRange("FINALIZE-MODULES") diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 83c6bfa64d..e6cbdf4758 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -701,7 +701,7 @@ def _load(): _r("precision", INT, {"output"}) _r("format", INT, {"output"}) _r("ib_force_stride", INT, {"output", "ib"}) - for n in ["parallel_io", "file_per_process", "run_time_info", "prim_vars_wrt", "cons_vars_wrt", "fft_wrt", "ib_state_wrt", "ib_force_wrt"]: + for n in ["parallel_io", "file_per_process", "run_time_info", "phase_timing_wrt", "prim_vars_wrt", "cons_vars_wrt", "fft_wrt", "ib_state_wrt", "ib_force_wrt"]: _r(n, LOG, {"output"}) for n in [ "schlieren_wrt", @@ -1403,6 +1403,7 @@ def _decl(targets: set, *names: str) -> None: "hyper_cleaning_speed", "hyper_cleaning_tau", "run_time_info", + "phase_timing_wrt", "bubble_model", "lag_params", "probe_wrt", diff --git a/toolchain/mfc/params/descriptions.py b/toolchain/mfc/params/descriptions.py index e5eaf3a1c1..7a2c593c78 100644 --- a/toolchain/mfc/params/descriptions.py +++ b/toolchain/mfc/params/descriptions.py @@ -115,6 +115,7 @@ "relativity": "Enable special relativity", # Output "run_time_info": "Output run-time information", + "phase_timing_wrt": "Append per-phase wall time (min/mean/max over ranks) to phase_time_data.dat", "prim_vars_wrt": "Write primitive variables", "cons_vars_wrt": "Write conservative variables", "probe_wrt": "Write probe data",