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/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/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_data_output.fpp b/src/simulation/m_data_output.fpp index 22bde22e1..f013dd92a 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,135 @@ 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) 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 + + ! 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) + 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 + !> 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_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_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_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 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 7d54dba97..32fced2bc 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 @@ -957,6 +1021,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 @@ -1102,10 +1167,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() @@ -1289,7 +1356,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(:) @@ -1311,11 +1378,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]') @@ -1331,16 +1400,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() @@ -1358,11 +1429,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]') @@ -1377,12 +1450,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() @@ -1443,6 +1518,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 @@ -1486,9 +1563,9 @@ contains ! Post all receives first, then sends 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) @@ -1503,9 +1580,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) @@ -1519,7 +1596,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) @@ -1531,12 +1608,13 @@ 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 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 @@ -1555,20 +1633,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 ib_gbl_idx_lookup = -1 - $:GPU_UPDATE(device='[ib_gbl_idx_lookup]') - - $:GPU_PARALLEL_LOOP(private='[i]') do i = 1, num_ibs ib_gbl_idx_lookup(patch_ib(i)%gbl_patch_id) = i end do - $:END_GPU_PARALLEL_LOOP() - $:GPU_UPDATE(host='[ib_gbl_idx_lookup]') + $:GPU_UPDATE(device='[ib_gbl_idx_lookup]') end subroutine s_update_ib_lookup 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 6e0843352..2fb07b1ee 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 @@ -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 4d10417ed..fb941098b 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,50 @@ 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. 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, 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 + + n_total_particles = 0 + do cloud_idx = 1, num_particle_clouds + n_total_particles = n_total_particles + particle_cloud(cloud_idx)%num_particles + end do + 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 + 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 + num_particle_cloud_ibs = ib_idx + + 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 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(