MFlowCode · danieljvickers · Mar 15, 2026 · Mar 15, 2026 · Mar 17, 2026 · Mar 17, 2026
diff --git a/docs/documentation/case.md b/docs/documentation/case.md
@@ -638,7 +638,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca
 | `alpha_wrt(i)`          | Logical | Add the volume fraction of fluid $i$ to the database	|
 | `gamma_wrt`             | Logical | Add the specific heat ratio function to the database	|
 | `heat_ratio_wrt`        | Logical | Add the specific heat ratio to the database	|
-| `ib_state_wrt`          | Logical | Write IB state and loads to a datafile at each time step |
+| `ib_state_wrt`          | Logical | Parameter to handle writing IB state on saves and outputting the state as a point mesh to SILO files. |
 | `pi_inf_wrt`            | Logical | Add the liquid stiffness function to the database |
 | `pres_inf_wrt`          | Logical | Add the liquid stiffness to the formatted database	 |
 | `c_wrt`                 | Logical | Add the sound speed to the database	 |
@@ -706,7 +706,7 @@ If `file_per_process` is true, then pre_process, simulation, and post_process mu
 
 - `probe_wrt` activates the output of state variables at coordinates specified by `probe(i)%[x;y,z]`.
 
-- `ib_state_wrt` activates the output of data specified by patch_ib(i)%force(:) (and torque, vel, angular_vel, angles, [x,y,z]_centroid) into a single binary datafile for all IBs at all timesteps. During post_processing, this file is converted into separate time histories for each IB.
+- `ib_state_wrt` is used to trigger post-processing of the IB state to be written out as a point mesh in the SILO files. When no IBs are moving, it also triggers force and torque calculation so that those values may be written to the output state files.
 
 - `output_partial_domain` activates the output of part of the domain specified by `[x,y,z]_output%%beg` and `[x,y,z]_output%%end`.
 This is useful for large domains where only a portion of the domain is of interest.

diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp
@@ -25,6 +25,7 @@ module m_constants
     integer, parameter  :: num_probes_max = 10                 !< Maximum number of flow probes in the simulation
     integer, parameter  :: num_patches_max = 10                !< Maximum number of IC patches
     integer, parameter  :: num_ib_patches_max = 50000          !< Maximum number of immersed boundary patches (patch_ib)
+    integer, parameter  :: num_local_ibs_max = 2000            !< Maximum number of immersed boundary patches (patch_ib)
     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)
     integer, parameter  :: max_sph_harm_degree = 5             !< Max degree L for 3D spherical harmonic patch (geometry 14)

diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp
@@ -202,12 +202,10 @@ module m_derived_types
 
     type :: t_model_array
         ! Original CPU-side fields (unchanged)
-        type(t_model), allocatable              :: model                    !< STL/OBJ geometry model
-        real(wp), allocatable, dimension(:,:,:) :: boundary_v               !< Boundary vertices
-        real(wp), allocatable, dimension(:,:)   :: interpolated_boundary_v  !< Interpolated boundary vertices
-        integer                                 :: boundary_edge_count      !< Number of boundary edges
-        integer                                 :: total_vertices           !< Total vertex count
-        integer                                 :: interpolate              !< Interpolation flag
+        type(t_model), allocatable              :: model                !< STL/OBJ geometry model
+        real(wp), allocatable, dimension(:,:,:) :: boundary_v           !< Boundary vertices
+        integer                                 :: boundary_edge_count  !< Number of boundary edges
+        integer                                 :: total_vertices       !< Total vertex count
 
         ! GPU-friendly flattened arrays
         integer                                 :: ntrs   !< Copy of model%ntrs
@@ -272,9 +270,10 @@ module m_derived_types
     end type ic_patch_parameters
 
     type ib_patch_parameters
-
         integer  :: geometry                            !< Type of geometry for the patch
+        integer  :: gbl_patch_id
         real(wp) :: x_centroid, y_centroid, z_centroid  !< Geometric center coordinates of the patch
+
         !> Centroid locations of intermediate steps in the time_stepper module
         real(wp)                 :: step_x_centroid, step_y_centroid, step_z_centroid
         real(wp), dimension(1:3) :: centroid_offset  !< offset of center of mass from computed cell center for odd-shaped IBs

diff --git a/src/common/m_model.fpp b/src/common/m_model.fpp
@@ -983,10 +983,11 @@ contains
         dx_local = minval(dx); dy_local = minval(dy)
         if (p /= 0) dz_local = minval(dz)
 
