Local Aware IBM

docs/documentation/case.md

@@ -312,6 +312,8 @@ This is enabled by adding ``'elliptic_smoothing': "T",`` and ``'elliptic_smoothi
 | Parameter            | Type    | Description |
 | ---:                 | :----:  | :---                |
 | `num_ibs`            | Integer | Number of immersed boundary patches |
+| `num_particle_beds`  | Integer | Number of particle bed specifications to generate immersed boundary patches from |
+| `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 +375,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                                    |
@@ -643,7 +647,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	 |
@@ -711,7 +715,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.

docs/module_categories.json

@@ -38,6 +38,7 @@
             "m_compute_cbc",
             "m_boundary_common",
             "m_ibm",
+            "m_particle_bed",
             "m_igr",
             "m_ib_patches",
             "m_compute_levelset"

examples/2D_mibm_particle_bed/case.py

@@ -0,0 +1,136 @@
+import json
+import math
+
+# 2D shock wave interacting with a bed of 20 free-floating circular particles.
+
+gam_a = 1.4
+
+# Shock parameters (Mach 1.5)
+mach_number = 1.5
+pre_shock_pressure = 1
+pre_shock_density = 1.4
+pre_shock_speed = 0.0
+post_shock_pressure = 2.4583
+post_shock_density = 2.6069
+post_shock_speed = 0.6944
+
+domain_size = 4.0
+wave_front = -1.5
+
+total_time = 1.5
+num_time_steps = 2000
+dt = float(total_time / num_time_steps)
+num_saves = 100
+steps_to_save = int(num_time_steps / num_saves)
+
+# Soft-sphere collision parameters (from 3D_mibm_sphere_head_on_collision)
+collision_time = 20.0 * dt
+
+# Particle bed parameters
+bed_x = 0.5
+bed_y = 0.0
+bed_lx = 2.0
+bed_ly = 3.5
+particle_radius = 0.15
+particle_mass = 0.25
+particle_min_spacing = 0.05
+
+print(
+    json.dumps(
+        {
+            # Logistics
+            "run_time_info": "T",
+            # Computational Domain Parameters
+            "x_domain%beg": -domain_size * 0.5,
+            "x_domain%end": domain_size * 0.5,
+            "y_domain%beg": -domain_size * 0.5,
+            "y_domain%end": domain_size * 0.5,
+            "cyl_coord": "F",
+            "m": 256,
+            "n": 256,
+            "p": 0,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": num_time_steps,
+            "t_step_save": steps_to_save,
+            # Simulation Algorithm Parameters
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "T",
+            "time_stepper": 3,
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "weno_Re_flux": "T",
+            "weno_avg": "T",
+            "avg_state": 2,
+            "mapped_weno": "T",
+            "null_weights": "F",
+            "mp_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "bc_x%beg": -17,
+            "bc_x%end": -8,
+            "bc_y%beg": -15,
+            "bc_y%end": -15,
+            # Immersed boundaries — all circles come from the particle bed
+            "ib": "T",
+            "num_ibs": 0,
+            "viscous": "T",
+            # Collision model (soft-sphere, from 3D_mibm_sphere_head_on_collision)
+            "collision_model": 1,
+            "coefficient_of_restitution": 0.9,
+            "collision_time": collision_time,
+            "ib_coefficient_of_friction": 0.1,
+            # Particle bed: 20 free-floating circles placed randomly in region
+            "num_particle_beds": 1,
+            "particle_bed(1)%x_centroid": bed_x,
+            "particle_bed(1)%y_centroid": bed_y,
+            "particle_bed(1)%z_centroid": 0.0,
+            "particle_bed(1)%length_x": bed_lx,
+            "particle_bed(1)%length_y": bed_ly,
+            "particle_bed(1)%length_z": 0.0,
+            "particle_bed(1)%num_particles": 20,
+            "particle_bed(1)%radius": particle_radius,
+            "particle_bed(1)%mass": particle_mass,
+            "particle_bed(1)%min_spacing": particle_min_spacing,
+            "particle_bed(1)%moving_ibm": 2,
+            "particle_bed(1)%seed": 42,
+            # Output
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "E_wrt": "T",
+            "ib_state_wrt": "F",
+            "parallel_io": "T",
+            # IC Patch 1: post-shock region (left of wave front)
+            "patch_icpp(1)%geometry": 3,
+            "patch_icpp(1)%x_centroid": 0.5 * wave_front - 0.25 * domain_size,
+            "patch_icpp(1)%y_centroid": 0.0,
+            "patch_icpp(1)%length_x": 0.5 * domain_size + wave_front,
+            "patch_icpp(1)%length_y": domain_size,
+            "patch_icpp(1)%vel(1)": post_shock_speed,
+            "patch_icpp(1)%vel(2)": 0.0,
+            "patch_icpp(1)%pres": post_shock_pressure,
+            "patch_icpp(1)%alpha_rho(1)": post_shock_density,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            # IC Patch 2: pre-shock region (right of wave front)
+            "patch_icpp(2)%geometry": 3,
+            "patch_icpp(2)%x_centroid": 0.5 * wave_front + 0.25 * domain_size,
+            "patch_icpp(2)%y_centroid": 0.0,
+            "patch_icpp(2)%length_x": 0.5 * domain_size - wave_front,
+            "patch_icpp(2)%length_y": domain_size,
+            "patch_icpp(2)%vel(1)": pre_shock_speed,
+            "patch_icpp(2)%vel(2)": 0.0,
+            "patch_icpp(2)%pres": pre_shock_pressure,
+            "patch_icpp(2)%alpha_rho(1)": pre_shock_density,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            # Fluid properties: air
+            "fluid_pp(1)%gamma": 1.0 / (gam_a - 1.0),
+            "fluid_pp(1)%pi_inf": 0,
+            "fluid_pp(1)%Re(1)": 2500000,
+        }
+    )
+)

examples/2D_mibm_shock_cylinder/case.py

@@ -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

src/common/include/case.fpp

@@ -4,10 +4,8 @@
 
 ! For pre-process.
 #:def analytical()
-
 #:enddef
 
 ! For moving immersed boundaries in simulation
 #:def mib_analytical()
-
 #:enddef

src/common/m_constants.fpp

@@ -16,15 +16,19 @@ module m_constants
     real(wp), parameter  :: verysmall = 1.e-12_wp       !< Very small number
     !> Radius cutoff to avoid division by zero for 3D spherical harmonic patch (geometry 14)
     real(wp), parameter :: small_radius = 1.e-32_wp
-    integer, parameter  :: num_stcls_min = 5                   !< Minimum # of stencils
-    integer, parameter  :: path_len = 400                      !< Maximum path length
-    integer, parameter  :: name_len = 50                       !< Maximum name length
-    integer, parameter  :: dflt_int = -100                     !< Default integer value
-    integer, parameter  :: fourier_rings = 5                   !< Fourier filter ring limit
-    integer, parameter  :: num_fluids_max = 10                 !< Maximum number of fluids in the simulation
-    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_stcls_min = 5             !< Minimum # of stencils
+    integer, parameter  :: path_len = 400                !< Maximum path length
+    integer, parameter  :: name_len = 50                 !< Maximum name length
+    integer, parameter  :: dflt_int = -100               !< Default integer value
+    integer, parameter  :: fourier_rings = 5             !< Fourier filter ring limit
+    integer, parameter  :: num_fluids_max = 10           !< Maximum number of fluids in the simulation
+    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 = 2050000  !< Maximum number of immersed boundary patches (patch_ib)
+    !> Max patches readable from the namelist; patch_ib grows to num_ib_patches_max only when particle beds are used
+    integer, parameter  :: num_ib_patches_max_namelist = 50000
+    integer, parameter  :: num_local_ibs_max = 2000            !< Maximum number of immersed boundary patches (patch_ib)
+    integer, parameter  :: num_particle_beds_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)
     integer, parameter  :: max_sph_harm_degree = 5             !< Max degree L for 3D spherical harmonic patch (geometry 14)

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

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)
 
+        @:ALLOCATE(models(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)

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.