Debug ibm stability

docs/documentation/case.md

@@ -456,13 +456,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

@@ -391,6 +391,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
@@ -520,6 +521,7 @@ module m_derived_types
         real(wp), dimension(3)       :: ip_loc         !< Physical location of the image point
         integer, dimension(3)        :: ip_grid        !< Top left grid point of IP
         real(wp), dimension(2, 2, 2) :: interp_coeffs  !< Interpolation Coefficients of image point
+        logical                      :: interp_valid   !< .false. if every image point stencil cell lies inside an IB
         integer                      :: ib_patch_id    !< ID of the IB Patch the ghost point is part of
         real(wp)                     :: levelset
         real(wp), dimension(1:3)     :: levelset_norm

src/common/m_finite_differences.fpp

@@ -30,14 +30,29 @@ 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
 
+        ! The stencil below reaches fd_number_in cells either side of i, and s_cc only exists over the
+        ! caller's buffer, so the loop cannot start closer than fd_number_in to either end of it. Without
+        ! this, a caller whose buffer is narrower than 2*fd_number_in reads off the front of s_cc: a
+        ! 4th-order hypoelastic case with buff_size = 2 starts at i = -2 and immediately asks for
+        ! s_cc(-4). Immersed-boundary cases are unaffected -- s_mfc_buff_size gives them buff_size >= 10,
+        ! so the clamp never binds and they still get coefficients out to -fd_number_in, which is the
+        ! whole point of widening the range above.
+        lB = max(lB, -local_buff_size + fd_number_in)
+        lE = min(lE, q + local_buff_size - fd_number_in)
+
         ! Computing the 1st order finite-difference coefficients
         if (fd_order_in == 1) then
             do i = lB, lE

src/post_process/m_data_output.fpp

@@ -454,31 +454,54 @@ contains
             ! Finally, the local quadrilateral mesh, either 2D or 3D, along with its offsets that indicate the presence and size of
             ! ghost zone layer(s), are put in the formatted database slave file.
 
+            ! Silo carries the mesh coordinates in their own datatype, separate from the flow variables, so the cell boundaries
+            ! are copied down to single precision here when that is what was asked for. Without this the mesh is always written as
+            ! DB_DOUBLE, which keeps downstream readers on a double-precision path regardless of `precision`.
+            if (precision == precision_single) then
+                x_cb_s = real(x_cb, sp)
+                if (n > 0) then
+                    y_cb_s = real(y_cb, sp)
+                    if (p > 0) z_cb_s = real(z_cb, sp)
+                end if
+            end if
+
             if (p > 0) then
                 err = DBMKOPTLIST(2, out%optlist)
                 err = DBADDIAOPT(out%optlist, DBOPT_LO_OFFSET, size(out%lo_offset), out%lo_offset)
                 err = DBADDIAOPT(out%optlist, DBOPT_HI_OFFSET, size(out%hi_offset), out%hi_offset)
-                if (grid_geometry == 3) then
-                    err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, y_cb, z_cb, x_cb, out%dims, 3, &
-                                  & DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr)
-                else
-                    err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, z_cb, out%dims, 3, &
-                                  & DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr)
-                end if
+                #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")]
+                    if (precision == ${PRECISION}$) then
+                        if (grid_geometry == 3) then
+                            err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, y_cb${SFX}$, z_cb${SFX}$, &
+                                          & x_cb${SFX}$, out%dims, 3, ${DBT}$, DB_COLLINEAR, out%optlist, ierr)
+                        else
+                            err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb${SFX}$, y_cb${SFX}$, &
+                                          & z_cb${SFX}$, out%dims, 3, ${DBT}$, DB_COLLINEAR, out%optlist, ierr)
+                        end if
+                    end if
+                #:endfor
                 err = DBFREEOPTLIST(out%optlist)
             else if (n > 0) then
                 err = DBMKOPTLIST(2, out%optlist)
                 err = DBADDIAOPT(out%optlist, DBOPT_LO_OFFSET, size(out%lo_offset), out%lo_offset)
                 err = DBADDIAOPT(out%optlist, DBOPT_HI_OFFSET, size(out%hi_offset), out%hi_offset)
-                err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, y_cb, DB_F77NULL, out%dims, 2, &
-                              & DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr)
+                #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")]
+                    if (precision == ${PRECISION}$) then
+                        err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb${SFX}$, y_cb${SFX}$, &
+                                      & DB_F77NULL, out%dims, 2, ${DBT}$, DB_COLLINEAR, out%optlist, ierr)
+                    end if
+                #:endfor
                 err = DBFREEOPTLIST(out%optlist)
             else
                 err = DBMKOPTLIST(2, out%optlist)
                 err = DBADDIAOPT(out%optlist, DBOPT_LO_OFFSET, size(out%lo_offset), out%lo_offset)
                 err = DBADDIAOPT(out%optlist, DBOPT_HI_OFFSET, size(out%hi_offset), out%hi_offset)
-                err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb, DB_F77NULL, DB_F77NULL, out%dims, &
-                              & 1, DB_DOUBLE, DB_COLLINEAR, out%optlist, ierr)
+                #:for PRECISION, SFX, DBT in [(1,'_s','DB_FLOAT'),(2,'',"DB_DOUBLE")]
+                    if (precision == ${PRECISION}$) then
+                        err = DBPUTQM(out%dbfile, 'rectilinear_grid', 16, 'x', 1, 'y', 1, 'z', 1, x_cb${SFX}$, DB_F77NULL, &
+                                      & DB_F77NULL, out%dims, 1, ${DBT}$, DB_COLLINEAR, out%optlist, ierr)
+                    end if
+                #:endfor
                 err = DBFREEOPTLIST(out%optlist)
             end if
         else if (format == format_binary) then
