Debug ibm stability

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.

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)

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

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

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

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

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)