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
@@ -312,6 +312,7 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | Parameter            | Type    | Description |
 | ---:                 | :----:  | :---                |
 | `num_ibs`            | Integer | Number of immersed boundary patches |
+| `ib_neighborhood_radius` | Integer | Parameter that controls the neighborhood size for IB detection. |
 | `geometry`           | Integer | Geometry configuration of the patch.|
 | `x[y,z]_centroid`    | Real    | Centroid of the applied geometry in the [x,y,z]-direction. |
 | `length_x[y,z]`      | Real    | Length, if applicable, in the [x,y,z]-direction. |
@@ -373,6 +374,8 @@ Additional details on this specification can be found in [NACA airfoil](https://
 
 - `ib_coefficient_of_friction` is the coefficient of friction used in IB collisions.
 
+- `ib_neighborhood_radius` controls the size of the neighborhood size. This value defaults to 1, which indicates that any given rank is aware of IB's up to 1 ranks away. This parameter is required to strong-scale a case when IB's eventually grow to be larger than one full processor domain wide.
+
 ### 5. Fluid Material's {#sec-fluid-materials}
 
 | Parameter | Type   | Description                                    |
@@ -638,7 +641,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 +709,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.

@@ -87,7 +87,7 @@
             "precision": 2,
             "prim_vars_wrt": "T",
             "E_wrt": "T",
-            "ib_state_wrt": "T",
+            "ib_state_wrt": "F",
             "parallel_io": "T",
             # Patch: Constant Tube filled with air
             # Specify the cylindrical air tube grid geometry

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,12 +983,13 @@ contains
         dx_local = minval(dx); dy_local = minval(dy)
         if (p /= 0) dz_local = minval(dz)
 
+        num_gbl_ibs = num_ibs
         allocate (stl_bounding_boxes(num_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
                 allocate (models(patch_id)%model)
-                print *, " * Reading model: " // trim(patch_ib(patch_id)%model_filepath)
+                if (proc_rank == 0) print *, " * Reading model: " // trim(patch_ib(patch_id)%model_filepath)
 
                 model = f_model_read(patch_ib(patch_id)%model_filepath)
                 params%scale(:) = patch_ib(patch_id)%model_scale(:)
@@ -1001,9 +1002,7 @@ contains
                     params%scale(:) = 1._wp
                 end if
 
-                if (proc_rank == 0) then
-                    print *, " * Transforming model."
-                end if
+                if (proc_rank == 0) print *, " * Transforming model."
 
                 ! Get the model center before transforming the model
                 bbox_old = f_create_bbox(model)

diff --git a/src/common/m_mpi_common.fpp b/src/common/m_mpi_common.fpp
@@ -94,6 +94,8 @@ contains
         proc_rank = 0
 #endif
 
+        $:GPU_UPDATE(device='[num_procs, proc_rank]')
+
     end subroutine s_mpi_initialize
 
     !> Set up MPI I/O data views and variable pointers for parallel file output.

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)