Optimized particle bed instantitation

Author: Daniel Vickers

src/common/m_model.fpp

@@ -983,7 +983,7 @@ contains
         dx_local = minval(dx); dy_local = minval(dy)
         if (p /= 0) dz_local = minval(dz)
 
-        num_gbl_ibs = num_ibs
+        @:ALLOCATE(models(num_ibs))
         allocate (stl_bounding_boxes(num_ibs,1:3,1:3))
 
         do patch_id = 1, num_ibs

src/simulation/m_collisions.fpp

@@ -39,8 +39,8 @@ contains
         spring_stiffness = (pi**2 + log(e)**2)/(collision_time**2)
         $:GPU_UPDATE(device='[damping_parameter, spring_stiffness]')
 
-        @:ALLOCATE(collision_lookup(num_ibs * 8, 4))
-        @:ALLOCATE(wall_overlap_distances(num_ibs, 6))
+        @:ALLOCATE(collision_lookup(num_local_ibs_max * 27 * 8, 4))
+        @:ALLOCATE(wall_overlap_distances(num_local_ibs_max*27, 6))
 
         wall_overlap_distances = 0
         $:GPU_UPDATE(device='[wall_overlap_distances]')

src/simulation/m_ibm.fpp

@@ -60,8 +60,6 @@ contains
             @:ALLOCATE(ib_markers%sf(-buff_size:m+buff_size, -buff_size:n+buff_size, 0:0))
         end if
 
-        @:ALLOCATE(models(num_ibs))
-
         @:ACC_SETUP_SFs(ib_markers)
 
         $:GPU_ENTER_DATA(copyin='[num_gps]')
@@ -1034,6 +1032,7 @@ contains
         end do
         $:END_GPU_PARALLEL_LOOP()
 
+        if (proc_rank == 0) print *, "s_apply_collision_forces"
         call s_apply_collision_forces(ghost_points, num_gps, ib_markers, forces, torques)
 
         ! reduce the forces across local neighborhood ranks

src/simulation/m_particle_bed.fpp

@@ -22,10 +22,11 @@ contains
     impure subroutine s_generate_particle_beds()
 
         integer               :: b, ib_idx, geom
-        integer               :: n_placed
+        integer               :: n_placed, n_total_placed
         integer(8)            :: n_attempts, max_attempts
         real(wp)              :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist
         real(wp)              :: rx, ry, rz, dist
+        real(wp)              :: t_start, t_end
         integer               :: seed
         logical               :: overlaps
         real(wp), allocatable :: placed(:,:)
@@ -38,6 +39,9 @@ contains
 
         if (num_particle_beds == 0) return
 
+        call cpu_time(t_start)
+        n_total_placed = 0
+
         do b = 1, num_particle_beds
             xmin = particle_bed(b)%x_centroid - 0.5_wp*particle_bed(b)%length_x
             xmax = particle_bed(b)%x_centroid + 0.5_wp*particle_bed(b)%length_x
@@ -58,7 +62,7 @@ contains
                 dz_hi = 1
             end if
 
-            max_attempts = int(particle_bed(b)%num_particles, 8)*10000_8
+            max_attempts = int(particle_bed(b)%num_particles, 8)*1000_8
             n_placed = 0
             n_attempts = 0
             seed = particle_bed(b)%seed
@@ -135,20 +139,33 @@ contains
                     patch_ib(ib_idx)%x_centroid = rx
                     patch_ib(ib_idx)%y_centroid = ry
                     patch_ib(ib_idx)%z_centroid = rz
+                    patch_ib(ib_idx)%angles(1) = 0._wp
+                    patch_ib(ib_idx)%angles(2) = 0._wp
+                    patch_ib(ib_idx)%angles(3) = 0._wp
+                    patch_ib(ib_idx)%vel(1) = 0._wp
+                    patch_ib(ib_idx)%vel(2) = 0._wp
+                    patch_ib(ib_idx)%vel(3) = 0._wp
+                    patch_ib(ib_idx)%angular_vel(1) = 0._wp
+                    patch_ib(ib_idx)%angular_vel(2) = 0._wp
+                    patch_ib(ib_idx)%angular_vel(3) = 0._wp
                     patch_ib(ib_idx)%radius = particle_bed(b)%radius
                     patch_ib(ib_idx)%mass = particle_bed(b)%mass
                     patch_ib(ib_idx)%moving_ibm = particle_bed(b)%moving_ibm
                 end if
             end do
 
             if (n_placed < particle_bed(b)%num_particles) then
-                print '("WARNING: particle_bed(",I0,"): placed ",I0," of ",I0," particles after ",I0," attempts")', &
-                    b, n_placed, particle_bed(b)%num_particles, n_attempts
+              print *, "Error :: Failed to place all IBs ib particle bed"
+              stop
             end if
 
+            n_total_placed = n_total_placed + n_placed
             deallocate (placed, hash_head, chain_next)
         end do
 
+        call cpu_time(t_end)
+        if (proc_rank == 0) print '(a,i0,a,f0.3,a)', 'Particle beds placed ', n_total_placed, ' particles in ', t_end - t_start, ' seconds.'
+
     end subroutine s_generate_particle_beds
 
     !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1).

src/simulation/m_start_up.fpp

@@ -922,6 +922,9 @@ contains
                 call s_read_ib_restart_data(n_start)
             else if (t_step_start > 0) then
                 call s_read_ib_restart_data(t_step_start)
+            else
+                ! particle bed generated on first tiem step
+                call s_generate_particle_beds()
             end if
             call s_instantiate_STL_models()
             call s_reduce_ib_patch_array()
@@ -1005,7 +1008,6 @@ contains
             call s_assign_default_values_to_user_inputs()
             call s_read_input_file()
             call s_check_input_file()
-            call s_generate_particle_beds()
 
             print '(" Simulating a ", A, " ", I0, "x", I0, "x", I0, " case on ", I0, " rank(s) ", A, ".")', &
             #:if not MFC_CASE_OPTIMIZATION
@@ -1217,7 +1219,6 @@ contains
         logical                                                  :: is_in_neighborhood, is_local
 
         ! do all set up for moving immersed boundaries
-
         moving_immersed_boundary_flag = .false.
         do i = 1, num_ibs
             if (patch_ib(i)%moving_ibm /= 0) then
@@ -1289,6 +1290,8 @@ contains
         integer, dimension(4) :: buf4, rbuf4
         integer, dimension(2) :: buf2, rbuf2
 
+        if (proc_rank == 0) print *, "Entering compute_ib_neighbor_ranks"
+
         ax = ib_neighborhood_radius
 
         if (allocated(ib_neighbor_ranks)) deallocate (ib_neighbor_ranks)
@@ -1432,6 +1435,8 @@ contains
         end if
 #endif
 
+        if (proc_rank == 0) print *, "Exiting compute_ib_neighbor_ranks"
+
     end subroutine s_compute_ib_neighbor_ranks
 
     subroutine get_neighbor_bounds()