-        allocate (stl_bounding_boxes(num_ibs,1:3,1:3))
+        allocate (stl_bounding_boxes(num_gbl_ibs,1:3,1:3))
 
         do patch_id = 1, num_ibs
             if (patch_ib(patch_id)%geometry == 5 .or. patch_ib(patch_id)%geometry == 12) then
+                print *, proc_rank, patch_id, num_ibs, patch_ib(patch_id)%geometry
                 allocate (models(patch_id)%model)
                 print *, " * Reading model: " // trim(patch_ib(patch_id)%model_filepath)
 

diff --git a/src/post_process/m_data_output.fpp b/src/post_process/m_data_output.fpp
@@ -1396,48 +1396,42 @@ contains
                 end do
 
                 close (file_unit)
-            end if
-
-            call MPI_BCAST(ib_data, nBodies*NFIELDS_PER_IB, mpi_p, 0, MPI_COMM_WORLD, ierr)
-
-            do i = 1, nBodies
-                force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4)
-                torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7)
-                vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10)
-                omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13)
-                angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16)
-                px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19)
-                ib_diameter(i) = ib_data(i, 20)*2.0_wp
-            end do
 
-            if (proc_rank == 0) then
-                do i = 1, num_procs
-                    write (meshnames(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:ib_bodies'
-                    meshtypes(i) = DB_POINTMESH
+                do i = 1, nBodies
+                    force_x(i) = ib_data(i, 2); force_y(i) = ib_data(i, 3); force_z(i) = ib_data(i, 4)
+                    torque_x(i) = ib_data(i, 5); torque_y(i) = ib_data(i, 6); torque_z(i) = ib_data(i, 7)
+                    vel_x(i) = ib_data(i, 8); vel_y(i) = ib_data(i, 9); vel_z(i) = ib_data(i, 10)
+                    omega_x(i) = ib_data(i, 11); omega_y(i) = ib_data(i, 12); omega_z(i) = ib_data(i, 13)
+                    angle_x(i) = ib_data(i, 14); angle_y(i) = ib_data(i, 15); angle_z(i) = ib_data(i, 16)
+                    px(i) = ib_data(i, 17); py(i) = ib_data(i, 18); pz(i) = ib_data(i, 19)
+                    ib_diameter(i) = ib_data(i, 20)*2.0_wp
                 end do
+
+                write (meshnames(1), '(A,I0,A)') '../p0/', t_step, '.silo:ib_bodies'
+                meshtypes(1) = DB_POINTMESH
                 err = DBSET2DSTRLEN(len(meshnames(1)))
-                err = DBPUTMMESH(dbroot, 'ib_bodies', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr)
+                err = DBPUTMMESH(dbroot, 'ib_bodies', 16, 1, meshnames, len_trim(meshnames), meshtypes, DB_F77NULL, ierr)
+
+                err = DBPUTPM(dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
+
+                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)
+                call s_write_ib_variable('ib_force_z', t_step, force_z, nBodies)
+                call s_write_ib_variable('ib_torque_x', t_step, torque_x, nBodies)
+                call s_write_ib_variable('ib_torque_y', t_step, torque_y, nBodies)
+                call s_write_ib_variable('ib_torque_z', t_step, torque_z, nBodies)
+                call s_write_ib_variable('ib_vel_x', t_step, vel_x, nBodies)
+                call s_write_ib_variable('ib_vel_y', t_step, vel_y, nBodies)
+                call s_write_ib_variable('ib_vel_z', t_step, vel_z, nBodies)
+                call s_write_ib_variable('ib_omega_x', t_step, omega_x, nBodies)
+                call s_write_ib_variable('ib_omega_y', t_step, omega_y, nBodies)
+                call s_write_ib_variable('ib_omega_z', t_step, omega_z, nBodies)
+                call s_write_ib_variable('ib_angle_x', t_step, angle_x, nBodies)
+                call s_write_ib_variable('ib_angle_y', t_step, angle_y, nBodies)
+                call s_write_ib_variable('ib_angle_z', t_step, angle_z, nBodies)
+                call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nBodies)
             end if
 
-            err = DBPUTPM(dbfile, 'ib_bodies', 9, 3, px, py, pz, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
-
-            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)
-            call s_write_ib_variable('ib_force_z', t_step, force_z, nBodies)
-            call s_write_ib_variable('ib_torque_x', t_step, torque_x, nBodies)
-            call s_write_ib_variable('ib_torque_y', t_step, torque_y, nBodies)
-            call s_write_ib_variable('ib_torque_z', t_step, torque_z, nBodies)
-            call s_write_ib_variable('ib_vel_x', t_step, vel_x, nBodies)
-            call s_write_ib_variable('ib_vel_y', t_step, vel_y, nBodies)
-            call s_write_ib_variable('ib_vel_z', t_step, vel_z, nBodies)
-            call s_write_ib_variable('ib_omega_x', t_step, omega_x, nBodies)
-            call s_write_ib_variable('ib_omega_y', t_step, omega_y, nBodies)
-            call s_write_ib_variable('ib_omega_z', t_step, omega_z, nBodies)
-            call s_write_ib_variable('ib_angle_x', t_step, angle_x, nBodies)
-            call s_write_ib_variable('ib_angle_y', t_step, angle_y, nBodies)
-            call s_write_ib_variable('ib_angle_z', t_step, angle_z, nBodies)
-            call s_write_ib_variable('ib_diameter', t_step, ib_diameter, nBodies)
-
             deallocate (ib_data, px, py, pz, force_x, force_y, force_z)
             deallocate (torque_x, torque_y, torque_z, vel_x, vel_y, vel_z)
             deallocate (omega_x, omega_y, omega_z, angle_x, angle_y, angle_z)
@@ -1450,23 +1444,18 @@ 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), dimension(num_procs) :: var_names
-        integer, dimension(num_procs)                   :: var_types
-        integer                                         :: ierr, i
-
-        if (proc_rank == 0) then
-            do i = 1, num_procs
-                write (var_names(i), '(A,I0,A,I0,A)') '../p', i - 1, '/', t_step, '.silo:' // trim(varname)
-                var_types(i) = DB_POINTVAR
-            end do
-            err = DBSET2DSTRLEN(len(var_names(1)))
-            err = DBPUTMVAR(dbroot, trim(varname), len_trim(varname), num_procs, var_names, len_trim(var_names), var_types, &
-                            & DB_F77NULL, ierr)
-        end if
+        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
+
+        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(dbroot, trim(varname), len_trim(varname), 1, var_name_entry, len_trim(var_name_entry), var_type_entry, &
+                        & DB_F77NULL, ierr)
 
         err = DBPUTPV1(dbfile, trim(varname), len_trim(varname), 'ib_bodies', 9, data, nBodies, DB_DOUBLE, DB_F77NULL, ierr)
 

