Add hemisphere shell particle cloud packing

docs/documentation/case.md

@@ -425,18 +425,25 @@ A particle cloud is a compact specification of a bed of identical circular (2D)
 | Parameter         | Type    | Description |
 | ---:              | :----:  | :---        |
 | `x[y,z]_centroid` | Real    | Centre of the cloud region in the [x,y,z]-direction. |
-| `length_x[y,z]`   | Real    | Extent of the cloud region in the [x,y,z]-direction. |
+| `length_x[y,z]`   | Real    | Extent of the cloud region in the [x,y,z]-direction for `geometry = 1`; ignored by `geometry = 2`. |
 | `num_particles`   | Integer | Number of particles to place in the region. |
 | `radius`          | Real    | Radius of every particle in the cloud. |
 | `mass`            | Real    | Mass of every particle in the cloud. |
 | `min_spacing`     | Real    | Minimum surface-to-surface gap between particles (centres are `2*radius + min_spacing` apart). |
+| `geometry`        | Integer | Shape of the cloud region. |
+| `shell_inner_radius` | Real | Inner radius for hemisphere-shell clouds (`geometry = 2`). |
+| `shell_outer_radius` | Real | Outer radius for hemisphere-shell clouds (`geometry = 2`). |
 | `moving_ibm`      | Integer | Motion flag applied to every particle (see `patch_ib(j)%%moving_ibm`). |
 | `seed`            | Integer | Random seed for reproducible placement (used by `packing_method = 1`). |
 | `packing_method`  | Integer | Algorithm used to place the particles. |
 
+- `geometry` selects the cloud region:
+  - `1` (box) uses `x[y,z]_centroid` and `length_x[y,z]` to define the region.
+  - `2` (hemisphere shell) uses `x[y,z]_centroid`, `shell_inner_radius`, and `shell_outer_radius` to define the region. Particle centres are sampled between `shell_inner_radius + radius` and `shell_outer_radius - radius`, and the flat hemisphere plane is kept clear by one particle radius.
 - `packing_method` selects how the `num_particles` are positioned within the cloud region:
   - `1` (rejection sampling) draws random positions and rejects any that violate `min_spacing`, producing a disordered bed. `seed` makes the placement reproducible.
   - `2` (lattice) places the particles on the optimally dense lattice for the geometry — a triangular lattice in 2D and a face-centered cubic lattice in 3D. The lattice spacing is derived from the particle density (`num_particles` over the region area/volume); if that spacing is below the required `2*radius + min_spacing`, the region is too dense and the run aborts.
+  - Hemisphere-shell clouds currently support rejection sampling only; `geometry = 2` with `packing_method = 2` is rejected during input validation.
 
 ### 5. Fluid Material's {#sec-fluid-materials}
 

src/common/m_derived_types.fpp

@@ -374,8 +374,11 @@ module m_derived_types
         real(wp) :: radius  !< Particle radius
         real(wp) :: mass  !< Particle mass
         real(wp) :: min_spacing  !< Minimum surface-to-surface gap (particle centers are 2*radius + min_spacing apart)
+        real(wp) :: shell_inner_radius  !< Inner radius for shell packing
+        real(wp) :: shell_outer_radius  !< Outer radius for shell packing
         integer  :: moving_ibm  !< Motion flag: 0=static, 1=moving (forces), 2=forced path
         integer  :: seed  !< Random seed for reproducible placement
+        integer  :: geometry  !< Cloud geometry: 1=box, 2=hemisphere shell
         integer  :: packing_method  !< Packing algorithm: 1=rejection sampling, 2=lattice
     end type particle_cloud_parameters
 

src/simulation/m_checker.fpp

@@ -183,19 +183,30 @@ contains
 
     end subroutine s_check_inputs_ib_injection
 
-    !> Checks that each active particle cloud has a valid packing_method specified
+    !> Checks that each active particle cloud has a valid geometry and packing_method specified
     impure subroutine s_check_inputs_particle_clouds
 
         integer          :: i
         character(len=5) :: idxStr
 
         do i = 1, num_particle_clouds
             call s_int_to_str(i, idxStr)
