Debug wall collisions and particle state writing (#1745)

Co-authored-by: Spencer Bryngelson <sbryngelson@gmail.com>

Author: danieljvickers

src/post_process/m_data_output.fpp

@@ -409,6 +409,7 @@ contains
         integer, dimension(num_procs)                   :: meshtypes
         integer                                         :: i
         integer                                         :: ierr
+        integer                                         :: extents_size
 
         if (format == format_silo) then
             ! For multidimensional data sets, the spatial extents of all of the grid(s) handled by the local processor(s) are
@@ -441,7 +442,8 @@ contains
 
                 err = DBSET2DSTRLEN(len(meshnames(1)))
                 err = DBMKOPTLIST(2, out%optlist)
-                err = DBADDIOPT(out%optlist, DBOPT_EXTENTS_SIZE, size(out%spatial_extents, 1))
+                extents_size = size(out%spatial_extents, 1)
+                err = DBADDIOPT(out%optlist, DBOPT_EXTENTS_SIZE, extents_size)
                 err = DBADDDOPT(out%optlist, DBOPT_EXTENTS, out%spatial_extents)
                 err = DBPUTMMESH(out%dbroot, 'rectilinear_grid', 16, num_procs, meshnames, len_trim(meshnames), meshtypes, &
                                  & out%optlist, ierr)
@@ -549,6 +551,7 @@ contains
         integer, dimension(num_procs)                   :: vartypes
         integer                                         :: i, j, k
         integer                                         :: ierr
+        integer                                         :: extents_size
 
         if (format == format_silo) then
             ! Determining the extents of the flow variable on each local process and gathering all this information on root process
@@ -567,7 +570,8 @@ contains
 
                 err = DBSET2DSTRLEN(len(varnames(1)))
                 err = DBMKOPTLIST(2, out%optlist)
-                err = DBADDIOPT(out%optlist, DBOPT_EXTENTS_SIZE, 2)
+                extents_size = size(out%data_extents, 1)
+                err = DBADDIOPT(out%optlist, DBOPT_EXTENTS_SIZE, extents_size)
                 err = DBADDDOPT(out%optlist, DBOPT_EXTENTS, out%data_extents)
                 err = DBPUTMVAR(out%dbroot, trim(varname), len_trim(varname), num_procs, varnames, len_trim(varnames), vartypes, &
                                 & out%optlist, ierr)

src/post_process/m_global_parameters.fpp

@@ -219,6 +219,23 @@ contains
             fluid_pp(i)%mu_bulk = dflt_real
         end do
 
+        num_particle_clouds = 0
+        do i = 1, num_particle_clouds_max
+            particle_cloud(i)%x_centroid = 0._wp
+            particle_cloud(i)%y_centroid = 0._wp
+            particle_cloud(i)%z_centroid = 0._wp
+            particle_cloud(i)%length_x = dflt_real
+            particle_cloud(i)%length_y = dflt_real
+            particle_cloud(i)%length_z = dflt_real
+            particle_cloud(i)%num_particles = 0
+            particle_cloud(i)%radius = dflt_real
+            particle_cloud(i)%mass = dflt_real
+            particle_cloud(i)%min_spacing = 0._wp
+            particle_cloud(i)%moving_ibm = 0
+            particle_cloud(i)%seed = 0
+            particle_cloud(i)%packing_method = dflt_int
+        end do
+
         ! Subgrid bubble parameters (bub_pp struct + scalar companions; bub_pp%R0ref is set in common
         ! via R0ref; the scalar companions are per-target manual declarations)
         bub_pp%R0ref = dflt_real; R0ref = dflt_real
@@ -323,6 +340,13 @@ contains
         ! defaulted on non-root ranks nor broadcast. Post-process never carries viscous stresses.
         viscous = .false.
 
+        ! Particle clouds expand into individual IB patches at simulation startup, so num_ibs as read
+        ! from the case file counts only the namelist patches. Match the global count the simulation
+        ! arrives at (s_reduce_ib_patch_array) so the IB state records can be read back.
+        do i = 1, num_particle_clouds
+            num_ibs = num_ibs + particle_cloud(i)%num_particles
+        end do
+
         ! Gamma/Pi_inf: force num_fluids=1 (post_process-specific side effect of the gamma-law model)
         if (model_eqns == model_eqns_gamma_law) num_fluids = 1
 

src/post_process/m_mpi_proxy.fpp

@@ -93,6 +93,18 @@ contains
             #:endfor
         #:endfor
 
+        ! manual: particle_cloud (runtime loop to num_particle_clouds; irregular member subset)
+        do i = 1, num_particle_clouds
+            #:for VAR in ['x_centroid', 'y_centroid', 'z_centroid', 'length_x', 'length_y', 'length_z', &
+                & 'radius', 'mass', 'min_spacing']
+                call MPI_BCAST(particle_cloud(i)%${VAR}$, 1, mpi_p, 0, MPI_COMM_WORLD, ierr)
+            #:endfor
+            call MPI_BCAST(particle_cloud(i)%num_particles, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%moving_ibm, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%seed, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+            call MPI_BCAST(particle_cloud(i)%packing_method, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+        end do
+
         ! manual: cfl_dt (runtime-computed logical), bc_io (BC-file existence)
         call MPI_BCAST(cfl_dt, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
         call MPI_BCAST(bc_io, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)

src/simulation/m_collisions.fpp

@@ -340,8 +340,6 @@ contains
         integer  :: gp_idx, i, j, k, patch_id
         real(wp) :: edge_location, overlap_distance
 
-        ! iterate over all ghost points to detect the one that is most-overlapping in each direction
-
         $:GPU_PARALLEL_LOOP(private='[patch_id, edge_location, overlap_distance]')
         do patch_id = 1, num_ibs
             #:for X, DIR, IDX in [('x', 1, 1), ('y', 2, 3), ('z', 3, 5)]

src/simulation/m_global_parameters.fpp

@@ -121,7 +121,7 @@ module m_global_parameters
     $:GPU_DECLARE(create='[bc_x%vb1, bc_x%vb2, bc_x%vb3, bc_x%ve1, bc_x%ve2, bc_x%ve3]')
     $:GPU_DECLARE(create='[bc_y%vb1, bc_y%vb2, bc_y%vb3, bc_y%ve1, bc_y%ve2, bc_y%ve3]')
     $:GPU_DECLARE(create='[bc_z%vb1, bc_z%vb2, bc_z%vb3, bc_z%ve1, bc_z%ve2, bc_z%ve3]')
-    $:GPU_DECLARE(create='[ib_bc_x%beg, ib_bc_y%beg, ib_bc_z%beg]')
+    $:GPU_DECLARE(create='[ib_bc_x%beg, ib_bc_x%end, ib_bc_y%beg, ib_bc_y%end, ib_bc_z%beg, ib_bc_z%end]')
 #elif defined(MFC_OpenMP)
     $:GPU_DECLARE(create='[bc_x, bc_y, bc_z]')
     $:GPU_DECLARE(create='[ib_bc_x, ib_bc_y, ib_bc_z]')

src/simulation/m_ibm.fpp

@@ -76,9 +76,9 @@ contains
         call nvtxStartRange("SETUP-IBM-MODULE")
 
         ! GPU routines require updated cell centers
-        $:GPU_UPDATE(device='[num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, ib_bc_x%beg, ib_bc_y%beg]')
+        $:GPU_UPDATE(device='[num_ibs, num_gbl_ibs, x_cc, y_cc, dx, dy, ib_bc_x%beg, ib_bc_x%end, ib_bc_y%beg, ib_bc_y%end]')
         if (p /= 0) then
-            $:GPU_UPDATE(device='[z_cc, dz, ib_bc_z%beg]')
+            $:GPU_UPDATE(device='[z_cc, dz, ib_bc_z%beg, ib_bc_z%end]')
         end if
         $:GPU_UPDATE(device='[patch_ib(1:num_ibs), glb_bounds]')
 

toolchain/mfc/params/definitions.py

@@ -1314,6 +1314,8 @@ def _decl(targets: set, *names: str) -> None:
     "avg_state",
     "alt_soundspeed",
     "mixture_err",
+    "num_particle_clouds",
+    "particle_cloud",
 )
 _nv(
     _PRE_SIM,
@@ -1395,10 +1397,8 @@ def _decl(targets: set, *names: str) -> None:
     "coefficient_of_restitution",
     "collision_time",
     "ib_coefficient_of_friction",
-    "num_particle_clouds",
     "ib_neighborhood_radius",
     "many_ib_patch_parallelism",
-    "particle_cloud",
     "tau_star",
     "cont_damage_s",
     "alpha_bar",