diff --git a/src/simulation/m_collisions.fpp b/src/simulation/m_collisions.fpp
@@ -19,7 +19,8 @@ module m_collisions
 
     implicit none
 
-    private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module
+    private; public :: s_apply_collision_forces, s_initialize_collisions_module, s_finalize_collisions_module, &
+        & f_local_rank_owns_collision
     ! overlap distances for computing collisions
     integer, allocatable, dimension(:,:)  :: collision_lookup
     real(wp), allocatable, dimension(:,:) :: wall_overlap_distances
@@ -405,35 +406,35 @@ contains
         logical                            :: owns_collision
         real(wp), dimension(3)             :: projected_location
 
-        #:if defined('MFC_MPI')
-            if (num_procs == 1) then
-                owns_collision = .true.
-            else
-                projected_location(:) = collision_location(:)
-
-                ! catch the edge case where th collision lies just outside the computational domain
-                #:for X, ID in [('x', 1), ('y', 2), ('z', 3)]
-                    if (num_dims >= ${ID}$) then
-                        if (ib_bc_${X}$%beg /= BC_PERIODIC) then
-                            ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
-                            if (collision_location(${ID}$) < ${X}$_domain%beg) then
-                                projected_location(${ID}$) = ${X}$_domain%beg
-                            else if (${X}$_domain%end < collision_location(${ID}$)) then
-                                projected_location(${ID}$) = ${X}$_domain%end - 1.0e-10_wp
-                            end if
+#ifdef MFC_MPI
+        if (num_procs == 1) then
+            owns_collision = .true.
+        else
+            projected_location(:) = collision_location(:)
+
+            ! catch the edge case where th collision lies just outside the computational domain
+            #:for X, ID in [('x', 1), ('y', 2), ('z', 3)]
+                if (num_dims >= ${ID}$) then
+                    if (ib_bc_${X}$%beg /= BC_PERIODIC) then
+                        ! if it is outside the domain in one direction, project it somewhere inside so at least one rank owns it
+                        if (collision_location(${ID}$) < ${X}$_domain%beg) then
+                            projected_location(${ID}$) = ${X}$_domain%beg
+                        else if (${X}$_domain%end < collision_location(${ID}$)) then
+                            projected_location(${ID}$) = ${X}$_domain%end - 1.0e-10_wp
                         end if
                     end if
-                #:endfor
+                end if
+            #:endfor
 