@@ -1349,6 +1372,10 @@ contains
         real(wp), dimension(:), allocatable             :: omega_x, omega_y, omega_z
         real(wp), dimension(:), allocatable             :: angle_x, angle_y, angle_z
         real(wp), dimension(:), allocatable             :: ib_diameter
+        real(sp), dimension(:), allocatable             :: px_s, py_s, pz_s
+        logical, dimension(:), allocatable              :: keep
+        integer                                         :: nKept
+        real(wp)                                        :: r_ib
 
         if (proc_rank == 0) then
             nBodies = num_ibs
@@ -1419,12 +1446,51 @@ contains
                     ib_diameter(i) = ib_data(i, 20)*2.0_wp
                 end do
 
+                ! When only part of the domain is written, the bodies outside that window are dropped so the point mesh matches
+                ! the cropped grid. A body is kept when its bounding sphere overlaps the window rather than when its centroid is
+                ! inside it, so one straddling the boundary still appears instead of vanishing at the edge.
+                if (output_partial_domain) then
+                    allocate (keep(nBodies))
+
+                    do i = 1, nBodies
+                        r_ib = 0.5_wp*ib_diameter(i)
+                        keep(i) = (px(i) + r_ib >= x_output%beg) .and. (px(i) - r_ib <= x_output%end)
+                        if (n > 0) keep(i) = keep(i) .and. (py(i) + r_ib >= y_output%beg) .and. (py(i) - r_ib <= y_output%end)
+                        if (p > 0) keep(i) = keep(i) .and. (pz(i) + r_ib >= z_output%beg) .and. (pz(i) - r_ib <= z_output%end)
+                    end do
+
+                    nKept = count(keep)
+
+                    if (nKept < nBodies) then
+                        #:for A in ['px','py','pz','ib_diameter','force_x','force_y','force_z','torque_x','torque_y','torque_z']
+                            ${A}$(1:nKept) = pack(${A}$(1:nBodies), keep)
+                        #:endfor
+                        #:for A in ['vel_x','vel_y','vel_z','omega_x','omega_y','omega_z','angle_x','angle_y','angle_z']
+                            ${A}$(1:nKept) = pack(${A}$(1:nBodies), keep)
+                        #:endfor
+                    end if
+
+                    nBodies = nKept
+                    deallocate (keep)
+                end if
+
                 write (meshnames(1), '(A,I0,A)') '../p0/', t_step, '.silo:ib_bodies'
                 meshtypes(1) = DB_POINTMESH
                 err = DBSET2DSTRLEN(len(meshnames(1)))
                 err = DBPUTMMESH(out%dbroot, 'ib_bodies', 16, 1, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr)
 
-                err = DBPUTPM(out%dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
+                ! Silo carries the point-mesh coordinates in their own datatype, so they need the same single-precision
+                ! treatment as the rectilinear mesh.
+                if (precision == precision_single) then
+                    allocate (px_s(nBodies), py_s(nBodies), pz_s(nBodies))
+                    px_s(1:nBodies) = real(px(1:nBodies), sp)
+                    py_s(1:nBodies) = real(py(1:nBodies), sp)
+                    pz_s(1:nBodies) = real(pz(1:nBodies), sp)
+                    err = DBPUTPM(out%dbfile, 'ib_bodies', 9, 3, px_s, py_s, pz_s, nBodies, DB_FLOAT, DB_F77NULL, ierr)
+                    deallocate (px_s, py_s, pz_s)
+                else
+                    err = DBPUTPM(out%dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
+                end if
 
                 call s_write_ib_variable('ib_force_x', t_step, force_x, nBodies)
                 call s_write_ib_variable('ib_force_y', t_step, force_y, nBodies)
@@ -1456,20 +1522,29 @@ contains
     !> Write a single IB point-variable to the Silo database slave and master files.
     subroutine s_write_ib_variable(varname, t_step, data, nBodies)
 
-        character(len=*), intent(in)       :: varname
-        integer, intent(in)                :: t_step
-        real(wp), dimension(:), intent(in) :: data
-        integer, intent(in)                :: nBodies
-        character(len=4*name_len)          :: var_name_entry
-        integer                            :: var_type_entry, ierr
+        character(len=*), intent(in)        :: varname
+        integer, intent(in)                 :: t_step
+        real(wp), dimension(:), intent(in)  :: data
+        integer, intent(in)                 :: nBodies
+        character(len=4*name_len)           :: var_name_entry
+        integer                             :: var_type_entry, ierr
+        real(sp), dimension(:), allocatable :: data_s
 
         write (var_name_entry, '(A,I0,A)') '../p0/', t_step, '.silo:' // trim(varname)
         var_type_entry = DB_POINTVAR
         err = DBSET2DSTRLEN(len(var_name_entry))
         err = DBPUTMVAR(out%dbroot, trim(varname), len_trim(varname), 1, var_name_entry, len_trim(var_name_entry), &
                         & var_type_entry, DB_F77NULL, ierr)
 
-        err = DBPUTPV1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
+        if (precision == precision_single) then
+            allocate (data_s(nBodies))
+            data_s(1:nBodies) = real(data(1:nBodies), sp)
+            err = DBPUTPV1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data_s, nBodies, DB_FLOAT, DB_F77NULL, &
+                           & ierr)
+            deallocate (data_s)
+        else
+            err = DBPUTPV1(out%dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
+        end if
 
     end subroutine s_write_ib_variable
 

src/post_process/m_derived_variables.fpp

@@ -33,18 +33,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