+            @:PROHIBIT(particle_cloud(i)%geometry /= 1 .and. particle_cloud(i)%geometry /= 2, &
+                       & "particle_cloud("//trim(idxStr) //")%geometry must be 1 (box) or 2 (hemisphere shell)")
             @:PROHIBIT(particle_cloud(i)%packing_method == dflt_int, &
                        & "particle_cloud("//trim(idxStr) &
                        & //")%packing_method must be specified (1 = rejection sampling, 2 = lattice)")
             @:PROHIBIT(particle_cloud(i)%packing_method /= 1 .and. particle_cloud(i)%packing_method /= 2, &
                        & "particle_cloud("//trim(idxStr) //")%packing_method must be 1 (rejection sampling) or 2 (lattice)")
+            @:PROHIBIT(particle_cloud(i)%geometry == 2 .and. particle_cloud(i)%shell_inner_radius < 0._wp, &
+                       & "particle_cloud("//trim(idxStr) //") hemisphere shell requires shell_inner_radius >= 0")
+            @:PROHIBIT(particle_cloud(i)%geometry == 2 &
+                       & .and. particle_cloud(i)%shell_outer_radius <= particle_cloud(i)%shell_inner_radius &
+                       & + 2._wp*particle_cloud(i)%radius, &
+                       & "particle_cloud("//trim(idxStr) &
+                       & //") hemisphere shell requires shell_outer_radius > shell_inner_radius + 2*radius")
+            @:PROHIBIT(particle_cloud(i)%geometry == 2 .and. particle_cloud(i)%packing_method == 2, &
+                       & "particle_cloud("//trim(idxStr) //") hemisphere-shell lattice packing is not implemented")
         end do
 
     end subroutine s_check_inputs_particle_clouds

src/simulation/m_global_parameters.fpp

@@ -643,8 +643,11 @@ contains
             particle_cloud(i)%radius = dflt_real
             particle_cloud(i)%mass = dflt_real
             particle_cloud(i)%min_spacing = 0._wp
+            particle_cloud(i)%shell_inner_radius = dflt_real
+            particle_cloud(i)%shell_outer_radius = dflt_real
             particle_cloud(i)%moving_ibm = 0
             particle_cloud(i)%seed = 0
+            particle_cloud(i)%geometry = 1
             particle_cloud(i)%packing_method = dflt_int
         end do
 

src/simulation/m_mpi_proxy.fpp

@@ -200,12 +200,13 @@ contains
         ! 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']
+                & 'radius', 'mass', 'min_spacing', 'shell_inner_radius', 'shell_outer_radius']
                 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)%geometry, 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
 

src/simulation/m_particle_cloud.fpp

@@ -43,11 +43,21 @@ contains
         ib_idx = 0  ! index into particle_cloud_ibs
 
         do cloud_idx = 1, num_particle_clouds
-            select case (particle_cloud(cloud_idx)%packing_method)
-            case (1)  ! random box packing method
-                call s_particle_cloud_random_box(cloud_idx, ib_idx, particle_cloud_ibs)
-            case (2)  ! lattice packing method
-                call s_particle_cloud_lattice(cloud_idx, ib_idx, particle_cloud_ibs)
+            select case (particle_cloud(cloud_idx)%geometry)
+            case (1)  ! box geometry
+                select case (particle_cloud(cloud_idx)%packing_method)
+                case (1)  ! rejection sampling method
+                    call s_particle_cloud_random_box(cloud_idx, ib_idx, particle_cloud_ibs)
+                case (2)  ! lattice packing method
+                    call s_particle_cloud_lattice(cloud_idx, ib_idx, particle_cloud_ibs)
+                end select
+            case (2)  ! hemisphere-shell geometry
+                select case (particle_cloud(cloud_idx)%packing_method)
+                case (1)  ! rejection sampling method
+                    call s_particle_cloud_random_hemi_shell(cloud_idx, ib_idx, particle_cloud_ibs)
+                case (2)
+                    call s_mpi_abort("Error :: Lattice packing is not implemented for hemisphere-shell particle clouds")
+                end select
             end select
         end do
 
@@ -65,12 +75,10 @@ contains
         integer                                                :: n_placed, geom, seed
         integer(8)                                             :: n_attempts, max_attempts
         real(wp)                                               :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist
-        real(wp)                                               :: rx, ry, rz, dist
-        logical                                                :: overlaps
+        real(wp)                                               :: rx, ry, rz
         real(wp), allocatable                                  :: placed(:,:)
         integer                                                :: hash_size, slot
-        integer                                                :: bx, by, bz, nbx, nby, nbz
-        integer                                                :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j
+        integer                                                :: bx, by, bz
         integer, allocatable                                   :: hash_head(:), chain_next(:)
 
         xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x
@@ -82,14 +90,10 @@ contains
 
         min_dist = 2._wp*particle_cloud(cloud_idx)%radius + particle_cloud(cloud_idx)%min_spacing
 
-        if (p == 0) then
+        if (num_dims < 3) then
             geom = 2  ! circle for 2D
-            dz_lo = 0
-            dz_hi = 0
         else
             geom = 8  ! sphere for 3D
-            dz_lo = -1
-            dz_hi = 1
         end if
 
         max_attempts = int(particle_cloud(cloud_idx)%num_particles, 8)*1000_8
@@ -113,50 +117,20 @@ contains
 
             rx = xmin + f_xorshift(seed)*(xmax - xmin)
             ry = ymin + f_xorshift(seed)*(ymax - ymin)
-            if (p == 0) then
+            if (num_dims < 3) then
                 rz = particle_cloud(cloud_idx)%z_centroid
             else
                 rz = zmin + f_xorshift(seed)*(zmax - zmin)
             end if
 
-            bx = int(floor(rx/min_dist))
-            by = int(floor(ry/min_dist))
-            bz = 0
-            if (p /= 0) bz = int(floor(rz/min_dist))
-
-            ! Check 3x3(x3) neighboring bins - O(1) average via hash lookup
-            overlaps = .false.
-            outer: do dx_b = -1, 1
-                do dy_b = -1, 1
-                    do dz_b = dz_lo, dz_hi
-                        nbx = bx + dx_b
-                        nby = by + dy_b
-                        nbz = bz + dz_b
-                        slot = f_bin_hash(nbx, nby, nbz, hash_size)
-                        j = hash_head(slot)
-                        do while (j > 0)
-                            if (p == 0) then
-                                dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2)
-                            else
-                                dist = sqrt((rx - placed(1, j))**2 + (ry - placed(2, j))**2 + (rz - placed(3, j))**2)
-                            end if
-                            if (dist < min_dist) then
-                                overlaps = .true.
-                                exit outer
-                            end if
-                            j = chain_next(j)
-                        end do
-                    end do
-                end do
-            end do outer
-
-            if (.not. overlaps) then
+            if (.not. f_cloud_particle_overlaps(rx, ry, rz, placed, hash_head, chain_next, hash_size, min_dist)) then
                 n_placed = n_placed + 1
                 placed(1, n_placed) = rx
                 placed(2, n_placed) = ry
                 placed(3, n_placed) = rz
 
                 ! Insert into hash grid as head of bucket chain
+                call s_get_cloud_bin(rx, ry, rz, min_dist, bx, by, bz)
                 slot = f_bin_hash(bx, by, bz, hash_size)
                 chain_next(n_placed) = hash_head(slot)
                 hash_head(slot) = n_placed
@@ -173,6 +147,100 @@ contains
 
     end subroutine s_particle_cloud_random_box
 
+    !> Generates a random distribution of particles in a hemisphere shell with a minimum spacing.
+    subroutine s_particle_cloud_random_hemi_shell(cloud_idx, ib_idx, particle_cloud_ibs)
+
+        integer, intent(in)                                    :: cloud_idx
+        integer, intent(inout)                                 :: ib_idx
+        type(ib_patch_parameters), intent(inout), dimension(:) :: particle_cloud_ibs
+        integer                                                :: n_placed, geom, seed
+        integer(8)                                             :: n_attempts, max_attempts
+        real(wp)                                               :: min_dist
+        real(wp)                                               :: rx, ry, rz, theta, phi, r_shell, rho, u
+        real(wp)                                               :: r_inner, r_outer, zdir
+        real(wp), allocatable                                  :: placed(:,:)
+        integer                                                :: hash_size, slot
+        integer                                                :: bx, by, bz
+        integer, allocatable                                   :: hash_head(:), chain_next(:)
+
+        r_inner = particle_cloud(cloud_idx)%shell_inner_radius + particle_cloud(cloud_idx)%radius
+        r_outer = particle_cloud(cloud_idx)%shell_outer_radius - particle_cloud(cloud_idx)%radius
+        min_dist = 2._wp*particle_cloud(cloud_idx)%radius + particle_cloud(cloud_idx)%min_spacing
+
+        if (r_inner < 0._wp .or. r_outer <= r_inner) then
+            call s_mpi_abort("Error :: Invalid hemisphere-shell radii for particle cloud packing")
+        end if
+
+        if (num_dims < 3) then
+            geom = 2
+        else
+            geom = 8
+        end if
+
+        max_attempts = int(particle_cloud(cloud_idx)%num_particles, 8)*1000_8
+        n_placed = 0
+        n_attempts = 0
+        seed = particle_cloud(cloud_idx)%seed
+        if (seed == 0) seed = 1 + cloud_idx*1013904223
+
+        allocate (placed(3, particle_cloud(cloud_idx)%num_particles))
+
+        hash_size = max(16, 4*particle_cloud(cloud_idx)%num_particles)
+        allocate (hash_head(hash_size))
+        allocate (chain_next(particle_cloud(cloud_idx)%num_particles))
+        hash_head = -1
+        chain_next = -1
+
+        do while (n_placed < particle_cloud(cloud_idx)%num_particles .and. n_attempts < max_attempts)
+            n_attempts = n_attempts + 1
+
+            if (num_dims < 3) then
+                theta = pi*f_xorshift(seed)
+                u = f_xorshift(seed)
+                r_shell = sqrt((r_outer**2._wp - r_inner**2._wp)*u + r_inner**2._wp)
+                rx = particle_cloud(cloud_idx)%x_centroid + r_shell*cos(theta)
+                ry = particle_cloud(cloud_idx)%y_centroid + r_shell*sin(theta)
+                rz = particle_cloud(cloud_idx)%z_centroid
+            else
+                phi = 2._wp*pi*f_xorshift(seed)
+                zdir = f_xorshift(seed)
+                rho = sqrt(max(0._wp, 1._wp - zdir**2._wp))
+                u = f_xorshift(seed)
+                r_shell = ((r_outer**3._wp - r_inner**3._wp)*u + r_inner**3._wp)**(1._wp/3._wp)
+                rx = particle_cloud(cloud_idx)%x_centroid + r_shell*rho*cos(phi)
+                ry = particle_cloud(cloud_idx)%y_centroid + r_shell*rho*sin(phi)
+                rz = particle_cloud(cloud_idx)%z_centroid + r_shell*zdir
+            end if
+
+            if (num_dims < 3) then
+                if (ry < particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius) cycle
+            else
+                if (rz < particle_cloud(cloud_idx)%z_centroid + particle_cloud(cloud_idx)%radius) cycle
+            end if
+
+            if (.not. f_cloud_particle_overlaps(rx, ry, rz, placed, hash_head, chain_next, hash_size, min_dist)) then
+                n_placed = n_placed + 1
+                placed(1, n_placed) = rx
+                placed(2, n_placed) = ry
+                placed(3, n_placed) = rz
+
+                call s_get_cloud_bin(rx, ry, rz, min_dist, bx, by, bz)
+                slot = f_bin_hash(bx, by, bz, hash_size)
+                chain_next(n_placed) = hash_head(slot)
+                hash_head(slot) = n_placed
+
+                call s_add_cloud_particle(cloud_idx, ib_idx, geom, rx, ry, rz, particle_cloud_ibs)
+            end if
+        end do
+
+        if (n_placed < particle_cloud(cloud_idx)%num_particles) then
+            call s_mpi_abort("Error :: Failed to place all particles in hemisphere-shell particle bed")
+        end if
+
+        deallocate (placed, hash_head, chain_next)
+
+    end subroutine s_particle_cloud_random_hemi_shell
+
     !> Places particles on the optimally dense lattice for the cloud region: a triangular lattice in 2D, a face-centered cubic
     !! lattice in 3D. The lattice spacing is set by the particle density (num_particles over the region area/volume); if that
     !! spacing falls below the required centre-to-centre distance (2*radius + min_spacing), the region is too dense and the run is
@@ -313,6 +381,69 @@ contains
 
     end subroutine s_add_cloud_particle
 
+    !> Convert a candidate particle centre to spatial-hash bin coordinates.
+    subroutine s_get_cloud_bin(px, py, pz, min_dist, bx, by, bz)
+
+        real(wp), intent(in) :: px, py, pz, min_dist
+        integer, intent(out) :: bx, by, bz
+
+        bx = int(floor(px/min_dist))
+        by = int(floor(py/min_dist))
+        bz = 0
+        if (num_dims == 3) bz = int(floor(pz/min_dist))
+
+    end subroutine s_get_cloud_bin
+
+    !> Check whether a candidate particle centre overlaps any already-placed particle in neighbouring spatial-hash bins.
+    function f_cloud_particle_overlaps(px, py, pz, placed, hash_head, chain_next, hash_size, min_dist) result(overlaps)
+
+        real(wp), intent(in)                 :: px, py, pz, min_dist
+        real(wp), intent(in), dimension(:,:) :: placed
+        integer, intent(in), dimension(:)    :: hash_head, chain_next
+        integer, intent(in)                  :: hash_size
+        logical                              :: overlaps
+        integer                              :: bx, by, bz, nbx, nby, nbz, slot
+        integer                              :: dx_b, dy_b, dz_b, dz_lo, dz_hi, j
+        real(wp)                             :: dist_sq, min_dist_sq
+
+        call s_get_cloud_bin(px, py, pz, min_dist, bx, by, bz)
+
+        dz_lo = -1
+        dz_hi = 1
+        if (num_dims < 3) then
+            dz_lo = 0
+            dz_hi = 0
+        end if
+
+        min_dist_sq = min_dist**2._wp
+        overlaps = .false.
+
+        do dx_b = -1, 1
+            do dy_b = -1, 1
+                do dz_b = dz_lo, dz_hi
+                    nbx = bx + dx_b
+                    nby = by + dy_b
+                    nbz = bz + dz_b
+                    slot = f_bin_hash(nbx, nby, nbz, hash_size)
+                    j = hash_head(slot)
+                    do while (j > 0)
+                        if (num_dims < 3) then
+                            dist_sq = (px - placed(1, j))**2._wp + (py - placed(2, j))**2._wp
+                        else
+                            dist_sq = (px - placed(1, j))**2._wp + (py - placed(2, j))**2._wp + (pz - placed(3, j))**2._wp
+                        end if
+                        if (dist_sq < min_dist_sq) then
+                            overlaps = .true.
+                            return
+                        end if
+                        j = chain_next(j)
+                    end do
+                end do
+            end do
+        end do
+
+    end function f_cloud_particle_overlaps
+
     !> Xorshift PRNG. Advances seed in-place and returns a value in [0, 1).
     function f_xorshift(seed) result(rval)