-                ! the object that contains the collision location owns the collisions
-                owns_collision = x_cb(-1) <= projected_location(1) .and. projected_location(1) < x_cb(m)
-                owns_collision = owns_collision .and. y_cb(-1) <= projected_location(2) .and. projected_location(2) < y_cb(n)
-                if (num_dims == 3) owns_collision = owns_collision .and. z_cb(-1) <= projected_location(3) &
-                    & .and. projected_location(3) < z_cb(p)
-            end if
-        #:else
-            owns_collision = .true.
-        #:endif
+            ! the object that contains the collision location owns the collisions
+            owns_collision = x_cb(-1) <= projected_location(1) .and. projected_location(1) < x_cb(m)
+            owns_collision = owns_collision .and. y_cb(-1) <= projected_location(2) .and. projected_location(2) < y_cb(n)
+            if (num_dims == 3) owns_collision = owns_collision .and. z_cb(-1) <= projected_location(3) .and. projected_location(3) &
+                & < z_cb(p)
+        end if
+#else
+        owns_collision = .true.
+#endif
 
     end function f_local_rank_owns_collision
 

diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
@@ -909,7 +909,7 @@ contains
         integer                              :: ifile, ierr
         integer, dimension(MPI_STATUS_SIZE)  :: status
         logical                              :: file_exist
-        integer                              :: i
+        integer                              :: i, ib_idx
         integer, parameter                   :: NFIELDS_PER_IB = 20
         real(wp)                             :: ib_buf(NFIELDS_PER_IB)
 
@@ -923,19 +923,6 @@ contains
         end if
         call s_mpi_barrier()
 
-        ! Divide num_ibs across num_procs
-        nibs_per_rank = num_ibs/num_procs
-        remainder = mod(num_ibs, num_procs)
-
-        ! Ranks < remainder get one extra IB
-        if (proc_rank < remainder) then
-            ib_start = proc_rank*(nibs_per_rank + 1) + 1
-            ib_end = ib_start + nibs_per_rank  ! nibs_per_rank + 1 total
-        else
-            ib_start = remainder*(nibs_per_rank + 1) + (proc_rank - remainder)*nibs_per_rank + 1
-            ib_end = ib_start + nibs_per_rank - 1
-        end if
-
         write (file_loc, '(A,I0,A)') '/restart_data/ib_state_', t_step, '.dat'
         file_loc = trim(case_dir) // trim(file_loc)
 
@@ -947,20 +934,21 @@ contains
 
         call MPI_FILE_OPEN(MPI_COMM_WORLD, file_loc, ior(MPI_MODE_WRONLY, MPI_MODE_CREATE), mpi_info_int, ifile, ierr)
 
-        do i = ib_start, ib_end
+        do i = 1, num_local_ibs
+            ib_idx = local_ib_patch_ids(i)
             ib_buf(1) = mytime
-            ib_buf(2:4) = patch_ib(i)%force(1:3)
-            ib_buf(5:7) = patch_ib(i)%torque(1:3)
-            ib_buf(8:10) = patch_ib(i)%vel(1:3)
-            ib_buf(11:13) = patch_ib(i)%angular_vel(1:3)
-            ib_buf(14:16) = patch_ib(i)%angles(1:3)
-            ib_buf(17) = patch_ib(i)%x_centroid
-            ib_buf(18) = patch_ib(i)%y_centroid
-            ib_buf(19) = patch_ib(i)%z_centroid
-            ib_buf(20) = patch_ib(i)%radius
+            ib_buf(2:4) = patch_ib(ib_idx)%force(1:3)
+            ib_buf(5:7) = patch_ib(ib_idx)%torque(1:3)
+            ib_buf(8:10) = patch_ib(ib_idx)%vel(1:3)
+            ib_buf(11:13) = patch_ib(ib_idx)%angular_vel(1:3)
+            ib_buf(14:16) = patch_ib(ib_idx)%angles(1:3)
+            ib_buf(17) = patch_ib(ib_idx)%x_centroid
+            ib_buf(18) = patch_ib(ib_idx)%y_centroid
+            ib_buf(19) = patch_ib(ib_idx)%z_centroid
+            ib_buf(20) = patch_ib(ib_idx)%radius
 
             ! Global IB index (i) determines position in file
-            disp = int(i - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK
+            disp = int(patch_ib(ib_idx)%gbl_patch_id - 1, MPI_OFFSET_KIND)*int(NFIELDS_PER_IB, MPI_OFFSET_KIND)*WP_MOK
 
             call MPI_FILE_WRITE_AT(ifile, disp, ib_buf, NFIELDS_PER_IB, mpi_p, status, ierr)
         end do