From 318e07c5c8956fede1c8a5f242289651a39ac5e8 Mon Sep 17 00:00:00 2001 From: "Daniel J. Vickers" Date: Tue, 25 Aug 2026 11:49:28 -0400 Subject: [PATCH 01/10] Better data restarting --- src/simulation/m_data_output.fpp | 154 ++++++++++++++++++++++++---- src/simulation/m_particle_cloud.fpp | 2 +- src/simulation/m_start_up.fpp | 61 +++++++++-- 3 files changed, 190 insertions(+), 27 deletions(-) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 22bde22e1..632d248ec 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -196,28 +196,32 @@ contains do l = 0, p do k = 0, n do j = 0, m - call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) + ! Cells inside/on an immersed boundary hold ghost-derived, non-physical state - + ! excluded here so they cannot spuriously trip a stability violation. + if ((.not. ib) .or. (ib_markers%sf(j, k, l) == 0)) then + call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) + + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv) + + if (any_non_newtonian) then + Re(1) = 0._wp + do fl = 1, num_fluids + if (is_non_newtonian(fl)) then + Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) + else + Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) + end if + end do + Re(1) = 1._wp/max(Re(1), sgm_eps) + end if - call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv) + call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl, vcfl, Rc, ccfl) - if (any_non_newtonian) then - Re(1) = 0._wp - do fl = 1, num_fluids - if (is_non_newtonian(fl)) then - Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) - else - Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) - end if - end do - Re(1) = 1._wp/max(Re(1), sgm_eps) + icfl_max_loc = max(icfl_max_loc, icfl) + vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous)) + ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension)) + Rc_min_loc = min(Rc_min_loc, merge(Rc, huge(1.0_wp), viscous)) end if - - call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl, vcfl, Rc, ccfl) - - icfl_max_loc = max(icfl_max_loc, icfl) - vcfl_max_loc = max(vcfl_max_loc, merge(vcfl, 0.0_wp, viscous)) - ccfl_max_loc = max(ccfl_max_loc, merge(ccfl, 0.0_wp, surface_tension)) - Rc_min_loc = min(Rc_min_loc, merge(Rc, huge(1.0_wp), viscous)) end do end do end do @@ -246,6 +250,13 @@ contains if (Rc_min_glb < Rc_min) Rc_min = Rc_min_glb end if + ! Any rank whose own local extremum violates the limit is, by construction of the + ! max-reduction above, a rank that actually contains the offending cell(s). + if ((.not. f_approx_equal(icfl_max_loc, icfl_max_loc)) .or. icfl_max_loc > 1._wp) then + call s_report_icfl_violation(q_prim_vf) + end if + call s_mpi_barrier() ! ensure diagnostic output above is flushed before any rank aborts below + if (proc_rank == 0) then write (3, '(13X,I9,13X,F10.6,13X,F10.6,13X,F10.6)', advance="no") t_step, dt, mytime, icfl_max_glb @@ -290,6 +301,111 @@ contains end subroutine s_write_run_time_information + !> Locate the grid cell responsible for an ICFL violation on this rank and report its state plus the nearest immersed-boundary + !! particles, to aid debugging stability failures in particle-laden high-Mach cases. + impure subroutine s_report_icfl_violation(q_prim_vf) + + type(scalar_field), dimension(sys_size), intent(in) :: q_prim_vf + real(wp), dimension(num_fluids) :: alpha + real(wp), dimension(num_vels) :: vel, vel_hit + real(wp), dimension(2) :: Re + real(wp) :: rho, vel_sum, pres, gamma, pi_inf, qv, c, H + real(wp) :: rho_hit, pres_hit, c_hit + real(wp) :: icfl, vcfl, Rc, ccfl, icfl_hit + integer :: i, j, k, l, fl, j_hit, k_hit, l_hit + real(wp) :: x_hit, y_hit, z_hit, dist + logical :: nan_hit + integer :: near1_id, near2_id + real(wp) :: near1_dist, near2_dist + + do i = 1, sys_size + $:GPU_UPDATE(host='[q_prim_vf(i)%sf(:, :, :)]') + end do + if (ib) then + $:GPU_UPDATE(host='[ib_markers%sf]') + end if + + icfl_hit = -huge(1._wp) + nan_hit = .false. + j_hit = 0; k_hit = 0; l_hit = 0 + + scan: do l = 0, p + do k = 0, n + do j = 0, m + if (ib) then + if (ib_markers%sf(j, k, l) /= 0) cycle + end if + + call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k, l) + call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, alpha, vel_sum, 0._wp, c, qv) + + if (any_non_newtonian) then + Re(1) = 0._wp + do fl = 1, num_fluids + if (is_non_newtonian(fl)) then + Re(1) = Re(1) + alpha(fl)*hb_mu_max(fl) + else + Re(1) = Re(1) + alpha(fl)*fluid_inv_re(fl) + end if + end do + Re(1) = 1._wp/max(Re(1), sgm_eps) + end if + + call s_compute_stability_from_dt(vel, c, rho, Re, j, k, l, icfl, vcfl, Rc, ccfl) + + if (.not. f_approx_equal(icfl, icfl)) then + nan_hit = .true. + j_hit = j; k_hit = k; l_hit = l + rho_hit = rho; pres_hit = pres; c_hit = c; vel_hit = vel + exit scan + else if (icfl > icfl_hit) then + icfl_hit = icfl + j_hit = j; k_hit = k; l_hit = l + rho_hit = rho; pres_hit = pres; c_hit = c; vel_hit = vel + end if + end do + end do + end do scan + + x_hit = x_cc(j_hit) + y_hit = 0._wp; if (n > 0) y_hit = y_cc(k_hit) + z_hit = 0._wp; if (p > 0) z_hit = z_cc(l_hit) + + print '(A,I0,A,I0,A,I0,A,I0,A)', 'ICFL violation on rank ', proc_rank, ': cell (j,k,l) = (', j_hit, ',', k_hit, ',', & + & l_hit, ')' + if (nan_hit) then + print '(A)', ' icfl = NaN' + else + print '(A,ES16.6)', ' icfl = ', icfl_hit + end if + print '(A,3(ES16.6,1X))', ' position = ', x_hit, y_hit, z_hit + print '(A,ES16.6,A,ES16.6,A,ES16.6)', ' rho, pres, c = ', rho_hit, ', ', pres_hit, ', ', c_hit + print '(A,3(ES16.6,1X))', ' velocity = ', vel_hit + if (ib) print '(A,I0)', ' ib_markers = ', ib_markers%sf(j_hit, k_hit, l_hit) + + if (ib .and. num_ibs > 0) then + near1_id = 0; near1_dist = huge(1._wp) + near2_id = 0; near2_dist = huge(1._wp) + do i = 1, num_ibs + dist = sqrt((x_hit - patch_ib(i)%x_centroid)**2 + (y_hit - patch_ib(i)%y_centroid)**2 + (z_hit & + & - patch_ib(i)%z_centroid)**2) + if (dist < near1_dist) then + near2_dist = near1_dist; near2_id = near1_id + near1_dist = dist; near1_id = i + else if (dist < near2_dist) then + near2_dist = dist; near2_id = i + end if + end do + if (near1_id > 0) print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' nearest particle id=', near1_id, ' dist=', & + & near1_dist, ' gap=', near1_dist - patch_ib(near1_id)%radius, ' vel=', patch_ib(near1_id)%vel + if (near2_id > 0) print '(A,I0,A,ES16.6,A,ES16.6)', ' 2nd nearest particle id=', near2_id, ' dist=', near2_dist, & + & ' gap=', near2_dist - patch_ib(near2_id)%radius + end if + + call flush (6) + + end subroutine s_report_icfl_violation + !> Write grid and conservative variable data files in serial format impure subroutine s_write_serial_data_files(q_cons_vf, q_T_sf, q_prim_vf, t_step, bc_type, beta) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 6e0843352..0911f5d2d 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -18,7 +18,7 @@ module m_particle_cloud private - public :: s_generate_particle_clouds + public :: s_generate_particle_clouds, s_add_cloud_particle contains diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 4d10417ed..760c02d11 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -907,12 +907,10 @@ contains 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 + call s_restart_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) 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 + call s_restart_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) else call s_generate_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) end if @@ -1142,8 +1140,14 @@ 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 + !> @brief Reads IB kinematic state from restart_data/ib_state.dat on restart. Rank 0 reads the last num_gbl_ibs records and !! broadcasts to all ranks. Overwrites patch_ib vel, angular_vel, angles, and centroid. + !! + !! num_ibs is still just the namelist-declared patch count at this point (particle-cloud patches haven't + !! been merged into patch_ib yet), so num_gbl_ibs is (re)computed directly from particle_cloud(:) here + !! rather than trusting the stale global - this is also what bounds the read/broadcast loops below. + !! s_restart_particle_clouds subsequently harvests the kinematic state written here into a + !! particle_cloud_ibs array for s_reduce_ib_patch_array, without re-running particle placement. impure subroutine s_read_ib_restart_data(t_step) integer, intent(in) :: t_step @@ -1155,6 +1159,11 @@ contains logical :: file_exist character(len=10) :: t_step_string + num_gbl_ibs = num_ibs + do i = 1, num_particle_clouds + num_gbl_ibs = num_gbl_ibs + particle_cloud(i)%num_particles + end do + if (file_per_process) then call s_int_to_str(t_step, t_step_string) @@ -1200,7 +1209,7 @@ contains 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)) - do i = 1, num_ibs + do i = 1, num_gbl_ibs read (file_unit, iostat=ios) ib_buf if (ios /= 0) call s_mpi_abort('Error reading IB state restart file') @@ -1216,7 +1225,7 @@ contains end if #ifdef MFC_MPI - do i = 1, num_ibs + do i = 1, num_gbl_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) @@ -1229,6 +1238,44 @@ contains end subroutine s_read_ib_restart_data + !> @brief Rebuilds particle_cloud_ibs after a restart from the kinematic state s_read_ib_restart_data just wrote into patch_ib + !! (keyed by gbl_patch_id), instead of re-running particle placement. Geometry (radius, mass, moving_ibm, etc) is uniform per + !! cloud and reconstructed directly from the particle_cloud(:) case parameters via s_add_cloud_particle - only + !! position/velocity/angles need to come from the checkpoint. The result feeds into the normal s_reduce_ib_patch_array call, so + !! neighborhood/local ownership is decided the usual way, just using the restart positions instead of freshly-packed ones. + impure subroutine s_restart_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) + + type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs + integer, intent(out) :: num_particle_cloud_ibs + integer :: cloud_idx, i, ib_idx, glbl_idx, geom + real(wp), dimension(3) :: vel, angular_vel, angles + + geom = merge(2, 8, num_dims < 3) ! circle for 2D, sphere for 3D - matches s_particle_cloud_rejection_pack/lattice + + num_particle_cloud_ibs = 0 + do cloud_idx = 1, num_particle_clouds + num_particle_cloud_ibs = num_particle_cloud_ibs + particle_cloud(cloud_idx)%num_particles + end do + allocate (particle_cloud_ibs(num_particle_cloud_ibs)) + + ib_idx = 0 + glbl_idx = num_ibs + do cloud_idx = 1, num_particle_clouds + do i = 1, particle_cloud(cloud_idx)%num_particles + glbl_idx = glbl_idx + 1 + vel = patch_ib(glbl_idx)%vel + angular_vel = patch_ib(glbl_idx)%angular_vel + angles = patch_ib(glbl_idx)%angles + call s_add_cloud_particle(cloud_idx, ib_idx, glbl_idx, geom, patch_ib(glbl_idx)%x_centroid, & + & patch_ib(glbl_idx)%y_centroid, patch_ib(glbl_idx)%z_centroid, particle_cloud_ibs) + particle_cloud_ibs(ib_idx)%vel = vel + particle_cloud_ibs(ib_idx)%angular_vel = angular_vel + particle_cloud_ibs(ib_idx)%angles = angles + end do + end do + + end subroutine s_restart_particle_clouds + !> @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 From 868e439795ea3ee09d0d6355d73d7fa81282de3b Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Thu, 27 Aug 2026 22:48:02 -0400 Subject: [PATCH 02/10] intermittent commit --- src/common/m_constants.fpp | 4 +-- src/common/m_derived_types.fpp | 1 + src/simulation/m_global_parameters.fpp | 1 + src/simulation/m_ibm.fpp | 2 ++ src/simulation/m_mpi_proxy.fpp | 1 + src/simulation/m_particle_cloud.fpp | 38 ++++++++++++++++++++------ src/simulation/m_start_up.fpp | 20 +++++++++----- 7 files changed, 49 insertions(+), 18 deletions(-) diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 8fc53f1ac..38cc3ae49 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 = 2000 !< Maximum number of immersed boundary patches (patch_ib) - integer, parameter :: num_ib_patches_max_namelist = 54000 + 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_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 2fd61794b..05d42b3bc 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -379,6 +379,7 @@ module m_derived_types integer :: moving_ibm !< Motion flag: 0=static, 1=moving (forces), 2=forced path integer :: seed !< Random seed for reproducible placement integer :: cloud_geometry !< Cloud region geometry: 1=box, 2=hemisphere shell + integer :: shell_axis !< Axis the hemisphere shell opens toward: 1=x, 2=y, 3=z (2D ignores 3) integer :: packing_method !< Packing algorithm: 1=rejection sampling, 2=lattice integer :: periodic !< Periodic overlap flag for box rejection packing: 0=off, 1=on end type particle_cloud_parameters diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 8bfc807d7..1a3ce704f 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -644,6 +644,7 @@ contains particle_cloud(i)%moving_ibm = 0 particle_cloud(i)%seed = 0 particle_cloud(i)%cloud_geometry = 1 + particle_cloud(i)%shell_axis = 3 particle_cloud(i)%packing_method = dflt_int particle_cloud(i)%periodic = 0 end do diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 7d54dba97..4b6003ff3 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1443,6 +1443,8 @@ contains ! check if in local domain if (f_local_rank_owns_location(centroid)) then local_output_idx = local_output_idx + 1 + @: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 end if end if diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index f8fe85a9b..1c566ee86 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -207,6 +207,7 @@ contains 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)%shell_axis, 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 diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 0911f5d2d..2fb07b1ee 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -182,8 +182,9 @@ contains !> Draws one rejection-sampling candidate centre (rx, ry, rz) for cloud_idx, advancing seed in place. For box geometry the !! candidate is uniform in the box and never rejected. For a hemisphere shell the candidate is uniform in the shell volume - 2D !! uses theta uniform on [0, pi] with the sqrt radial CDF; 3D uses uniform phi, uniform cos(polar) on [0, 1], and the cube-root - !! radial CDF - and reject is set when it lands within one particle radius of the flat face (the plane at y_centroid in 2D, - !! z_centroid in 3D), a hard geometric cut applied after sampling that preserves uniformity over the remaining region. + !! radial CDF - and reject is set when it lands within one particle radius of the flat face (the plane through the centroid + !! perpendicular to shell_axis: 1=x, 2=y, 3=z; 2D has no z-axis so any value other than 1 falls back to y, matching the + !! pre-shell_axis default), a hard geometric cut applied after sampling that preserves uniformity over the remaining region. subroutine s_sample_cloud_candidate(cloud_idx, seed, rx, ry, rz, reject) integer, intent(in) :: cloud_idx @@ -219,20 +220,39 @@ contains theta = pi*f_xorshift(seed) u = f_xorshift(seed) r_shell = sqrt((r_outer**2 - r_inner**2)*u + r_inner**2) - rx = particle_cloud(cloud_idx)%x_centroid + r_shell*cos(theta) - ry = particle_cloud(cloud_idx)%y_centroid + r_shell*sin(theta) rz = particle_cloud(cloud_idx)%z_centroid - if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + if (particle_cloud(cloud_idx)%shell_axis == 1) then + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*sin(theta) + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*cos(theta) + if (rx < particle_cloud(cloud_idx)%x_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + else + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*cos(theta) + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*sin(theta) + if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + end if else phi = 2._wp*pi*f_xorshift(seed) zdir = f_xorshift(seed) rho = sqrt(max(0._wp, 1._wp - zdir**2)) u = f_xorshift(seed) r_shell = ((r_outer**3 - r_inner**3)*u + r_inner**3)**(1._wp/3._wp) - rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi) - ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*sin(phi) - rz = particle_cloud(cloud_idx)%z_centroid + r_shell*zdir - if (rz < particle_cloud(cloud_idx)%z_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + select case (particle_cloud(cloud_idx)%shell_axis) + case (1) ! opens toward +x + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*zdir + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*cos(phi) + rz = particle_cloud(cloud_idx)%z_centroid + r_shell*rho*sin(phi) + if (rx < particle_cloud(cloud_idx)%x_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + case (2) ! opens toward +y + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*zdir + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi) + rz = particle_cloud(cloud_idx)%z_centroid + r_shell*rho*sin(phi) + if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + case default ! 3: opens toward +z + rz = particle_cloud(cloud_idx)%z_centroid + r_shell*zdir + rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi) + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*sin(phi) + if (rz < particle_cloud(cloud_idx)%z_centroid + particle_cloud(cloud_idx)%radius) reject = .true. + end select end if case default call s_mpi_abort("Particle cloud geometry is not a known cloud geometry of MFC. Exiting.") diff --git a/src/simulation/m_start_up.fpp b/src/simulation/m_start_up.fpp index 760c02d11..fb941098b 100644 --- a/src/simulation/m_start_up.fpp +++ b/src/simulation/m_start_up.fpp @@ -1241,28 +1241,33 @@ contains !> @brief Rebuilds particle_cloud_ibs after a restart from the kinematic state s_read_ib_restart_data just wrote into patch_ib !! (keyed by gbl_patch_id), instead of re-running particle placement. Geometry (radius, mass, moving_ibm, etc) is uniform per !! cloud and reconstructed directly from the particle_cloud(:) case parameters via s_add_cloud_particle - only - !! position/velocity/angles need to come from the checkpoint. The result feeds into the normal s_reduce_ib_patch_array call, so - !! neighborhood/local ownership is decided the usual way, just using the restart positions instead of freshly-packed ones. + !! position/velocity/angles need to come from the checkpoint. Every rank walks the full global particle list but only keeps + !! entries this rank's IB neighborhood owns (f_neighborhood_ranks_own_location, same test s_reduce_particle_cloud_ibs applies + !! after fresh packing) - s_reduce_ib_patch_array trusts particle_cloud_ibs to already be neighborhood-filtered and does not + !! recheck it, so skipping this filter here would silently inflate num_ibs on every rank to the full global particle count. impure subroutine s_restart_particle_clouds(particle_cloud_ibs, num_particle_cloud_ibs) type(ib_patch_parameters), allocatable, intent(out), dimension(:) :: particle_cloud_ibs integer, intent(out) :: num_particle_cloud_ibs - integer :: cloud_idx, i, ib_idx, glbl_idx, geom - real(wp), dimension(3) :: vel, angular_vel, angles + integer :: cloud_idx, i, ib_idx, glbl_idx, geom, n_total_particles + real(wp), dimension(3) :: vel, angular_vel, angles, centroid geom = merge(2, 8, num_dims < 3) ! circle for 2D, sphere for 3D - matches s_particle_cloud_rejection_pack/lattice - num_particle_cloud_ibs = 0 + n_total_particles = 0 do cloud_idx = 1, num_particle_clouds - num_particle_cloud_ibs = num_particle_cloud_ibs + particle_cloud(cloud_idx)%num_particles + n_total_particles = n_total_particles + particle_cloud(cloud_idx)%num_particles end do - allocate (particle_cloud_ibs(num_particle_cloud_ibs)) + allocate (particle_cloud_ibs(min(num_ib_patches_max_namelist, n_total_particles))) ib_idx = 0 glbl_idx = num_ibs do cloud_idx = 1, num_particle_clouds do i = 1, particle_cloud(cloud_idx)%num_particles glbl_idx = glbl_idx + 1 + centroid = [patch_ib(glbl_idx)%x_centroid, patch_ib(glbl_idx)%y_centroid, 0._wp] + if (num_dims == 3) centroid(3) = patch_ib(glbl_idx)%z_centroid + if (.not. f_neighborhood_ranks_own_location(centroid)) cycle vel = patch_ib(glbl_idx)%vel angular_vel = patch_ib(glbl_idx)%angular_vel angles = patch_ib(glbl_idx)%angles @@ -1273,6 +1278,7 @@ contains particle_cloud_ibs(ib_idx)%angles = angles end do end do + num_particle_cloud_ibs = ib_idx end subroutine s_restart_particle_clouds From 6fb3e87a226ba4b15049143413853ef472721f94 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Thu, 27 Aug 2026 23:43:39 -0400 Subject: [PATCH 03/10] Resolved race condition on ghost point correction --- src/simulation/m_ibm.fpp | 94 +++++++++++++++++++++++++++++++++------- 1 file changed, 79 insertions(+), 15 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 4b6003ff3..a7f43bfd2 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -184,6 +184,18 @@ contains type(ghost_point) :: gp type(ghost_point) :: innerp + ! Per-ghost-point image-point interpolation results, stashed between the interpolation + ! kernel and the correction kernel below so that the correction kernel never reads + ! q_prim_vf (or pb_in/mv_in) at a cell another ghost point's correction may have already + ! overwritten this stage - i.e. so the two kernels never race on those shared fields. + real(wp), allocatable :: alpha_rho_IP_buf(:,:), alpha_IP_buf(:,:) + real(wp), allocatable :: pres_IP_buf(:), c_IP_buf(:) + real(wp), allocatable :: vel_IP_buf(:,:) + real(wp), allocatable :: r_IP_buf(:,:), v_IP_buf(:,:), pb_IP_buf(:,:), mv_IP_buf(:,:) + real(wp), allocatable :: nmom_IP_buf(:,:) + real(wp), allocatable :: presb_IP_buf(:,:), massv_IP_buf(:,:) + real(wp), allocatable :: Ys_IP_buf(:,:) + ! set the Moving IBM interior conservative variables $:GPU_PARALLEL_LOOP(private='[i, j, k, patch_id, rho]', collapse=3) do l = 0, p @@ -219,23 +231,22 @@ contains $:END_GPU_PARALLEL_LOOP() if (num_gps > 0) then - $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, & - & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & - & innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id, & - & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff]') + @:ALLOCATE(alpha_rho_IP_buf(1:num_fluids, 1:num_gps), alpha_IP_buf(1:num_fluids, 1:num_gps), & + & pres_IP_buf(1:num_gps), c_IP_buf(1:num_gps), vel_IP_buf(1:3, 1:num_gps), & + & r_IP_buf(1:nb, 1:num_gps), v_IP_buf(1:nb, 1:num_gps), pb_IP_buf(1:nb, 1:num_gps), & + & mv_IP_buf(1:nb, 1:num_gps), nmom_IP_buf(1:nb*nmom, 1:num_gps), & + & presb_IP_buf(1:nb*nnode, 1:num_gps), massv_IP_buf(1:nb*nnode, 1:num_gps), & + & Ys_IP_buf(1:num_species, 1:num_gps)) + + ! Phase 1: interpolate image-point primitives for every ghost point from the pre-correction field only. Kept in its + ! own kernel (rather than fused with phase 2 below) so the kernel-launch boundary between them guarantees every + ! interpolation here happens-before any q_prim_vf/q_cons_vf write in phase 2 - otherwise, with densely-packed ghost + ! regions, one ghost point's image-point stencil can land on a cell that is itself another ghost point being + ! corrected in the same parallel loop, racing the read against that write. + $:GPU_PARALLEL_LOOP(private='[i, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP, r_IP, v_IP, pb_IP, mv_IP, & + & nmom_IP, presb_IP, massv_IP, Ys_IP]') do i = 1, num_gps gp = ghost_points(i) - j = gp%loc(1) - k = gp%loc(2) - l = gp%loc(3) - patch_id = ghost_points(i)%ib_patch_id - - ! Calculate physical location of GP - if (p > 0) then - physical_loc = [x_cc(j), y_cc(k), z_cc(l)] - else - physical_loc = [x_cc(j), y_cc(k), 0._wp] - end if ! Interpolate primitive variables at image point associated w/ GP if (bubbles_euler .and. .not. qbmm) then @@ -253,6 +264,56 @@ contains call s_interpolate_image_point(q_prim_vf, gp, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, c_IP) end if + alpha_rho_IP_buf(:, i) = alpha_rho_IP(1:num_fluids) + alpha_IP_buf(:, i) = alpha_IP(1:num_fluids) + pres_IP_buf(i) = pres_IP + vel_IP_buf(:, i) = vel_IP + c_IP_buf(i) = c_IP + r_IP_buf(:, i) = r_IP(1:nb) + v_IP_buf(:, i) = v_IP(1:nb) + pb_IP_buf(:, i) = pb_IP(1:nb) + mv_IP_buf(:, i) = mv_IP(1:nb) + nmom_IP_buf(:, i) = nmom_IP(1:nb*nmom) + presb_IP_buf(:, i) = presb_IP(1:nb*nnode) + massv_IP_buf(:, i) = massv_IP(1:nb*nnode) + Ys_IP_buf(:, i) = Ys_IP(1:num_species) + end do + $:END_GPU_PARALLEL_LOOP() + + ! Phase 2: apply the buffered image-point results as ghost-point corrections to q_prim_vf/q_cons_vf. + $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, & + & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & + & innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id, & + & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff]') + do i = 1, num_gps + gp = ghost_points(i) + j = gp%loc(1) + k = gp%loc(2) + l = gp%loc(3) + patch_id = ghost_points(i)%ib_patch_id + + ! Calculate physical location of GP + if (p > 0) then + physical_loc = [x_cc(j), y_cc(k), z_cc(l)] + else + physical_loc = [x_cc(j), y_cc(k), 0._wp] + end if + + ! Recover the image-point interpolation computed for this ghost point in phase 1 + alpha_rho_IP(1:num_fluids) = alpha_rho_IP_buf(:, i) + alpha_IP(1:num_fluids) = alpha_IP_buf(:, i) + pres_IP = pres_IP_buf(i) + vel_IP = vel_IP_buf(:, i) + c_IP = c_IP_buf(i) + r_IP(1:nb) = r_IP_buf(:, i) + v_IP(1:nb) = v_IP_buf(:, i) + pb_IP(1:nb) = pb_IP_buf(:, i) + mv_IP(1:nb) = mv_IP_buf(:, i) + nmom_IP(1:nb*nmom) = nmom_IP_buf(:, i) + presb_IP(1:nb*nnode) = presb_IP_buf(:, i) + massv_IP(1:nb*nnode) = massv_IP_buf(:, i) + Ys_IP(1:num_species) = Ys_IP_buf(:, i) + ! Injecting (burning) surface: replace the mirrored ghost composition with pure ! injected fuel at the local pressure and the ambient (image-point) temperature. ! Setting a consistent injected density here (rather than reusing the heavy ambient @@ -449,6 +510,9 @@ contains end if end do $:END_GPU_PARALLEL_LOOP() + + @:DEALLOCATE(alpha_rho_IP_buf, alpha_IP_buf, pres_IP_buf, c_IP_buf, vel_IP_buf, r_IP_buf, v_IP_buf, pb_IP_buf, & + & mv_IP_buf, nmom_IP_buf, presb_IP_buf, massv_IP_buf, Ys_IP_buf) end if end subroutine s_ibm_correct_state From 2198d0e245878d21cbe45b38a94b55b40c40b84d Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Fri, 28 Aug 2026 02:02:42 -0400 Subject: [PATCH 04/10] Found a bug I labeled with a TODO and more improvements to ownership handoff, which ic clearly the broken function --- src/simulation/m_ib_patches.fpp | 67 +++++++++++++++++++++++---------- src/simulation/m_ibm.fpp | 2 + 2 files changed, 50 insertions(+), 19 deletions(-) diff --git a/src/simulation/m_ib_patches.fpp b/src/simulation/m_ib_patches.fpp index 0a845c8fe..ab291f736 100644 --- a/src/simulation/m_ib_patches.fpp +++ b/src/simulation/m_ib_patches.fpp @@ -228,38 +228,53 @@ contains ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - ! perform the interior check for the patch geometry of this IB + ! perform the interior check for the patch geometry of this IB. Writes to ib_markers use + ! an atomic max (not a plain assignment) because this loop is parallel over patch_id: the + ! soft-sphere collision model allows particles to physically interpenetrate by design, so + ! two different patches can both claim the same overlapping cell here. A plain write would + ! be a data race with a nondeterministic winner; the atomic max makes the higher + ! encoded_patch_id win consistently every time, regardless of thread scheduling. if (patch_ib(patch_id)%geometry == 8) then ! sphere geometry radius = patch_ib(patch_id)%radius - if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), & - & radius)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_sphere(xyz_local(1), xyz_local(2), xyz_local(3), radius)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 9) then ! cuboid geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, & & patch_ib(patch_id)%length_z] - if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), & - & length)) ib_markers%sf(i, j, k) = encoded_patch_id + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 10) then ! cylinder geometry radius = patch_ib(patch_id)%radius if (f_is_inside_cylinder(xyz_local(2), xyz_local(3), xyz_local(1), radius, & - & patch_ib(patch_id)%length_x)) ib_markers%sf(i, j, k) = encoded_patch_id + & patch_ib(patch_id)%length_x)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 11) then ! 3D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), xyz_local(3), & - & patch_ib(patch_id)%length_z, airfoil_id)) ib_markers%sf(i, j, & - & k) = encoded_patch_id + & patch_ib(patch_id)%length_z, airfoil_id)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 12) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset model_id = patch_ib(patch_id)%model_id eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local) if (eta > stl_models(model_id)%model_threshold) then - ib_markers%sf(i, j, k) = encoded_patch_id + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, k) = max(ib_markers%sf(i, j, k), encoded_patch_id) end if end if end do @@ -296,36 +311,50 @@ contains ! rotate the frame into the IB's coordinates xyz_local = matmul(patch_ib(patch_id)%rotation_matrix_inverse, xyz_local) - ! perform the interior check for the patch geometry of this IB + ! perform the interior check for the patch geometry of this IB. Writes to ib_markers use an + ! atomic max (not a plain assignment) because this loop is parallel over patch_id: the + ! soft-sphere collision model allows particles to physically interpenetrate by design, so two + ! different patches can both claim the same overlapping cell here. A plain write would be a data + ! race with a nondeterministic winner; the atomic max makes the higher encoded_patch_id win + ! consistently every time, regardless of thread scheduling. if (patch_ib(patch_id)%geometry == 2) then ! circular geometries radius = patch_ib(patch_id)%radius - if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) ib_markers%sf(i, & - & j, 0) = encoded_patch_id + if (f_is_inside_cylinder(xyz_local(1), xyz_local(2), 0._wp, radius, 0._wp)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 3) then ! rectangular geometries length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) ib_markers%sf(i, j, & - & 0) = encoded_patch_id + if (f_is_inside_cuboid(xyz_local(1), xyz_local(2), xyz_local(3), length)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 4) then ! 2D airfoil geometry airfoil_id = patch_ib(patch_id)%airfoil_id xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset - if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, & - & airfoil_id)) ib_markers%sf(i, j, 0) = encoded_patch_id + if (f_is_inside_airfoil(xyz_local(1), xyz_local(2), 0._wp, 0._wp, airfoil_id)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if else if (patch_ib(patch_id)%geometry == 5) then ! STL model geometry xyz_local = xyz_local - patch_ib(patch_id)%centroid_offset model_id = patch_ib(patch_id)%model_id eta = f_model_is_inside(gpu_ntrs(model_id), model_id, xyz_local) if (eta > stl_models(model_id)%model_threshold) then - ib_markers%sf(i, j, 0) = encoded_patch_id + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) end if else if (patch_ib(patch_id)%geometry == 6) then ! ellipse geometry length = [patch_ib(patch_id)%length_x, patch_ib(patch_id)%length_y, 0._wp] - if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) ib_markers%sf(i, j, & - & 0) = encoded_patch_id + if (f_is_inside_ellipse(xyz_local(1), xyz_local(2), length)) then + $:GPU_ATOMIC(atomic='update') + ib_markers%sf(i, j, 0) = max(ib_markers%sf(i, j, 0), encoded_patch_id) + end if end if end do end do diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index a7f43bfd2..6dfa10cda 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1550,6 +1550,7 @@ contains pack_pos = storage_size(0)/8 + new_count*patch_bytes ! Post all receives first, then sends + ! TODO :: THIS NEEDS TO ITERATE OVER -ib_neighborhood_radius to ib_neighborhood_radius, not -1 to 1 nreqs = 0 nbr_idx = 0 do dz = merge(-1, 0, num_dims == 3), merge(1, 0, num_dims == 3) @@ -1597,6 +1598,7 @@ contains num_ibs = num_ibs + 1 @:ASSERT(num_ibs <= size(patch_ib), 'patch_ib overflow in neighborhood handoff') patch_ib(num_ibs) = tmp_patch + ib_gbl_idx_lookup(tmp_patch%gbl_patch_id) = num_ibs end if end do end do From 627873ff6b8d45d1416f92b4072114320cc76107 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Fri, 28 Aug 2026 11:50:26 -0400 Subject: [PATCH 05/10] Intermittent commit --- src/simulation/m_ibm.fpp | 235 +++++++++++++++++++++++++---- src/simulation/m_time_steppers.fpp | 1 + 2 files changed, 208 insertions(+), 28 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 6dfa10cda..6109a702b 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -1166,10 +1166,12 @@ contains end do ! apply the summed forces - $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') do i = 1, num_ibs - patch_ib(i)%force(:) = forces(i,:) - patch_ib(i)%torque(:) = torques(i,:) + do l = 1, 3 + patch_ib(i)%force(l) = forces(i,l) + patch_ib(i)%torque(l) = torques(i,l) + end do end do $:END_GPU_PARALLEL_LOOP() @@ -1353,7 +1355,7 @@ contains real(wp), dimension(num_ibs, 3), intent(inout) :: forces, torques #ifdef MFC_MPI - integer :: i, j, k, pack_pos, unpack_pos, buf_size, ierr + 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(:) @@ -1375,11 +1377,13 @@ contains 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 - $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') do i = 1, num_ibs send_ids(i) = patch_ib(i)%gbl_patch_id - send_ft(1:3,i) = forces(i,:) - send_ft(4:6,i) = torques(i,:) + 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() $:GPU_UPDATE(host='[send_ids, send_ft]') @@ -1395,16 +1399,18 @@ contains 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_PARALLEL_LOOP(private='[i, j]', copyin='[recv_ft, recv_ids]', copy='[forces, torques, & + $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques, & & recv_forces_snap, recv_torques_snap]') 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 - forces(j,:) = forces(j,:) + recv_ft(1:3,i) - recv_forces_snap(j,:) - torques(j,:) = torques(j,:) + recv_ft(4:6,i) - recv_torques_snap(j,:) - recv_forces_snap(j,:) = recv_ft(1:3,i) - recv_torques_snap(j,:) = recv_ft(4:6,i) + 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() @@ -1422,11 +1428,13 @@ contains do k = 1, min(2*ib_neighborhood_radius, num_procs_${X}$ - 1) pack_pos = 0 - $:GPU_PARALLEL_LOOP(private='[i]', copyin='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') do i = 1, num_ibs send_ids(i) = patch_ib(i)%gbl_patch_id - send_ft(1:3,i) = forces(i,:) - send_ft(4:6,i) = torques(i,:) + 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() $:GPU_UPDATE(host='[send_ids, send_ft]') @@ -1441,12 +1449,14 @@ contains 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_PARALLEL_LOOP(private='[i, j]', copyin='[recv_ft, recv_ids]', copy='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques]') do i = 1, recv_count call s_get_neighborhood_idx(recv_ids(i), j) if (j > 0) then - forces(j,:) = recv_ft(1:3,i) - torques(j,:) = recv_ft(4:6,i) + 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() @@ -1550,12 +1560,11 @@ contains pack_pos = storage_size(0)/8 + new_count*patch_bytes ! Post all receives first, then sends - ! TODO :: THIS NEEDS TO ITERATE OVER -ib_neighborhood_radius to ib_neighborhood_radius, not -1 to 1 nreqs = 0 nbr_idx = 0 - do dz = merge(-1, 0, num_dims == 3), merge(1, 0, num_dims == 3) - do dy = -1, 1 - do dx = -1, 1 + do dz = merge(-ib_neighborhood_radius, 0, num_dims == 3), merge(ib_neighborhood_radius, 0, num_dims == 3) + do dy = -ib_neighborhood_radius, ib_neighborhood_radius + do dx = -ib_neighborhood_radius, ib_neighborhood_radius if (dx == 0 .and. dy == 0 .and. dz == 0) cycle nbr_idx = nbr_idx + 1 tag = 200 + (dx + 1)*9 + (dy + 1)*3 + (dz + 1) @@ -1570,9 +1579,9 @@ contains end do end do - do dz = merge(-1, 0, num_dims == 3), merge(1, 0, num_dims == 3) - do dy = -1, 1 - do dx = -1, 1 + do dz = merge(-ib_neighborhood_radius, 0, num_dims == 3), merge(ib_neighborhood_radius, 0, num_dims == 3) + do dy = -ib_neighborhood_radius, ib_neighborhood_radius + do dx = -ib_neighborhood_radius, ib_neighborhood_radius if (dx == 0 .and. dy == 0 .and. dz == 0) cycle tag = 200 + (dx + 1)*9 + (dy + 1)*3 + (dz + 1) send_neighbor = ib_neighbor_ranks(dx, dy, dz) @@ -1586,7 +1595,7 @@ contains call MPI_WAITALL(nreqs, requests, MPI_STATUSES_IGNORE, ierr) ! Unpack all received buffers - do nbr_idx = 1, merge(26, 8, num_dims == 3) + do nbr_idx = 1, ((2*ib_neighborhood_radius+1)**num_dims) - 1 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) @@ -1611,6 +1620,171 @@ contains end subroutine s_handoff_ib_ownership + !> TEMPORARY DEBUG INSTRUMENTATION (remove once the cross-rank IB divergence bug is found). Gathers every rank's tracked + !! patch_ib states and, for every gbl_patch_id two or more ranks both track, logs a full side-by-side dump to + !! ib_divergence_rank.log whenever their dynamic (kinematic/force) fields disagree at all. + subroutine s_debug_log_ib_divergence(t_step) + + integer, intent(in) :: t_step + +#ifdef MFC_MPI + integer :: ierr, patch_bytes, i, j, r, unpack_pos, unit_num + integer, dimension(0:num_procs - 1) :: rank_counts, rank_counts_bytes, rank_displs_bytes + character(len=1), allocatable :: send_buf(:), recv_buf(:) + type(ib_patch_parameters) :: other_patch + character(len=64) :: fname + logical :: mismatch + integer, dimension(9) :: topo_local + integer, dimension(9, 0:num_procs - 1) :: topo_all + integer :: r2, jr, roster_unpack_pos + logical :: found + type(ib_patch_parameters) :: roster_patch + + if (num_procs == 1) return + + ! Gather each rank's Cartesian coords + flow-field boundary neighbor ranks (bc_x/y/z%beg/end) so a divergence dump can + ! show the real adjacency graph instead of an assumed one. + topo_local = -999 + topo_local(1:num_dims) = proc_coords(1:num_dims) + topo_local(4) = bc_x%beg; topo_local(5) = bc_x%end + if (num_dims >= 2) then + topo_local(6) = bc_y%beg; topo_local(7) = bc_y%end + end if + if (num_dims >= 3) then + topo_local(8) = bc_z%beg; topo_local(9) = bc_z%end + end if + call MPI_ALLGATHER(topo_local, 9, MPI_INTEGER, topo_all, 9, MPI_INTEGER, MPI_COMM_WORLD, ierr) + + call MPI_ALLGATHER(num_ibs, 1, MPI_INTEGER, rank_counts, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) + + patch_bytes = storage_size(patch_ib(1))/8 + rank_counts_bytes = rank_counts*patch_bytes + rank_displs_bytes(0) = 0 + do r = 1, num_procs - 1 + rank_displs_bytes(r) = rank_displs_bytes(r - 1) + rank_counts_bytes(r - 1) + end do + + allocate (send_buf(max(1, num_ibs*patch_bytes))) + allocate (recv_buf(max(1, sum(rank_counts_bytes)))) + + unpack_pos = 0 + do i = 1, num_ibs + call MPI_PACK(patch_ib(i), patch_bytes, MPI_BYTE, send_buf, size(send_buf), unpack_pos, MPI_COMM_WORLD, ierr) + end do + + call MPI_ALLGATHERV(send_buf, num_ibs*patch_bytes, MPI_PACKED, recv_buf, rank_counts_bytes, rank_displs_bytes, & + & MPI_PACKED, MPI_COMM_WORLD, ierr) + + write (fname, '(A,I0,A)') 'ib_divergence_rank', proc_rank, '.log' + unit_num = 900 + proc_rank + + do r = 0, num_procs - 1 + if (r == proc_rank) cycle + unpack_pos = rank_displs_bytes(r) + do j = 1, rank_counts(r) + call MPI_UNPACK(recv_buf, size(recv_buf), unpack_pos, other_patch, patch_bytes, MPI_BYTE, MPI_COMM_WORLD, ierr) + + do i = 1, num_ibs + if (patch_ib(i)%gbl_patch_id /= other_patch%gbl_patch_id) cycle + + mismatch = (patch_ib(i)%x_centroid /= other_patch%x_centroid) .or. & + & (patch_ib(i)%y_centroid /= other_patch%y_centroid) .or. & + & (patch_ib(i)%z_centroid /= other_patch%z_centroid) .or. & + & any(patch_ib(i)%vel /= other_patch%vel) .or. & + & any(patch_ib(i)%angular_vel /= other_patch%angular_vel) .or. & + & any(patch_ib(i)%angles /= other_patch%angles) .or. & + & any(patch_ib(i)%force /= other_patch%force) .or. & + & any(patch_ib(i)%torque /= other_patch%torque) .or. & + & (patch_ib(i)%step_x_centroid /= other_patch%step_x_centroid) .or. & + & (patch_ib(i)%step_y_centroid /= other_patch%step_y_centroid) .or. & + & (patch_ib(i)%step_z_centroid /= other_patch%step_z_centroid) .or. & + & any(patch_ib(i)%step_vel /= other_patch%step_vel) .or. & + & any(patch_ib(i)%step_angular_vel /= other_patch%step_angular_vel) .or. & + & any(patch_ib(i)%step_angles /= other_patch%step_angles) + + if (mismatch) then + open (unit=unit_num, file=fname, position='append', action='write', status='unknown') + write (unit_num, '(A)') '====================================================================' + write (unit_num, '(A,I0,A,ES23.15,A,I0,A,I0,A,I0)') 'DIVERGENCE t_step=', t_step, ' mytime=', mytime, & + & ' gbl_patch_id=', patch_ib(i)%gbl_patch_id, ' rank_A=', proc_rank, ' rank_B=', r + write (unit_num, '(A,I0,A,I0,A,I0,A,I0,A,I0)') 'num_procs_x=', num_procs_x, ' num_procs_y=', & + & num_procs_y, ' num_procs_z=', num_procs_z, ' ib_neighborhood_radius=', ib_neighborhood_radius, & + & ' num_dims=', num_dims + write (unit_num, '(A,I0,A,3I3,A,4I5,A,2I5)') '--- rank ', proc_rank, & + & ' coords=', topo_all(1:3, proc_rank), ' bc_x(beg,end)/bc_y(beg,end)=', & + & topo_all(4:5, proc_rank), topo_all(6:7, proc_rank), ' bc_z(beg,end)=', topo_all(8:9, proc_rank) + call s_debug_write_ib_state(unit_num, patch_ib(i)) + write (unit_num, '(A,I0,A,3I3,A,4I5,A,2I5)') '--- rank ', r, & + & ' coords=', topo_all(1:3, r), ' bc_x(beg,end)/bc_y(beg,end)=', & + & topo_all(4:5, r), topo_all(6:7, r), ' bc_z(beg,end)=', topo_all(8:9, r) + call s_debug_write_ib_state(unit_num, other_patch) + + ! Full roster: every rank's tracking status for this gbl_patch_id, not just the mismatching pair, so a + ! rank silently sending 0s (or not tracking at all) shows up instead of being inferred from absence. + write (unit_num, '(A)') '--- full roster for this gbl_patch_id ---' + do r2 = 0, num_procs - 1 + if (r2 == proc_rank) then + found = .false. + do jr = 1, num_ibs + if (patch_ib(jr)%gbl_patch_id == patch_ib(i)%gbl_patch_id) then + found = .true. + write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, ' TRACKS coords=', topo_all(1:3, r2) + call s_debug_write_ib_state(unit_num, patch_ib(jr)) + exit + end if + end do + if (.not. found) write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, & + & ' NOT TRACKED coords=', topo_all(1:3, r2) + else + found = .false. + roster_unpack_pos = rank_displs_bytes(r2) + do jr = 1, rank_counts(r2) + call MPI_UNPACK(recv_buf, size(recv_buf), roster_unpack_pos, roster_patch, patch_bytes, & + & MPI_BYTE, MPI_COMM_WORLD, ierr) + if (roster_patch%gbl_patch_id == patch_ib(i)%gbl_patch_id) then + found = .true. + write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, ' TRACKS coords=', topo_all(1:3, r2) + call s_debug_write_ib_state(unit_num, roster_patch) + exit + end if + end do + if (.not. found) write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, & + & ' NOT TRACKED coords=', topo_all(1:3, r2) + end if + end do + + close (unit_num) + end if + end do + end do + end do + + deallocate (send_buf, recv_buf) +#endif + + end subroutine s_debug_log_ib_divergence + + !> TEMPORARY DEBUG INSTRUMENTATION helper for s_debug_log_ib_divergence: dumps every dynamic field of an IB patch state. + subroutine s_debug_write_ib_state(unit_num, patch) + + integer, intent(in) :: unit_num + type(ib_patch_parameters), intent(in) :: patch + + write (unit_num, '(A,3ES23.15)') ' centroid = ', patch%x_centroid, patch%y_centroid, patch%z_centroid + write (unit_num, '(A,3ES23.15)') ' step_centroid = ', patch%step_x_centroid, patch%step_y_centroid, & + & patch%step_z_centroid + write (unit_num, '(A,3ES23.15)') ' vel = ', patch%vel + write (unit_num, '(A,3ES23.15)') ' step_vel = ', patch%step_vel + write (unit_num, '(A,3ES23.15)') ' angular_vel = ', patch%angular_vel + write (unit_num, '(A,3ES23.15)') ' step_angular_vel = ', patch%step_angular_vel + write (unit_num, '(A,3ES23.15)') ' angles = ', patch%angles + write (unit_num, '(A,3ES23.15)') ' step_angles = ', patch%step_angles + write (unit_num, '(A,3ES23.15)') ' force = ', patch%force + write (unit_num, '(A,3ES23.15)') ' torque = ', patch%torque + write (unit_num, '(A,ES23.15)') ' moment = ', patch%moment + + end subroutine s_debug_write_ib_state + subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx) $:GPU_ROUTINE(parallelism='[seq]') @@ -1627,9 +1801,14 @@ contains integer :: i - ib_gbl_idx_lookup = -1 - $:GPU_UPDATE(device='[ib_gbl_idx_lookup]') + ! reset the lookup + $:GPU_PARALLEL_LOOP(private='[i]') + do i = 1, num_gbl_ibs + ib_gbl_idx_lookup(i) = -1 + end do + $:END_GPU_PARALLEL_LOOP() + ! populate the table $:GPU_PARALLEL_LOOP(private='[i]') do i = 1, num_ibs ib_gbl_idx_lookup(patch_ib(i)%gbl_patch_id) = i diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index daf1c2a73..5583872a0 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -589,6 +589,7 @@ contains if (moving_immersed_boundary_flag) then call s_wrap_periodic_ibs() ! wraps the positions of IBs to the local proc call s_handoff_ib_ownership() ! recomputes which ranks own which IBs and communicate to neighbors + call s_debug_log_ib_divergence(t_step) ! TEMPORARY: logs any cross-rank IB state mismatch else if (ib_state_wrt) then call s_compute_ib_forces(q_prim_vf, fluid_pp) end if From 88877c1249eb22468b57fed2cf29fb5414635d2d Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Sun, 30 Aug 2026 05:03:26 -0400 Subject: [PATCH 06/10] Error found, num_ibs was never updated to GPU, causing out of bounds memory writes which created particle drift --- src/simulation/m_ibm.fpp | 232 +++++++++++++++++++++++++---- src/simulation/m_time_steppers.fpp | 2 + 2 files changed, 209 insertions(+), 25 deletions(-) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 6109a702b..6ee16f8f3 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -48,6 +48,15 @@ module m_ibm real(wp), allocatable :: send_ft(:,:), recv_ft(:,:) real(wp), allocatable :: recv_forces_snap(:,:), recv_torques_snap(:,:) + ! TEMPORARY DEBUG INSTRUMENTATION state (remove once the IB force-communication bug is found): + ! dbg_t_step is set once per timestep by m_time_steppers; dbg_divergence_count/dbg_last_divergent_t_step track how many + ! distinct timesteps s_debug_log_ib_divergence has found a mismatch on, to gate the wire trace and stop the run once + ! two measurements are captured. + integer, save :: dbg_t_step = -1 + integer, save :: dbg_divergence_count = 0 + integer, save :: dbg_last_divergent_t_step = -1 + integer, parameter :: dbg_track_gbl_id = 352 ! gbl_patch_id currently under investigation + contains !> Allocates memory for the variables in the IBM module @@ -1358,9 +1367,20 @@ contains 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(:) + character(len=64) :: fname_trace + integer :: unit_trace + logical :: trace_active if (num_procs == 1) return + ! TEMPORARY DEBUG INSTRUMENTATION (remove once found): trace the wire values for gbl_patch_id=2 at every send/recv, + ! gated to the known failure window only - this does file open/close per match and is too slow to run unconditionally. + write (fname_trace, '(A,I0,A)') 'ib_wire_trace_rank', proc_rank, '.log' + unit_trace = 950 + proc_rank + ! Only trace once the roster check (s_debug_log_ib_divergence) has confirmed at least one divergent timestep, so we + ! don't pay the per-send/recv file I/O cost before we know we're near a real event. + trace_active = (dbg_divergence_count >= 1) + 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)) @@ -1377,7 +1397,10 @@ contains 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 - $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') + ! TEMPORARY: GPU offload disabled on this whole subroutine to test whether it's the source of the y/z + ! force-communication corruption bug (see ib_wire_trace evidence) - these loops are tiny (num_ibs + ! entries), no perf concern. + ! $: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 @@ -1385,8 +1408,19 @@ contains send_ft(l + 3,i) = torques(i,l) end do end do - $:END_GPU_PARALLEL_LOOP() - $:GPU_UPDATE(host='[send_ids, send_ft]') + ! $:END_GPU_PARALLEL_LOOP() + ! $:GPU_UPDATE(host='[send_ids, send_ft]') + if (trace_active) then + do i = 1, num_ibs + if (send_ids(i) == dbg_track_gbl_id) then + open (unit=unit_trace, file=fname_trace, position='append', action='write', status='unknown') + write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') 'SEND phase=ACCUM axis=${X}$ k=', & + & k, ' from_rank=', proc_rank, ' to_rank=', send_neighbor, ' t_step=', dbg_t_step, & + & ' ft=', send_ft(:,i) + close (unit_trace) + end if + end do + end if 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) @@ -1399,8 +1433,20 @@ contains 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_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques, & - & recv_forces_snap, recv_torques_snap]') + if (trace_active) then + do i = 1, recv_count + if (recv_ids(i) == dbg_track_gbl_id) then + open (unit=unit_trace, file=fname_trace, position='append', action='write', & + & status='unknown') + write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & + & 'RECV phase=ACCUM axis=${X}$ k=', k, ' at_rank=', proc_rank, ' from_rank=', & + & recv_neighbor, ' t_step=', dbg_t_step, ' ft=', recv_ft(:,i) + close (unit_trace) + end if + end do + end if + ! $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques, & + ! & recv_forces_snap, recv_torques_snap]') do i = 1, recv_count call s_get_neighborhood_idx(recv_ids(i), j) if (j > 0) then @@ -1413,7 +1459,21 @@ contains end do end if end do - $:END_GPU_PARALLEL_LOOP() + ! $:END_GPU_PARALLEL_LOOP() + ! TEMPORARY: checkpoint forces/torques for dbg_track_gbl_id (resolved fresh via lookup, not a + ! hardcoded local index) right after this receive loop, independent of any id-matching logic in + ! the loop above, to see if it changed even though nothing should have touched it. + if (trace_active) then + call s_get_neighborhood_idx(dbg_track_gbl_id, j) + if (j > 0) then + open (unit=unit_trace, file=fname_trace, position='append', action='write', status='unknown') + write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & + & 'CHECKPOINT phase=ACCUM after=RECV_LOOP axis=${X}$ k=', k, ' at_rank=', proc_rank, & + & ' num_ibs=', num_ibs, ' local_idx=', j, ' forces+torques=', forces(j,:), torques(j,:) + call s_debug_write_ib_lookup_state(unit_trace) + close (unit_trace) + end if + end if end if tag = tag + 2 end do @@ -1428,7 +1488,7 @@ contains do k = 1, min(2*ib_neighborhood_radius, num_procs_${X}$ - 1) pack_pos = 0 - $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') + ! $: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 @@ -1436,8 +1496,19 @@ contains send_ft(l + 3,i) = torques(i,l) end do end do - $:END_GPU_PARALLEL_LOOP() - $:GPU_UPDATE(host='[send_ids, send_ft]') + ! $:END_GPU_PARALLEL_LOOP() + ! $:GPU_UPDATE(host='[send_ids, send_ft]') + if (trace_active) then + do i = 1, num_ibs + if (send_ids(i) == dbg_track_gbl_id) then + open (unit=unit_trace, file=fname_trace, position='append', action='write', status='unknown') + write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & + & 'SEND phase=BACKPROP axis=${X}$ k=', k, ' from_rank=', proc_rank, ' to_rank=', & + & send_neighbor, ' t_step=', dbg_t_step, ' ft=', send_ft(:,i) + close (unit_trace) + end if + end do + end if 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) @@ -1449,7 +1520,34 @@ contains 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_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques]') + if (trace_active) then + do i = 1, recv_count + if (recv_ids(i) == dbg_track_gbl_id) then + open (unit=unit_trace, file=fname_trace, position='append', action='write', & + & status='unknown') + write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & + & 'RECV phase=BACKPROP axis=${X}$ k=', k, ' at_rank=', proc_rank, ' from_rank=', & + & recv_neighbor, ' t_step=', dbg_t_step, ' ft=', recv_ft(:,i) + close (unit_trace) + end if + end do + ! TEMPORARY: dump every recv_ids(i) -> local index j resolution this receive will act on, to + ! check for an ib_gbl_idx_lookup collision (a different received id resolving to + ! gbl_patch_id=2's slot). + do i = 1, recv_count + call s_get_neighborhood_idx(recv_ids(i), j) + if (j > 0) then + open (unit=unit_trace, file=fname_trace, position='append', action='write', & + & status='unknown') + write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & + & 'RESOLVE phase=BACKPROP axis=${X}$ k=', k, ' at_rank=', proc_rank, & + & ' recv_ids(i)=', recv_ids(i), ' local_j=', j, ' local_gbl_id=', & + & patch_ib(j)%gbl_patch_id, ' i=', i, ' ft=', recv_ft(:,i) + close (unit_trace) + end if + end do + end if + ! $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques]') do i = 1, recv_count call s_get_neighborhood_idx(recv_ids(i), j) if (j > 0) then @@ -1459,7 +1557,23 @@ contains end do end if end do - $:END_GPU_PARALLEL_LOOP() + ! $:END_GPU_PARALLEL_LOOP() + ! TEMPORARY: checkpoint forces/torques for dbg_track_gbl_id (resolved fresh via lookup, not a + ! hardcoded local index) right after this receive loop, independent of any id-matching logic in + ! the loop above, to see if it changed even though nothing should have touched it. + if (trace_active) then + call s_get_neighborhood_idx(dbg_track_gbl_id, j) + if (j > 0) then + open (unit=unit_trace, file=fname_trace, position='append', action='write', & + & status='unknown') + write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & + & 'CHECKPOINT phase=BACKPROP after=RECV_LOOP axis=${X}$ k=', k, ' at_rank=', & + & proc_rank, ' num_ibs=', num_ibs, ' local_idx=', j, ' forces+torques=', & + & forces(j,:), torques(j,:) + call s_debug_write_ib_lookup_state(unit_trace) + close (unit_trace) + end if + end if end if tag = tag + 2 end do @@ -1608,12 +1722,30 @@ contains @:ASSERT(num_ibs <= size(patch_ib), 'patch_ib overflow in neighborhood handoff') patch_ib(num_ibs) = tmp_patch ib_gbl_idx_lookup(tmp_patch%gbl_patch_id) = num_ibs + + ! TEMPORARY DEBUG INSTRUMENTATION: log every fresh broadcast receipt of dbg_track_gbl_id, ungated + ! by trace_active, so we can see the full history of when this rank last got an updated copy and + ! what it contained at that moment - not just near a detected divergence. + if (tmp_patch%gbl_patch_id == dbg_track_gbl_id) then + block + character(len=64) :: dbg_bfname + integer :: dbg_bunit + write (dbg_bfname, '(A,I0,A)') 'ib_broadcast_trace_rank', proc_rank, '.log' + dbg_bunit = 990 + proc_rank + open (unit=dbg_bunit, file=dbg_bfname, position='append', action='write', & + & status='unknown') + write (dbg_bunit, '(A,I0,A,I0,A,I0)') 'NEW_BROADCAST t_step=', dbg_t_step, ' at_rank=', & + & proc_rank, ' from_rank=', recv_neighbor_list(nbr_idx) + call s_debug_write_ib_state(dbg_bunit, tmp_patch) + close (dbg_bunit) + end block + end if end if end do end do deallocate (send_buf, recv_bufs) - $:GPU_UPDATE(device='[patch_ib]') + $:GPU_UPDATE(device='[patch_ib, num_ibs]') call s_update_ib_lookup() end if #endif @@ -1639,9 +1771,12 @@ contains integer :: r2, jr, roster_unpack_pos logical :: found type(ib_patch_parameters) :: roster_patch + logical :: mismatch_found_this_call if (num_procs == 1) return + mismatch_found_this_call = .false. + ! Gather each rank's Cartesian coords + flow-field boundary neighbor ranks (bc_x/y/z%beg/end) so a divergence dump can ! show the real adjacency graph instead of an assumed one. topo_local = -999 @@ -1703,10 +1838,11 @@ contains & any(patch_ib(i)%step_angles /= other_patch%step_angles) if (mismatch) then + mismatch_found_this_call = .true. open (unit=unit_num, file=fname, position='append', action='write', status='unknown') write (unit_num, '(A)') '====================================================================' - write (unit_num, '(A,I0,A,ES23.15,A,I0,A,I0,A,I0)') 'DIVERGENCE t_step=', t_step, ' mytime=', mytime, & - & ' gbl_patch_id=', patch_ib(i)%gbl_patch_id, ' rank_A=', proc_rank, ' rank_B=', r + write (unit_num, '(A,I0,A,I0,A,I0,A,I0)') 'DIVERGENCE t_step=', t_step, ' gbl_patch_id=', & + & patch_ib(i)%gbl_patch_id, ' rank_A=', proc_rank, ' rank_B=', r write (unit_num, '(A,I0,A,I0,A,I0,A,I0,A,I0)') 'num_procs_x=', num_procs_x, ' num_procs_y=', & & num_procs_y, ' num_procs_z=', num_procs_z, ' ib_neighborhood_radius=', ib_neighborhood_radius, & & ' num_dims=', num_dims @@ -1760,6 +1896,16 @@ contains end do deallocate (send_buf, recv_buf) + + ! TEMPORARY DEBUG INSTRUMENTATION: stop after two distinct divergent timesteps have been captured, so the log + ! files don't grow unbounded once the mechanism is confirmed reproducible. + if (mismatch_found_this_call .and. t_step /= dbg_last_divergent_t_step) then + dbg_last_divergent_t_step = t_step + dbg_divergence_count = dbg_divergence_count + 1 + if (dbg_divergence_count >= 2) then + call s_mpi_abort('TEMPORARY DEBUG INSTRUMENTATION: stopping after capturing 2 divergent timesteps') + end if + end if #endif end subroutine s_debug_log_ib_divergence @@ -1785,6 +1931,47 @@ contains end subroutine s_debug_write_ib_state + !> TEMPORARY DEBUG INSTRUMENTATION: dumps local_ib_patch_ids (this rank's locally-owned patches) and every non-negative + !! entry of ib_gbl_idx_lookup (every patch this rank tracks, local+shadow), plus explicit self-consistency and + !! collision checks, so a stale or duplicate lookup entry is directly visible instead of inferred. + subroutine s_debug_write_ib_lookup_state(unit_num) + + integer, intent(in) :: unit_num + integer :: gi, gi2, jloc + + write (unit_num, '(A,I0)') 'local_ib_patch_ids: num_local_ibs=', num_local_ibs + do jloc = 1, num_local_ibs + write (unit_num, '(A,I0,A,I0,A,I0)') ' local_owner_idx=', jloc, ' -> patch_ib_idx=', local_ib_patch_ids(jloc), & + & ' gbl_patch_id=', patch_ib(local_ib_patch_ids(jloc))%gbl_patch_id + end do + + write (unit_num, '(A,I0,A,I0)') 'ib_gbl_idx_lookup non-negative entries: num_gbl_ibs=', num_gbl_ibs, ' num_ibs=', & + & num_ibs + do gi = 1, num_gbl_ibs + if (ib_gbl_idx_lookup(gi) > 0) then + write (unit_num, '(A,I0,A,I0)') ' gbl_id=', gi, ' -> local_idx=', ib_gbl_idx_lookup(gi) + end if + end do + + do jloc = 1, num_ibs + if (ib_gbl_idx_lookup(patch_ib(jloc)%gbl_patch_id) /= jloc) then + write (unit_num, '(A,I0,A,I0,A,I0)') ' INCONSISTENT: patch_ib(', jloc, ')%gbl_patch_id=', & + & patch_ib(jloc)%gbl_patch_id, ' but lookup points to local_idx=', & + & ib_gbl_idx_lookup(patch_ib(jloc)%gbl_patch_id) + end if + end do + do gi = 1, num_gbl_ibs + if (ib_gbl_idx_lookup(gi) <= 0) cycle + do gi2 = gi + 1, num_gbl_ibs + if (ib_gbl_idx_lookup(gi2) == ib_gbl_idx_lookup(gi)) then + write (unit_num, '(A,I0,A,I0,A,I0)') ' COLLISION: gbl_id=', gi, ' and gbl_id=', gi2, & + & ' both map to local_idx=', ib_gbl_idx_lookup(gi) + end if + end do + end do + + end subroutine s_debug_write_ib_lookup_state + subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx) $:GPU_ROUTINE(parallelism='[seq]') @@ -1797,25 +1984,20 @@ contains end subroutine s_get_neighborhood_idx + !> Rebuilds ib_gbl_idx_lookup from patch_ib on the host and mirrors the result to the device. Built on the host because + !! patch_ib is already host-current at every call site (compaction and neighbor-unpack are host-side Fortran), and + !! because the populate step is a scatter write (target index = patch_ib(i)%gbl_patch_id, not the loop index) that was + !! found to leave stale entries behind under GPU offload after compaction shrinks num_ibs. subroutine s_update_ib_lookup() integer :: i - ! reset the lookup - $:GPU_PARALLEL_LOOP(private='[i]') - do i = 1, num_gbl_ibs - ib_gbl_idx_lookup(i) = -1 - end do - $:END_GPU_PARALLEL_LOOP() - - ! populate the table - $:GPU_PARALLEL_LOOP(private='[i]') + ib_gbl_idx_lookup = -1 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]') + $:GPU_UPDATE(device='[ib_gbl_idx_lookup]') end subroutine s_update_ib_lookup diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 5583872a0..7881e193c 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -452,6 +452,8 @@ contains call cpu_time(start) call nvtxStartRange("TIMESTEP") + dbg_t_step = t_step ! TEMPORARY DEBUG INSTRUMENTATION: see m_ibm + ! Adaptive dt: initial stage if (adap_dt) call s_adaptive_dt_bubble(1) From 6780e3a90536d61b9c6d721470b20824f7540981 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Sun, 30 Aug 2026 05:40:08 -0400 Subject: [PATCH 07/10] Fixed out of bounds memory read when integrating forces --- src/common/m_finite_differences.fpp | 12 ++++++++---- src/post_process/m_derived_variables.fpp | 10 ++++++---- src/simulation/m_bubbles_EL.fpp | 7 ++++--- src/simulation/m_derived_variables.fpp | 14 +++++++++----- src/simulation/m_hypoelastic.fpp | 7 ++++--- 5 files changed, 31 insertions(+), 19 deletions(-) diff --git a/src/common/m_finite_differences.fpp b/src/common/m_finite_differences.fpp index 5639890da..63ac43cb6 100644 --- a/src/common/m_finite_differences.fpp +++ b/src/common/m_finite_differences.fpp @@ -30,12 +30,16 @@ contains real(wp), dimension(-local_buff_size:q + local_buff_size), intent(in) :: s_cc integer :: i !< Generic loop iterator + ! Coefficients always extend at least fd_number_in beyond the interior on each side, so a stencil centered on a + ! ghost-adjacent cell (e.g. an immersed boundary near a domain boundary) has a real coefficient to read instead of + ! reading past the caller's allocation. offset_s, when given, widens this further (never narrows it) for callers + ! that need more than fd_number_in of margin. if (present(offset_s)) then - lB = -offset_s%beg - lE = q + offset_s%end + lB = -max(fd_number_in, offset_s%beg) + lE = q + max(fd_number_in, offset_s%end) else - lB = 0 - lE = q + lB = -fd_number_in + lE = q + fd_number_in end if ! Computing the 1st order finite-difference coefficients diff --git a/src/post_process/m_derived_variables.fpp b/src/post_process/m_derived_variables.fpp index 05f1fc180..b1d1c9a04 100644 --- a/src/post_process/m_derived_variables.fpp +++ b/src/post_process/m_derived_variables.fpp @@ -32,18 +32,20 @@ contains allocate (fd%gm_rho_sf(-offset_x%beg:m + offset_x%end,-offset_y%beg:n + offset_y%end,-offset_z%beg:p + offset_z%end)) end if - ! Allocate FD coefficients (up to 4th order; higher orders need extension) + ! Allocate FD coefficients (up to 4th order; higher orders need extension). s_compute_finite_difference_coefficients + ! always extends at least fd_number beyond the interior on each side, widened further by offset_x/y/z when those + ! are larger (multi-block Silo ghost zones); the allocation must cover whichever bound ends up wider. if (omega_wrt(2) .or. omega_wrt(3) .or. qm_wrt .or. schlieren_wrt .or. liutex_wrt) then - allocate (fd%fd_coeff_x(-fd_number:fd_number,-offset_x%beg:m + offset_x%end)) + allocate (fd%fd_coeff_x(-fd_number:fd_number,-max(fd_number, offset_x%beg):m + max(fd_number, offset_x%end))) end if if (omega_wrt(1) .or. omega_wrt(3) .or. qm_wrt .or. liutex_wrt .or. (n > 0 .and. schlieren_wrt)) then - allocate (fd%fd_coeff_y(-fd_number:fd_number,-offset_y%beg:n + offset_y%end)) + allocate (fd%fd_coeff_y(-fd_number:fd_number,-max(fd_number, offset_y%beg):n + max(fd_number, offset_y%end))) end if if (omega_wrt(1) .or. omega_wrt(2) .or. qm_wrt .or. liutex_wrt .or. (p > 0 .and. schlieren_wrt)) then - allocate (fd%fd_coeff_z(-fd_number:fd_number,-offset_z%beg:p + offset_z%end)) + allocate (fd%fd_coeff_z(-fd_number:fd_number,-max(fd_number, offset_z%beg):p + max(fd_number, offset_z%end))) end if end subroutine s_initialize_derived_variables_module diff --git a/src/simulation/m_bubbles_EL.fpp b/src/simulation/m_bubbles_EL.fpp index 62a512ef1..93d285258 100644 --- a/src/simulation/m_bubbles_EL.fpp +++ b/src/simulation/m_bubbles_EL.fpp @@ -185,18 +185,19 @@ contains if (lag_params%vel_model > 0 .and. lag_params%pressure_force) then @:ALLOCATE(grad_p_x(0:m, 0:n, 0:p)) - @:ALLOCATE(fd_coeff_x_pgrad(-fd_number:fd_number, 0:m)) + ! s_compute_finite_difference_coefficients always extends fd_number beyond the interior on each side + @:ALLOCATE(fd_coeff_x_pgrad(-fd_number:fd_number,-fd_number:m + fd_number)) call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x_pgrad, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_x_pgrad]') if (n > 0) then @:ALLOCATE(grad_p_y(0:m, 0:n, 0:p)) - @:ALLOCATE(fd_coeff_y_pgrad(-fd_number:fd_number, 0:n)) + @:ALLOCATE(fd_coeff_y_pgrad(-fd_number:fd_number,-fd_number:n + fd_number)) call s_compute_finite_difference_coefficients(n, y_cc, fd_coeff_y_pgrad, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_y_pgrad]') end if if (p > 0) then @:ALLOCATE(grad_p_z(0:m, 0:n, 0:p)) - @:ALLOCATE(fd_coeff_z_pgrad(-fd_number:fd_number, 0:p)) + @:ALLOCATE(fd_coeff_z_pgrad(-fd_number:fd_number,-fd_number:p + fd_number)) call s_compute_finite_difference_coefficients(p, z_cc, fd_coeff_z_pgrad, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_z_pgrad]') end if diff --git a/src/simulation/m_derived_variables.fpp b/src/simulation/m_derived_variables.fpp index 6a4b92ddc..c822e0962 100644 --- a/src/simulation/m_derived_variables.fpp +++ b/src/simulation/m_derived_variables.fpp @@ -38,14 +38,17 @@ contains ! higher than fourth-order accuracy coefficients are wanted, the formulae required to compute these coefficients will have ! to be implemented in the subroutine s_compute_finite_difference_coefficients. - ! Allocating centered finite-difference coefficients + ! Allocating centered finite-difference coefficients. The coefficient (second) index is extended by fd_number beyond + ! the interior on each side: s_compute_ib_forces evaluates the viscous-stress stencil centered on ghost-adjacent + ! cells (i+l for l in -fd_number:fd_number) when an IB sits near a domain boundary, so the coefficient array must + ! cover those centers too, not just the interior 0:m. if (probe_wrt .or. ib) then - @:ALLOCATE(fd_coeff_x(-fd_number:fd_number, 0:m)) + @:ALLOCATE(fd_coeff_x(-fd_number:fd_number,-fd_number:m + fd_number)) if (n > 0) then - @:ALLOCATE(fd_coeff_y(-fd_number:fd_number, 0:n)) + @:ALLOCATE(fd_coeff_y(-fd_number:fd_number,-fd_number:n + fd_number)) end if if (p > 0) then - @:ALLOCATE(fd_coeff_z(-fd_number:fd_number, 0:p)) + @:ALLOCATE(fd_coeff_z(-fd_number:fd_number,-fd_number:p + fd_number)) end if @:ALLOCATE(accel_mag(0:m, 0:n, 0:p)) @@ -69,7 +72,8 @@ contains call s_open_probe_files() call s_open_com_files() end if - ! Computing centered finite difference coefficients + ! Computing centered finite difference coefficients (s_compute_finite_difference_coefficients always extends + ! fd_number beyond the interior on each side; the allocation above matches) call s_compute_finite_difference_coefficients(m, x_cc, fd_coeff_x, buff_size, fd_number, fd_order) $:GPU_UPDATE(device='[fd_coeff_x]') diff --git a/src/simulation/m_hypoelastic.fpp b/src/simulation/m_hypoelastic.fpp index 3dce0a5ae..1b9560a03 100644 --- a/src/simulation/m_hypoelastic.fpp +++ b/src/simulation/m_hypoelastic.fpp @@ -57,12 +57,13 @@ contains end do $:GPU_UPDATE(device='[Gs_hypo]') - @:ALLOCATE(fd_coeff_x_hypo(-fd_number:fd_number, 0:m)) + ! s_compute_finite_difference_coefficients always extends fd_number beyond the interior on each side + @:ALLOCATE(fd_coeff_x_hypo(-fd_number:fd_number,-fd_number:m + fd_number)) if (n > 0) then - @:ALLOCATE(fd_coeff_y_hypo(-fd_number:fd_number, 0:n)) + @:ALLOCATE(fd_coeff_y_hypo(-fd_number:fd_number,-fd_number:n + fd_number)) end if if (p > 0) then - @:ALLOCATE(fd_coeff_z_hypo(-fd_number:fd_number, 0:p)) + @:ALLOCATE(fd_coeff_z_hypo(-fd_number:fd_number,-fd_number:p + fd_number)) end if ! Computing centered finite difference coefficients From 2fe42f0d4066258f0ffb5f05c8638f5a1cc1a665 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Sun, 30 Aug 2026 07:08:46 -0400 Subject: [PATCH 08/10] Found identical GPU update issue with num of ghost points --- src/simulation/m_data_output.fpp | 32 +++++++++++++++++++--- src/simulation/m_ibm.fpp | 47 +++++++++++++++++++++++++++++--- 2 files changed, 71 insertions(+), 8 deletions(-) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 632d248ec..0319802b8 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -396,12 +396,36 @@ contains near2_dist = dist; near2_id = i end if end do - if (near1_id > 0) print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' nearest particle id=', near1_id, ' dist=', & - & near1_dist, ' gap=', near1_dist - patch_ib(near1_id)%radius, ' vel=', patch_ib(near1_id)%vel - if (near2_id > 0) print '(A,I0,A,ES16.6,A,ES16.6)', ' 2nd nearest particle id=', near2_id, ' dist=', near2_dist, & - & ' gap=', near2_dist - patch_ib(near2_id)%radius + if (near1_id > 0) then + print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' nearest particle id=', near1_id, ' dist=', & + & near1_dist, ' gap=', near1_dist - patch_ib(near1_id)%radius, ' vel=', patch_ib(near1_id)%vel + print '(A,3(ES16.6,1X))', ' centroid = ', patch_ib(near1_id)%x_centroid, patch_ib(near1_id)%y_centroid, & + & patch_ib(near1_id)%z_centroid + print '(A,3(ES16.6,1X))', ' angular_vel = ', patch_ib(near1_id)%angular_vel + print '(A,3(ES16.6,1X))', ' force = ', patch_ib(near1_id)%force + print '(A,3(ES16.6,1X))', ' torque = ', patch_ib(near1_id)%torque + print '(A,I0,A,ES16.6,A,ES16.6)', ' moving_ibm = ', patch_ib(near1_id)%moving_ibm, ' mass=', & + & patch_ib(near1_id)%mass, ' moment=', patch_ib(near1_id)%moment + end if + if (near2_id > 0) then + print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' 2nd nearest particle id=', near2_id, ' dist=', & + & near2_dist, ' gap=', near2_dist - patch_ib(near2_id)%radius, ' vel=', patch_ib(near2_id)%vel + print '(A,3(ES16.6,1X))', ' centroid = ', patch_ib(near2_id)%x_centroid, patch_ib(near2_id)%y_centroid, & + & patch_ib(near2_id)%z_centroid + print '(A,3(ES16.6,1X))', ' angular_vel = ', patch_ib(near2_id)%angular_vel + print '(A,3(ES16.6,1X))', ' force = ', patch_ib(near2_id)%force + end if end if + ! TEMPORARY DEBUG INSTRUMENTATION: dump a small x-neighborhood around the violating cell (including into the ghost/ + ! halo region on either side) to check for a sharp discontinuity right at a processor boundary versus a smoothly + ! diverging field, since ICFL blowups have been observed specifically near rank boundaries. + print '(A)', ' x-neighborhood (dj, rho, pres, vel) around violating cell:' + do j = max(-buff_size, j_hit - 3), min(m + buff_size, j_hit + 3) + call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k_hit, l_hit) + print '(A,I0,A,ES16.6,A,ES16.6,A,3(ES16.6,1X))', ' dj=', j - j_hit, ' rho=', rho, ' pres=', pres, ' vel=', vel + end do + call flush (6) end subroutine s_report_icfl_violation diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 6ee16f8f3..8916160ae 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -193,6 +193,14 @@ contains type(ghost_point) :: gp type(ghost_point) :: innerp + ! TEMPORARY DEBUG INSTRUMENTATION: test whether clamping the moving-IB pressure-correction denominator (below) away + ! from zero prevents the ICFL blowups observed near strongly-coupled, high-force IBs. dbg_denom is the per-thread + ! (unclamped) denominator value; dbg_clamp_count/dbg_min_abs_denom are reduced across all ghost points this call. + real(wp) :: dbg_denom + integer :: dbg_clamp_count + real(wp) :: dbg_min_abs_denom + real(wp), parameter :: dbg_denom_floor = 0.05_wp + ! Per-ghost-point image-point interpolation results, stashed between the interpolation ! kernel and the correction kernel below so that the correction kernel never reads ! q_prim_vf (or pb_in/mv_in) at a cell another ghost point's correction may have already @@ -240,6 +248,9 @@ contains $:END_GPU_PARALLEL_LOOP() if (num_gps > 0) then + dbg_clamp_count = 0 + dbg_min_abs_denom = huge(1._wp) + @:ALLOCATE(alpha_rho_IP_buf(1:num_fluids, 1:num_gps), alpha_IP_buf(1:num_fluids, 1:num_gps), & & pres_IP_buf(1:num_gps), c_IP_buf(1:num_gps), vel_IP_buf(1:3, 1:num_gps), & & r_IP_buf(1:nb, 1:num_gps), v_IP_buf(1:nb, 1:num_gps), pb_IP_buf(1:nb, 1:num_gps), & @@ -293,7 +304,7 @@ contains $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, & & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & & innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id, & - & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff]') + & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, dbg_denom]', copy='[dbg_clamp_count, dbg_min_abs_denom]') do i = 1, num_gps gp = ghost_points(i) j = gp%loc(1) @@ -357,9 +368,20 @@ contains $:GPU_LOOP(parallelism='[seq]') do q = 1, num_fluids ! Pressure correction for moving IB: accounts for acceleration of IB surface - q_prim_vf(eqn_idx%E)%sf(j, k, l) = q_prim_vf(eqn_idx%E)%sf(j, k, & - & l) + pres_IP/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & - & *dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, gp%levelset_norm)) + dbg_denom = 1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & + & *dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, gp%levelset_norm) + + ! TEMPORARY DEBUG INSTRUMENTATION: test whether flooring |dbg_denom| away from zero (preserving its + ! sign) prevents the observed ICFL blowups near strongly-coupled, high-force IBs. + $:GPU_ATOMIC(atomic='update') + dbg_min_abs_denom = min(dbg_min_abs_denom, abs(dbg_denom)) + if (abs(dbg_denom) < dbg_denom_floor) then + $:GPU_ATOMIC(atomic='update') + dbg_clamp_count = dbg_clamp_count + 1 + dbg_denom = sign(dbg_denom_floor, dbg_denom) + end if + + q_prim_vf(eqn_idx%E)%sf(j, k, l) = q_prim_vf(eqn_idx%E)%sf(j, k, l) + pres_IP/dbg_denom end do end if @@ -520,6 +542,22 @@ contains end do $:END_GPU_PARALLEL_LOOP() + ! TEMPORARY DEBUG INSTRUMENTATION: report whenever the pressure-correction denominator clamp above actually + ! engaged this call, so we can confirm whether it's firing (and how close to zero it was) before/around an ICFL + ! blowup. + if (dbg_clamp_count > 0) then + block + character(len=64) :: dbg_fname + integer :: dbg_unit + write (dbg_fname, '(A,I0,A)') 'ib_pres_clamp_rank', proc_rank, '.log' + dbg_unit = 980 + proc_rank + open (unit=dbg_unit, file=dbg_fname, position='append', action='write', status='unknown') + write (dbg_unit, '(A,I0,A,I0,A,ES16.6)') 't_step=', dbg_t_step, ' clamp_count=', dbg_clamp_count, & + & ' min_abs_denom=', dbg_min_abs_denom + close (dbg_unit) + end block + end if + @:DEALLOCATE(alpha_rho_IP_buf, alpha_IP_buf, pres_IP_buf, c_IP_buf, vel_IP_buf, r_IP_buf, v_IP_buf, pb_IP_buf, & & mv_IP_buf, nmom_IP_buf, presb_IP_buf, massv_IP_buf, Ys_IP_buf) end if @@ -1030,6 +1068,7 @@ contains call nvtxStartRange("COMPUTE-GHOST-POINTS") ! recalculate the ghost point locations and coefficients call s_find_num_ghost_points(num_gps) + $:GPU_UPDATE(device='[num_gps]') call s_find_ghost_points(ghost_points) call nvtxEndRange From a4f3fab1a5efa98dd13f4fbb2cf9039fcf73de3e Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Sun, 30 Aug 2026 07:33:30 -0400 Subject: [PATCH 09/10] Cleaning up and running final tests. Things are looking very stable --- src/simulation/m_data_output.fpp | 6 +- src/simulation/m_ibm.fpp | 420 ++--------------------------- src/simulation/m_time_steppers.fpp | 3 - 3 files changed, 18 insertions(+), 411 deletions(-) diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp index 0319802b8..f013dd92a 100644 --- a/src/simulation/m_data_output.fpp +++ b/src/simulation/m_data_output.fpp @@ -417,9 +417,9 @@ contains end if end if - ! TEMPORARY DEBUG INSTRUMENTATION: dump a small x-neighborhood around the violating cell (including into the ghost/ - ! halo region on either side) to check for a sharp discontinuity right at a processor boundary versus a smoothly - ! diverging field, since ICFL blowups have been observed specifically near rank boundaries. + ! Dump a small x-neighborhood around the violating cell (reaching into the ghost/halo region on either side) to + ! distinguish a sharp discontinuity at a processor boundary - the signature of stale or corrupted halo/IB state - + ! from a smoothly diverging field, which indicates a genuine physical/numerical instability. print '(A)', ' x-neighborhood (dj, rho, pres, vel) around violating cell:' do j = max(-buff_size, j_hit - 3), min(m + buff_size, j_hit + 3) call s_compute_enthalpy(q_prim_vf, pres, rho, gamma, pi_inf, Re, H, alpha, vel, vel_sum, qv, j, k_hit, l_hit) diff --git a/src/simulation/m_ibm.fpp b/src/simulation/m_ibm.fpp index 8916160ae..32fced2bc 100644 --- a/src/simulation/m_ibm.fpp +++ b/src/simulation/m_ibm.fpp @@ -48,15 +48,6 @@ module m_ibm real(wp), allocatable :: send_ft(:,:), recv_ft(:,:) real(wp), allocatable :: recv_forces_snap(:,:), recv_torques_snap(:,:) - ! TEMPORARY DEBUG INSTRUMENTATION state (remove once the IB force-communication bug is found): - ! dbg_t_step is set once per timestep by m_time_steppers; dbg_divergence_count/dbg_last_divergent_t_step track how many - ! distinct timesteps s_debug_log_ib_divergence has found a mismatch on, to gate the wire trace and stop the run once - ! two measurements are captured. - integer, save :: dbg_t_step = -1 - integer, save :: dbg_divergence_count = 0 - integer, save :: dbg_last_divergent_t_step = -1 - integer, parameter :: dbg_track_gbl_id = 352 ! gbl_patch_id currently under investigation - contains !> Allocates memory for the variables in the IBM module @@ -193,14 +184,6 @@ contains type(ghost_point) :: gp type(ghost_point) :: innerp - ! TEMPORARY DEBUG INSTRUMENTATION: test whether clamping the moving-IB pressure-correction denominator (below) away - ! from zero prevents the ICFL blowups observed near strongly-coupled, high-force IBs. dbg_denom is the per-thread - ! (unclamped) denominator value; dbg_clamp_count/dbg_min_abs_denom are reduced across all ghost points this call. - real(wp) :: dbg_denom - integer :: dbg_clamp_count - real(wp) :: dbg_min_abs_denom - real(wp), parameter :: dbg_denom_floor = 0.05_wp - ! Per-ghost-point image-point interpolation results, stashed between the interpolation ! kernel and the correction kernel below so that the correction kernel never reads ! q_prim_vf (or pb_in/mv_in) at a cell another ghost point's correction may have already @@ -248,9 +231,6 @@ contains $:END_GPU_PARALLEL_LOOP() if (num_gps > 0) then - dbg_clamp_count = 0 - dbg_min_abs_denom = huge(1._wp) - @:ALLOCATE(alpha_rho_IP_buf(1:num_fluids, 1:num_gps), alpha_IP_buf(1:num_fluids, 1:num_gps), & & pres_IP_buf(1:num_gps), c_IP_buf(1:num_gps), vel_IP_buf(1:3, 1:num_gps), & & r_IP_buf(1:nb, 1:num_gps), v_IP_buf(1:nb, 1:num_gps), pb_IP_buf(1:nb, 1:num_gps), & @@ -304,7 +284,7 @@ contains $:GPU_PARALLEL_LOOP(private='[i, physical_loc, dyn_pres, alpha_rho_IP, alpha_IP, pres_IP, vel_IP, vel_g, vel_norm_IP, & & r_IP, v_IP, pb_IP, mv_IP, nmom_IP, presb_IP, massv_IP, rho, gamma, pi_inf, Re_K, G_K, Gs, gp, & & innerp, norm, buf, radial_vector, rotation_velocity, j, k, l, q, qv_K, c_IP, nbub, patch_id, & - & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff, dbg_denom]', copy='[dbg_clamp_count, dbg_min_abs_denom]') + & Ys_IP, T_IP, mw_IP, e_IP, v_blow_eff]') do i = 1, num_gps gp = ghost_points(i) j = gp%loc(1) @@ -368,20 +348,9 @@ contains $:GPU_LOOP(parallelism='[seq]') do q = 1, num_fluids ! Pressure correction for moving IB: accounts for acceleration of IB surface - dbg_denom = 1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & - & *dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, gp%levelset_norm) - - ! TEMPORARY DEBUG INSTRUMENTATION: test whether flooring |dbg_denom| away from zero (preserving its - ! sign) prevents the observed ICFL blowups near strongly-coupled, high-force IBs. - $:GPU_ATOMIC(atomic='update') - dbg_min_abs_denom = min(dbg_min_abs_denom, abs(dbg_denom)) - if (abs(dbg_denom) < dbg_denom_floor) then - $:GPU_ATOMIC(atomic='update') - dbg_clamp_count = dbg_clamp_count + 1 - dbg_denom = sign(dbg_denom_floor, dbg_denom) - end if - - q_prim_vf(eqn_idx%E)%sf(j, k, l) = q_prim_vf(eqn_idx%E)%sf(j, k, l) + pres_IP/dbg_denom + q_prim_vf(eqn_idx%E)%sf(j, k, l) = q_prim_vf(eqn_idx%E)%sf(j, k, & + & l) + pres_IP/(1._wp - 2._wp*abs(gp%levelset*alpha_rho_IP(q)/pres_IP) & + & *dot_product(patch_ib(patch_id)%force/patch_ib(patch_id)%mass, gp%levelset_norm)) end do end if @@ -542,22 +511,6 @@ contains end do $:END_GPU_PARALLEL_LOOP() - ! TEMPORARY DEBUG INSTRUMENTATION: report whenever the pressure-correction denominator clamp above actually - ! engaged this call, so we can confirm whether it's firing (and how close to zero it was) before/around an ICFL - ! blowup. - if (dbg_clamp_count > 0) then - block - character(len=64) :: dbg_fname - integer :: dbg_unit - write (dbg_fname, '(A,I0,A)') 'ib_pres_clamp_rank', proc_rank, '.log' - dbg_unit = 980 + proc_rank - open (unit=dbg_unit, file=dbg_fname, position='append', action='write', status='unknown') - write (dbg_unit, '(A,I0,A,I0,A,ES16.6)') 't_step=', dbg_t_step, ' clamp_count=', dbg_clamp_count, & - & ' min_abs_denom=', dbg_min_abs_denom - close (dbg_unit) - end block - end if - @:DEALLOCATE(alpha_rho_IP_buf, alpha_IP_buf, pres_IP_buf, c_IP_buf, vel_IP_buf, r_IP_buf, v_IP_buf, pb_IP_buf, & & mv_IP_buf, nmom_IP_buf, presb_IP_buf, massv_IP_buf, Ys_IP_buf) end if @@ -1406,20 +1359,9 @@ contains 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(:) - character(len=64) :: fname_trace - integer :: unit_trace - logical :: trace_active if (num_procs == 1) return - ! TEMPORARY DEBUG INSTRUMENTATION (remove once found): trace the wire values for gbl_patch_id=2 at every send/recv, - ! gated to the known failure window only - this does file open/close per match and is too slow to run unconditionally. - write (fname_trace, '(A,I0,A)') 'ib_wire_trace_rank', proc_rank, '.log' - unit_trace = 950 + proc_rank - ! Only trace once the roster check (s_debug_log_ib_divergence) has confirmed at least one divergent timestep, so we - ! don't pay the per-send/recv file I/O cost before we know we're near a real event. - trace_active = (dbg_divergence_count >= 1) - 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)) @@ -1436,10 +1378,7 @@ contains 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 - ! TEMPORARY: GPU offload disabled on this whole subroutine to test whether it's the source of the y/z - ! force-communication corruption bug (see ib_wire_trace evidence) - these loops are tiny (num_ibs - ! entries), no perf concern. - ! $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') + $: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 @@ -1447,19 +1386,8 @@ contains send_ft(l + 3,i) = torques(i,l) end do end do - ! $:END_GPU_PARALLEL_LOOP() - ! $:GPU_UPDATE(host='[send_ids, send_ft]') - if (trace_active) then - do i = 1, num_ibs - if (send_ids(i) == dbg_track_gbl_id) then - open (unit=unit_trace, file=fname_trace, position='append', action='write', status='unknown') - write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') 'SEND phase=ACCUM axis=${X}$ k=', & - & k, ' from_rank=', proc_rank, ' to_rank=', send_neighbor, ' t_step=', dbg_t_step, & - & ' ft=', send_ft(:,i) - close (unit_trace) - end if - end do - end if + $:END_GPU_PARALLEL_LOOP() + $: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) @@ -1472,20 +1400,8 @@ contains 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) - if (trace_active) then - do i = 1, recv_count - if (recv_ids(i) == dbg_track_gbl_id) then - open (unit=unit_trace, file=fname_trace, position='append', action='write', & - & status='unknown') - write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & - & 'RECV phase=ACCUM axis=${X}$ k=', k, ' at_rank=', proc_rank, ' from_rank=', & - & recv_neighbor, ' t_step=', dbg_t_step, ' ft=', recv_ft(:,i) - close (unit_trace) - end if - end do - end if - ! $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques, & - ! & recv_forces_snap, recv_torques_snap]') + $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques, & + & recv_forces_snap, recv_torques_snap]') do i = 1, recv_count call s_get_neighborhood_idx(recv_ids(i), j) if (j > 0) then @@ -1498,21 +1414,7 @@ contains end do end if end do - ! $:END_GPU_PARALLEL_LOOP() - ! TEMPORARY: checkpoint forces/torques for dbg_track_gbl_id (resolved fresh via lookup, not a - ! hardcoded local index) right after this receive loop, independent of any id-matching logic in - ! the loop above, to see if it changed even though nothing should have touched it. - if (trace_active) then - call s_get_neighborhood_idx(dbg_track_gbl_id, j) - if (j > 0) then - open (unit=unit_trace, file=fname_trace, position='append', action='write', status='unknown') - write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & - & 'CHECKPOINT phase=ACCUM after=RECV_LOOP axis=${X}$ k=', k, ' at_rank=', proc_rank, & - & ' num_ibs=', num_ibs, ' local_idx=', j, ' forces+torques=', forces(j,:), torques(j,:) - call s_debug_write_ib_lookup_state(unit_trace) - close (unit_trace) - end if - end if + $:END_GPU_PARALLEL_LOOP() end if tag = tag + 2 end do @@ -1527,7 +1429,7 @@ contains do k = 1, min(2*ib_neighborhood_radius, num_procs_${X}$ - 1) pack_pos = 0 - ! $:GPU_PARALLEL_LOOP(private='[i, l]', copyin='[forces, torques]') + $: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 @@ -1535,19 +1437,8 @@ contains send_ft(l + 3,i) = torques(i,l) end do end do - ! $:END_GPU_PARALLEL_LOOP() - ! $:GPU_UPDATE(host='[send_ids, send_ft]') - if (trace_active) then - do i = 1, num_ibs - if (send_ids(i) == dbg_track_gbl_id) then - open (unit=unit_trace, file=fname_trace, position='append', action='write', status='unknown') - write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & - & 'SEND phase=BACKPROP axis=${X}$ k=', k, ' from_rank=', proc_rank, ' to_rank=', & - & send_neighbor, ' t_step=', dbg_t_step, ' ft=', send_ft(:,i) - close (unit_trace) - end if - end do - end if + $:END_GPU_PARALLEL_LOOP() + $: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) @@ -1559,34 +1450,7 @@ contains 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) - if (trace_active) then - do i = 1, recv_count - if (recv_ids(i) == dbg_track_gbl_id) then - open (unit=unit_trace, file=fname_trace, position='append', action='write', & - & status='unknown') - write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & - & 'RECV phase=BACKPROP axis=${X}$ k=', k, ' at_rank=', proc_rank, ' from_rank=', & - & recv_neighbor, ' t_step=', dbg_t_step, ' ft=', recv_ft(:,i) - close (unit_trace) - end if - end do - ! TEMPORARY: dump every recv_ids(i) -> local index j resolution this receive will act on, to - ! check for an ib_gbl_idx_lookup collision (a different received id resolving to - ! gbl_patch_id=2's slot). - do i = 1, recv_count - call s_get_neighborhood_idx(recv_ids(i), j) - if (j > 0) then - open (unit=unit_trace, file=fname_trace, position='append', action='write', & - & status='unknown') - write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & - & 'RESOLVE phase=BACKPROP axis=${X}$ k=', k, ' at_rank=', proc_rank, & - & ' recv_ids(i)=', recv_ids(i), ' local_j=', j, ' local_gbl_id=', & - & patch_ib(j)%gbl_patch_id, ' i=', i, ' ft=', recv_ft(:,i) - close (unit_trace) - end if - end do - end if - ! $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques]') + $:GPU_PARALLEL_LOOP(private='[i, j, l]', copyin='[recv_ft, recv_ids]', copy='[forces, torques]') do i = 1, recv_count call s_get_neighborhood_idx(recv_ids(i), j) if (j > 0) then @@ -1596,23 +1460,7 @@ contains end do end if end do - ! $:END_GPU_PARALLEL_LOOP() - ! TEMPORARY: checkpoint forces/torques for dbg_track_gbl_id (resolved fresh via lookup, not a - ! hardcoded local index) right after this receive loop, independent of any id-matching logic in - ! the loop above, to see if it changed even though nothing should have touched it. - if (trace_active) then - call s_get_neighborhood_idx(dbg_track_gbl_id, j) - if (j > 0) then - open (unit=unit_trace, file=fname_trace, position='append', action='write', & - & status='unknown') - write (unit_trace, '(A,I0,A,I0,A,I0,A,I0,A,6ES16.8)') & - & 'CHECKPOINT phase=BACKPROP after=RECV_LOOP axis=${X}$ k=', k, ' at_rank=', & - & proc_rank, ' num_ibs=', num_ibs, ' local_idx=', j, ' forces+torques=', & - & forces(j,:), torques(j,:) - call s_debug_write_ib_lookup_state(unit_trace) - close (unit_trace) - end if - end if + $:END_GPU_PARALLEL_LOOP() end if tag = tag + 2 end do @@ -1761,24 +1609,6 @@ contains @:ASSERT(num_ibs <= size(patch_ib), 'patch_ib overflow in neighborhood handoff') patch_ib(num_ibs) = tmp_patch ib_gbl_idx_lookup(tmp_patch%gbl_patch_id) = num_ibs - - ! TEMPORARY DEBUG INSTRUMENTATION: log every fresh broadcast receipt of dbg_track_gbl_id, ungated - ! by trace_active, so we can see the full history of when this rank last got an updated copy and - ! what it contained at that moment - not just near a detected divergence. - if (tmp_patch%gbl_patch_id == dbg_track_gbl_id) then - block - character(len=64) :: dbg_bfname - integer :: dbg_bunit - write (dbg_bfname, '(A,I0,A)') 'ib_broadcast_trace_rank', proc_rank, '.log' - dbg_bunit = 990 + proc_rank - open (unit=dbg_bunit, file=dbg_bfname, position='append', action='write', & - & status='unknown') - write (dbg_bunit, '(A,I0,A,I0,A,I0)') 'NEW_BROADCAST t_step=', dbg_t_step, ' at_rank=', & - & proc_rank, ' from_rank=', recv_neighbor_list(nbr_idx) - call s_debug_write_ib_state(dbg_bunit, tmp_patch) - close (dbg_bunit) - end block - end if end if end do end do @@ -1791,226 +1621,6 @@ contains end subroutine s_handoff_ib_ownership - !> TEMPORARY DEBUG INSTRUMENTATION (remove once the cross-rank IB divergence bug is found). Gathers every rank's tracked - !! patch_ib states and, for every gbl_patch_id two or more ranks both track, logs a full side-by-side dump to - !! ib_divergence_rank.log whenever their dynamic (kinematic/force) fields disagree at all. - subroutine s_debug_log_ib_divergence(t_step) - - integer, intent(in) :: t_step - -#ifdef MFC_MPI - integer :: ierr, patch_bytes, i, j, r, unpack_pos, unit_num - integer, dimension(0:num_procs - 1) :: rank_counts, rank_counts_bytes, rank_displs_bytes - character(len=1), allocatable :: send_buf(:), recv_buf(:) - type(ib_patch_parameters) :: other_patch - character(len=64) :: fname - logical :: mismatch - integer, dimension(9) :: topo_local - integer, dimension(9, 0:num_procs - 1) :: topo_all - integer :: r2, jr, roster_unpack_pos - logical :: found - type(ib_patch_parameters) :: roster_patch - logical :: mismatch_found_this_call - - if (num_procs == 1) return - - mismatch_found_this_call = .false. - - ! Gather each rank's Cartesian coords + flow-field boundary neighbor ranks (bc_x/y/z%beg/end) so a divergence dump can - ! show the real adjacency graph instead of an assumed one. - topo_local = -999 - topo_local(1:num_dims) = proc_coords(1:num_dims) - topo_local(4) = bc_x%beg; topo_local(5) = bc_x%end - if (num_dims >= 2) then - topo_local(6) = bc_y%beg; topo_local(7) = bc_y%end - end if - if (num_dims >= 3) then - topo_local(8) = bc_z%beg; topo_local(9) = bc_z%end - end if - call MPI_ALLGATHER(topo_local, 9, MPI_INTEGER, topo_all, 9, MPI_INTEGER, MPI_COMM_WORLD, ierr) - - call MPI_ALLGATHER(num_ibs, 1, MPI_INTEGER, rank_counts, 1, MPI_INTEGER, MPI_COMM_WORLD, ierr) - - patch_bytes = storage_size(patch_ib(1))/8 - rank_counts_bytes = rank_counts*patch_bytes - rank_displs_bytes(0) = 0 - do r = 1, num_procs - 1 - rank_displs_bytes(r) = rank_displs_bytes(r - 1) + rank_counts_bytes(r - 1) - end do - - allocate (send_buf(max(1, num_ibs*patch_bytes))) - allocate (recv_buf(max(1, sum(rank_counts_bytes)))) - - unpack_pos = 0 - do i = 1, num_ibs - call MPI_PACK(patch_ib(i), patch_bytes, MPI_BYTE, send_buf, size(send_buf), unpack_pos, MPI_COMM_WORLD, ierr) - end do - - call MPI_ALLGATHERV(send_buf, num_ibs*patch_bytes, MPI_PACKED, recv_buf, rank_counts_bytes, rank_displs_bytes, & - & MPI_PACKED, MPI_COMM_WORLD, ierr) - - write (fname, '(A,I0,A)') 'ib_divergence_rank', proc_rank, '.log' - unit_num = 900 + proc_rank - - do r = 0, num_procs - 1 - if (r == proc_rank) cycle - unpack_pos = rank_displs_bytes(r) - do j = 1, rank_counts(r) - call MPI_UNPACK(recv_buf, size(recv_buf), unpack_pos, other_patch, patch_bytes, MPI_BYTE, MPI_COMM_WORLD, ierr) - - do i = 1, num_ibs - if (patch_ib(i)%gbl_patch_id /= other_patch%gbl_patch_id) cycle - - mismatch = (patch_ib(i)%x_centroid /= other_patch%x_centroid) .or. & - & (patch_ib(i)%y_centroid /= other_patch%y_centroid) .or. & - & (patch_ib(i)%z_centroid /= other_patch%z_centroid) .or. & - & any(patch_ib(i)%vel /= other_patch%vel) .or. & - & any(patch_ib(i)%angular_vel /= other_patch%angular_vel) .or. & - & any(patch_ib(i)%angles /= other_patch%angles) .or. & - & any(patch_ib(i)%force /= other_patch%force) .or. & - & any(patch_ib(i)%torque /= other_patch%torque) .or. & - & (patch_ib(i)%step_x_centroid /= other_patch%step_x_centroid) .or. & - & (patch_ib(i)%step_y_centroid /= other_patch%step_y_centroid) .or. & - & (patch_ib(i)%step_z_centroid /= other_patch%step_z_centroid) .or. & - & any(patch_ib(i)%step_vel /= other_patch%step_vel) .or. & - & any(patch_ib(i)%step_angular_vel /= other_patch%step_angular_vel) .or. & - & any(patch_ib(i)%step_angles /= other_patch%step_angles) - - if (mismatch) then - mismatch_found_this_call = .true. - open (unit=unit_num, file=fname, position='append', action='write', status='unknown') - write (unit_num, '(A)') '====================================================================' - write (unit_num, '(A,I0,A,I0,A,I0,A,I0)') 'DIVERGENCE t_step=', t_step, ' gbl_patch_id=', & - & patch_ib(i)%gbl_patch_id, ' rank_A=', proc_rank, ' rank_B=', r - write (unit_num, '(A,I0,A,I0,A,I0,A,I0,A,I0)') 'num_procs_x=', num_procs_x, ' num_procs_y=', & - & num_procs_y, ' num_procs_z=', num_procs_z, ' ib_neighborhood_radius=', ib_neighborhood_radius, & - & ' num_dims=', num_dims - write (unit_num, '(A,I0,A,3I3,A,4I5,A,2I5)') '--- rank ', proc_rank, & - & ' coords=', topo_all(1:3, proc_rank), ' bc_x(beg,end)/bc_y(beg,end)=', & - & topo_all(4:5, proc_rank), topo_all(6:7, proc_rank), ' bc_z(beg,end)=', topo_all(8:9, proc_rank) - call s_debug_write_ib_state(unit_num, patch_ib(i)) - write (unit_num, '(A,I0,A,3I3,A,4I5,A,2I5)') '--- rank ', r, & - & ' coords=', topo_all(1:3, r), ' bc_x(beg,end)/bc_y(beg,end)=', & - & topo_all(4:5, r), topo_all(6:7, r), ' bc_z(beg,end)=', topo_all(8:9, r) - call s_debug_write_ib_state(unit_num, other_patch) - - ! Full roster: every rank's tracking status for this gbl_patch_id, not just the mismatching pair, so a - ! rank silently sending 0s (or not tracking at all) shows up instead of being inferred from absence. - write (unit_num, '(A)') '--- full roster for this gbl_patch_id ---' - do r2 = 0, num_procs - 1 - if (r2 == proc_rank) then - found = .false. - do jr = 1, num_ibs - if (patch_ib(jr)%gbl_patch_id == patch_ib(i)%gbl_patch_id) then - found = .true. - write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, ' TRACKS coords=', topo_all(1:3, r2) - call s_debug_write_ib_state(unit_num, patch_ib(jr)) - exit - end if - end do - if (.not. found) write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, & - & ' NOT TRACKED coords=', topo_all(1:3, r2) - else - found = .false. - roster_unpack_pos = rank_displs_bytes(r2) - do jr = 1, rank_counts(r2) - call MPI_UNPACK(recv_buf, size(recv_buf), roster_unpack_pos, roster_patch, patch_bytes, & - & MPI_BYTE, MPI_COMM_WORLD, ierr) - if (roster_patch%gbl_patch_id == patch_ib(i)%gbl_patch_id) then - found = .true. - write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, ' TRACKS coords=', topo_all(1:3, r2) - call s_debug_write_ib_state(unit_num, roster_patch) - exit - end if - end do - if (.not. found) write (unit_num, '(A,I0,A,3I3)') 'rank ', r2, & - & ' NOT TRACKED coords=', topo_all(1:3, r2) - end if - end do - - close (unit_num) - end if - end do - end do - end do - - deallocate (send_buf, recv_buf) - - ! TEMPORARY DEBUG INSTRUMENTATION: stop after two distinct divergent timesteps have been captured, so the log - ! files don't grow unbounded once the mechanism is confirmed reproducible. - if (mismatch_found_this_call .and. t_step /= dbg_last_divergent_t_step) then - dbg_last_divergent_t_step = t_step - dbg_divergence_count = dbg_divergence_count + 1 - if (dbg_divergence_count >= 2) then - call s_mpi_abort('TEMPORARY DEBUG INSTRUMENTATION: stopping after capturing 2 divergent timesteps') - end if - end if -#endif - - end subroutine s_debug_log_ib_divergence - - !> TEMPORARY DEBUG INSTRUMENTATION helper for s_debug_log_ib_divergence: dumps every dynamic field of an IB patch state. - subroutine s_debug_write_ib_state(unit_num, patch) - - integer, intent(in) :: unit_num - type(ib_patch_parameters), intent(in) :: patch - - write (unit_num, '(A,3ES23.15)') ' centroid = ', patch%x_centroid, patch%y_centroid, patch%z_centroid - write (unit_num, '(A,3ES23.15)') ' step_centroid = ', patch%step_x_centroid, patch%step_y_centroid, & - & patch%step_z_centroid - write (unit_num, '(A,3ES23.15)') ' vel = ', patch%vel - write (unit_num, '(A,3ES23.15)') ' step_vel = ', patch%step_vel - write (unit_num, '(A,3ES23.15)') ' angular_vel = ', patch%angular_vel - write (unit_num, '(A,3ES23.15)') ' step_angular_vel = ', patch%step_angular_vel - write (unit_num, '(A,3ES23.15)') ' angles = ', patch%angles - write (unit_num, '(A,3ES23.15)') ' step_angles = ', patch%step_angles - write (unit_num, '(A,3ES23.15)') ' force = ', patch%force - write (unit_num, '(A,3ES23.15)') ' torque = ', patch%torque - write (unit_num, '(A,ES23.15)') ' moment = ', patch%moment - - end subroutine s_debug_write_ib_state - - !> TEMPORARY DEBUG INSTRUMENTATION: dumps local_ib_patch_ids (this rank's locally-owned patches) and every non-negative - !! entry of ib_gbl_idx_lookup (every patch this rank tracks, local+shadow), plus explicit self-consistency and - !! collision checks, so a stale or duplicate lookup entry is directly visible instead of inferred. - subroutine s_debug_write_ib_lookup_state(unit_num) - - integer, intent(in) :: unit_num - integer :: gi, gi2, jloc - - write (unit_num, '(A,I0)') 'local_ib_patch_ids: num_local_ibs=', num_local_ibs - do jloc = 1, num_local_ibs - write (unit_num, '(A,I0,A,I0,A,I0)') ' local_owner_idx=', jloc, ' -> patch_ib_idx=', local_ib_patch_ids(jloc), & - & ' gbl_patch_id=', patch_ib(local_ib_patch_ids(jloc))%gbl_patch_id - end do - - write (unit_num, '(A,I0,A,I0)') 'ib_gbl_idx_lookup non-negative entries: num_gbl_ibs=', num_gbl_ibs, ' num_ibs=', & - & num_ibs - do gi = 1, num_gbl_ibs - if (ib_gbl_idx_lookup(gi) > 0) then - write (unit_num, '(A,I0,A,I0)') ' gbl_id=', gi, ' -> local_idx=', ib_gbl_idx_lookup(gi) - end if - end do - - do jloc = 1, num_ibs - if (ib_gbl_idx_lookup(patch_ib(jloc)%gbl_patch_id) /= jloc) then - write (unit_num, '(A,I0,A,I0,A,I0)') ' INCONSISTENT: patch_ib(', jloc, ')%gbl_patch_id=', & - & patch_ib(jloc)%gbl_patch_id, ' but lookup points to local_idx=', & - & ib_gbl_idx_lookup(patch_ib(jloc)%gbl_patch_id) - end if - end do - do gi = 1, num_gbl_ibs - if (ib_gbl_idx_lookup(gi) <= 0) cycle - do gi2 = gi + 1, num_gbl_ibs - if (ib_gbl_idx_lookup(gi2) == ib_gbl_idx_lookup(gi)) then - write (unit_num, '(A,I0,A,I0,A,I0)') ' COLLISION: gbl_id=', gi, ' and gbl_id=', gi2, & - & ' both map to local_idx=', ib_gbl_idx_lookup(gi) - end if - end do - end do - - end subroutine s_debug_write_ib_lookup_state - subroutine s_get_neighborhood_idx(gbl_idx, neighborhood_idx) $:GPU_ROUTINE(parallelism='[seq]') diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp index 7881e193c..daf1c2a73 100644 --- a/src/simulation/m_time_steppers.fpp +++ b/src/simulation/m_time_steppers.fpp @@ -452,8 +452,6 @@ contains call cpu_time(start) call nvtxStartRange("TIMESTEP") - dbg_t_step = t_step ! TEMPORARY DEBUG INSTRUMENTATION: see m_ibm - ! Adaptive dt: initial stage if (adap_dt) call s_adaptive_dt_bubble(1) @@ -591,7 +589,6 @@ contains if (moving_immersed_boundary_flag) then call s_wrap_periodic_ibs() ! wraps the positions of IBs to the local proc call s_handoff_ib_ownership() ! recomputes which ranks own which IBs and communicate to neighbors - call s_debug_log_ib_divergence(t_step) ! TEMPORARY: logs any cross-rank IB state mismatch else if (ib_state_wrt) then call s_compute_ib_forces(q_prim_vf, fluid_pp) end if From 2683fcfa7264ab7af349fdaf99353812cb70b724 Mon Sep 17 00:00:00 2001 From: Daniel Vickers Date: Sun, 30 Aug 2026 14:52:19 -0400 Subject: [PATCH 10/10] I forgot to commit the docs and toolchain changes that make the hemispherical shepp direction work --- docs/documentation/case.md | 3 +- toolchain/mfc/case_validator.py | 49 ++++++++++++++++------------- toolchain/mfc/params/definitions.py | 1 + 3 files changed, 31 insertions(+), 22 deletions(-) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index 83f5a50fd..0972420f7 100644 --- a/docs/documentation/case.md +++ b/docs/documentation/case.md @@ -431,13 +431,14 @@ A particle cloud is a compact specification of a bed of identical circular (2D) | `cloud_geometry` | Integer | Shape of the cloud region. | | `shell_inner_radius` | Real | Inner radius for hemisphere-shell clouds (`cloud_geometry = 2`). | | `shell_outer_radius` | Real | Outer radius for hemisphere-shell clouds (`cloud_geometry = 2`). | +| `shell_axis` | Integer | Axis the hemisphere-shell cloud opens toward (`cloud_geometry = 2`). | | `moving_ibm` | Integer | Motion flag applied to every particle (see `patch_ib(j)%%moving_ibm`). | | `seed` | Integer | Random seed for reproducible placement (used by `packing_method = 1`). | | `packing_method` | Integer | Algorithm used to place the particles. | - `cloud_geometry` selects the cloud region: - `1` (box) uses `x[y,z]_centroid` and `length_x[y,z]` to define the region. - - `2` uses `x[y,z]_centroid`, `shell_inner_radius`, and `shell_outer_radius` to define a half-annulus in 2D and a hemisphere shell in 3D. Particle centres are sampled between `shell_inner_radius + radius` and `shell_outer_radius - radius`, and the flat plane is kept clear by one particle radius. The flat face is fixed at `y_centroid` in 2D and `z_centroid` in 3D; the filled region opens toward positive `y` in 2D and positive `z` in 3D. The full shell extent (`x[y,z]_centroid +/- shell_outer_radius` on the open side, and one particle radius of clearance on the flat-face side) must lie inside the computational domain; a hemisphere shell also requires at least two dimensions (`n > 0`). + - `2` uses `x[y,z]_centroid`, `shell_inner_radius`, `shell_outer_radius`, and `shell_axis` to define a half-annulus in 2D and a hemisphere shell in 3D. Particle centres are sampled between `shell_inner_radius + radius` and `shell_outer_radius - radius`, and the flat plane is kept clear by one particle radius. `shell_axis` (`1`=x, `2`=y, `3`=z; default `3`) selects which axis the shell opens toward from its flat face at that axis's centroid; in 2D there is no z-axis, so any value other than `1` opens toward `+y` (matching the fixed behavior before `shell_axis` existed). The open axis needs one particle radius of clearance on its flat-face side and the full `shell_outer_radius` on its open side; the other axis (2D) or two axes (3D) need the full shell extent (`centroid +/- shell_outer_radius`) inside the domain. A hemisphere shell also requires at least two dimensions (`n > 0`). - `packing_method` selects how the `num_particles` are positioned within the cloud region: - `1` (rejection sampling) draws random positions and rejects any that violate `min_spacing`, producing a disordered bed. `seed` makes the placement reproducible. - `2` (lattice) places the particles on the optimally dense lattice for the geometry — a triangular lattice in 2D and a face-centered cubic lattice in 3D. The lattice spacing is derived from the particle density (`num_particles` over the region area/volume); if that spacing is below the required `2*radius + min_spacing`, the region is too dense and the run aborts. diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 36604bac6..79293c353 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -785,6 +785,11 @@ def check_ibm(self): geometry == 2 and packing_method == 2, f"particle_cloud({i}) hemisphere-shell lattice packing is not implemented", ) + shell_axis = self.get(f"particle_cloud({i})%shell_axis", 3) + self.prohibit( + geometry == 2 and shell_axis not in [1, 2, 3], + f"particle_cloud({i})%shell_axis must be 1 (x), 2 (y), or 3 (z)", + ) if geometry == 2 and shell_outer_radius is not None and self._is_numeric(shell_outer_radius): x_centroid = self.get(f"particle_cloud({i})%x_centroid", None) y_centroid = self.get(f"particle_cloud({i})%y_centroid", None) @@ -795,32 +800,34 @@ def check_ibm(self): y_end = self.get("y_domain%end", None) z_beg = self.get("z_domain%beg", None) z_end = self.get("z_domain%end", None) - - if all(self._is_numeric(v) for v in [x_centroid, x_beg, x_end]): - self.prohibit( - x_centroid - shell_outer_radius < x_beg or x_centroid + shell_outer_radius > x_end, - f"particle_cloud({i}) hemisphere shell x-extent must lie within x_domain", - ) - if n > 0 and all(self._is_numeric(v) for v in [y_centroid, y_beg, y_end, radius]): - if p > 0: + # 2D has no z-axis; shell_axis values other than 1 (x) fall back to y, matching the + # fixed +y orientation used before shell_axis existed (see s_sample_cloud_candidate). + open_axis = shell_axis if (p > 0 or shell_axis == 1) else 2 + + axes = [ + (1, "x", x_centroid, x_beg, x_end), + (2, "y", y_centroid, y_beg, y_end), + (3, "z", z_centroid, z_beg, z_end), + ] + for axis_id, name, centroid, beg, end in axes: + if axis_id == 2 and n == 0: + continue + if axis_id == 3 and p == 0: + continue + if not all(self._is_numeric(v) for v in [centroid, beg, end, radius]): + continue + if axis_id == open_axis: + # the flat face sits at the centroid and the shell opens toward +axis; require + # one particle radius of standoff so no particle surface sits on the domain wall. self.prohibit( - y_centroid - shell_outer_radius < y_beg or y_centroid + shell_outer_radius > y_end, - f"particle_cloud({i}) hemisphere shell y-extent must lie within y_domain", + centroid - radius < beg or centroid + shell_outer_radius > end, + f"particle_cloud({i}) hemisphere shell must clear {name}_domain by one particle radius", ) else: - # 2D half-annulus opens toward +y from the flat face at y_centroid; require one - # particle radius of standoff so no particle surface sits on the domain wall. self.prohibit( - y_centroid - radius < y_beg or y_centroid + shell_outer_radius > y_end, - f"particle_cloud({i}) half-annulus must clear y_domain by one particle radius", + centroid - shell_outer_radius < beg or centroid + shell_outer_radius > end, + f"particle_cloud({i}) hemisphere shell {name}-extent must lie within {name}_domain", ) - if p > 0 and all(self._is_numeric(v) for v in [z_centroid, z_beg, z_end, radius]): - # 3D hemisphere shell opens toward +z from the flat face at z_centroid; require one - # particle radius of standoff so no particle surface sits on the domain wall. - self.prohibit( - z_centroid - radius < z_beg or z_centroid + shell_outer_radius > z_end, - f"particle_cloud({i}) hemisphere shell must clear z_domain by one particle radius", - ) num_ib_airfoils_max = get_fortran_constants().get("num_ib_airfoils_max", 5) num_stl_models_max = get_fortran_constants().get("num_stl_models_max", 10) diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 5f218f5d0..b233ad8e9 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -1011,6 +1011,7 @@ def _load(): _pb_attrs["moving_ibm"] = (INT, _pb_tags) _pb_attrs["seed"] = (INT, _pb_tags) _pb_attrs["cloud_geometry"] = (INT, _pb_tags) + _pb_attrs["shell_axis"] = (INT, _pb_tags) _pb_attrs["packing_method"] = (INT, _pb_tags) _pb_attrs["periodic"] = (INT, _pb_tags) REGISTRY.register_family(