From 81e9038632be28d403678d26f0710975ce57fcf6 Mon Sep 17 00:00:00 2001 From: BC Kim Date: Wed, 22 Jul 2026 01:24:39 -0400 Subject: [PATCH] Add hemisphere shell particle cloud packing --- .../3D_mibm_particle_cloud_hemi_shell/case.py | 101 ++++++++++++ src/common/m_derived_types.fpp | 4 +- src/simulation/m_checker.fpp | 11 +- src/simulation/m_global_parameters.fpp | 2 + src/simulation/m_mpi_proxy.fpp | 2 +- src/simulation/m_particle_cloud.fpp | 145 ++++++++++++++++++ toolchain/mfc/case_validator.py | 6 +- toolchain/mfc/params/definitions.py | 9 +- 8 files changed, 267 insertions(+), 13 deletions(-) create mode 100644 examples/3D_mibm_particle_cloud_hemi_shell/case.py diff --git a/examples/3D_mibm_particle_cloud_hemi_shell/case.py b/examples/3D_mibm_particle_cloud_hemi_shell/case.py new file mode 100644 index 0000000000..2de4f4e4d7 --- /dev/null +++ b/examples/3D_mibm_particle_cloud_hemi_shell/case.py @@ -0,0 +1,101 @@ +import json + +gam_a = 1.4 + +domain_size = 4.0 +particle_radius = 0.08 + +print( + json.dumps( + { + # Logistics + "run_time_info": "T", + # Computational Domain Parameters + "x_domain%beg": -0.5 * domain_size, + "x_domain%end": 0.5 * domain_size, + "y_domain%beg": -0.5 * domain_size, + "y_domain%end": 0.5 * domain_size, + "z_domain%beg": 0.0, + "z_domain%end": domain_size, + "cyl_coord": "F", + "m": 48, + "n": 48, + "p": 48, + "dt": 1.0e-4, + "t_step_start": 0, + "t_step_stop": 1, + "t_step_save": 1, + # Simulation Algorithm Parameters + "num_patches": 1, + "model_eqns": "5eq", + "alt_soundspeed": "F", + "num_fluids": 1, + "mpp_lim": "F", + "mixture_err": "T", + "time_stepper": "rk3", + "weno_order": 5, + "weno_eps": 1.0e-16, + "weno_avg": "T", + "avg_state": "arithmetic", + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "T", + "riemann_solver": "hllc", + "wave_speeds": "direct", + "fd_order": 2, + "bc_x%beg": -3, + "bc_x%end": -3, + "bc_y%beg": -3, + "bc_y%end": -3, + "bc_z%beg": -15, + "bc_z%end": -3, + # Immersed boundaries: all spheres come from the particle cloud + "ib": "T", + "num_ibs": 0, + "viscous": "T", + "many_ib_patch_parallelism": "T", + # Hemi-shell particle cloud + "num_particle_clouds": 1, + "particle_cloud(1)%x_centroid": 0.0, + "particle_cloud(1)%y_centroid": 0.0, + "particle_cloud(1)%z_centroid": 0.0, + "particle_cloud(1)%length_x": domain_size, + "particle_cloud(1)%length_y": domain_size, + "particle_cloud(1)%length_z": domain_size, + "particle_cloud(1)%num_particles": 40, + "particle_cloud(1)%radius": particle_radius, + "particle_cloud(1)%mass": 1.0, + "particle_cloud(1)%min_spacing": 0.02, + "particle_cloud(1)%shell_inner_radius": 0.75, + "particle_cloud(1)%shell_outer_radius": 1.5, + "particle_cloud(1)%moving_ibm": 0, + "particle_cloud(1)%seed": 42, + "particle_cloud(1)%packing_method": 3, + # Output + "format": "silo", + "precision": "double", + "prim_vars_wrt": "T", + "E_wrt": "T", + "ib_state_wrt": "T", + "parallel_io": "T", + # IC Patch 1: air + "patch_icpp(1)%geometry": 9, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%z_centroid": 0.5 * domain_size, + "patch_icpp(1)%length_x": domain_size, + "patch_icpp(1)%length_y": domain_size, + "patch_icpp(1)%length_z": domain_size, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%vel(3)": 0.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": 1.0, + "patch_icpp(1)%alpha(1)": 1.0, + # Fluid properties: air + "fluid_pp(1)%gamma": 1.0 / (gam_a - 1.0), + "fluid_pp(1)%pi_inf": 0.0, + "fluid_pp(1)%Re(1)": 2.5e6, + } + ) +) diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index 5aadd664a0..d412a1676c 100644 --- a/src/common/m_derived_types.fpp +++ b/src/common/m_derived_types.fpp @@ -371,9 +371,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 :: packing_method !< Packing algorithm: 1=rejection sampling, 2=lattice + integer :: packing_method !< Packing algorithm: 1=rejection sampling, 2=lattice, 3=hemisphere shell end type particle_cloud_parameters !> Derived type annexing the physical parameters (PP) of the fluids. These include the specific heat ratio function and liquid diff --git a/src/simulation/m_checker.fpp b/src/simulation/m_checker.fpp index 7ebe1eed27..15d1e7a6ea 100644 --- a/src/simulation/m_checker.fpp +++ b/src/simulation/m_checker.fpp @@ -142,9 +142,14 @@ contains call s_int_to_str(i, idxStr) @: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)") + & //")%packing_method must be specified (1 = rejection sampling, 2 = lattice, 3 = hemisphere shell)") + @:PROHIBIT(particle_cloud(i)%packing_method /= 1 .and. particle_cloud(i)%packing_method /= 2 & + & .and. particle_cloud(i)%packing_method /= 3, & + & "particle_cloud("//trim(idxStr) & + & //")%packing_method must be 1 (rejection sampling), 2 (lattice), or 3 (hemisphere shell)") + @:PROHIBIT(particle_cloud(i)%packing_method == 3 & + & .and. particle_cloud(i)%shell_outer_radius <= particle_cloud(i)%shell_inner_radius, & + & "particle_cloud("//trim(idxStr) //") hemisphere shell requires shell_outer_radius > shell_inner_radius") end do end subroutine s_check_inputs_particle_clouds diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index 02f7c72915..f69d0db7d3 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/src/simulation/m_global_parameters.fpp @@ -647,6 +647,8 @@ 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)%packing_method = dflt_int diff --git a/src/simulation/m_mpi_proxy.fpp b/src/simulation/m_mpi_proxy.fpp index 7300ef33dc..c2a23ed496 100644 --- a/src/simulation/m_mpi_proxy.fpp +++ b/src/simulation/m_mpi_proxy.fpp @@ -198,7 +198,7 @@ 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) diff --git a/src/simulation/m_particle_cloud.fpp b/src/simulation/m_particle_cloud.fpp index 3a771093d3..1fd21abdcb 100644 --- a/src/simulation/m_particle_cloud.fpp +++ b/src/simulation/m_particle_cloud.fpp @@ -48,6 +48,8 @@ contains 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) + case (3) ! random hemisphere-shell packing method + call s_particle_cloud_random_hemi_shell(cloud_idx, ib_idx, particle_cloud_ibs) end select end do @@ -173,6 +175,149 @@ 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) :: xmin, xmax, ymin, ymax, zmin, zmax, min_dist + real(wp) :: rx, ry, rz, dist, theta, phi, r_shell, rho, u + real(wp) :: r_inner, r_outer, xdir, ydir, zdir + logical :: overlaps, inside_bounds + 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, allocatable :: hash_head(:), chain_next(:) + + xmin = particle_cloud(cloud_idx)%x_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_x + xmax = particle_cloud(cloud_idx)%x_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_x + ymin = particle_cloud(cloud_idx)%y_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_y + ymax = particle_cloud(cloud_idx)%y_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_y + zmin = particle_cloud(cloud_idx)%z_centroid - 0.5_wp*particle_cloud(cloud_idx)%length_z + zmax = particle_cloud(cloud_idx)%z_centroid + 0.5_wp*particle_cloud(cloud_idx)%length_z + + 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 (p == 0) then + geom = 2 + dz_lo = 0 + dz_hi = 0 + else + geom = 8 + dz_lo = -1 + dz_hi = 1 + 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 (p == 0) 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)) + xdir = rho*cos(phi) + ydir = rho*sin(phi) + 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*xdir + ry = particle_cloud(cloud_idx)%y_centroid + r_shell*ydir + rz = particle_cloud(cloud_idx)%z_centroid + r_shell*zdir + end if + + inside_bounds = rx >= xmin + particle_cloud(cloud_idx)%radius .and. rx <= xmax - particle_cloud(cloud_idx)%radius & + & .and. ry >= ymin + particle_cloud(cloud_idx)%radius .and. ry <= ymax - particle_cloud(cloud_idx)%radius + if (p == 0) then + inside_bounds = inside_bounds .and. ry >= particle_cloud(cloud_idx)%y_centroid + particle_cloud(cloud_idx)%radius + else + inside_bounds = inside_bounds .and. rz >= zmin + particle_cloud(cloud_idx)%radius .and. rz <= zmax & + & - particle_cloud(cloud_idx)%radius .and. rz >= particle_cloud(cloud_idx)%z_centroid & + & + particle_cloud(cloud_idx)%radius + end if + if (.not. inside_bounds) cycle + + bx = int(floor(rx/min_dist)) + by = int(floor(ry/min_dist)) + bz = 0 + if (p /= 0) bz = int(floor(rz/min_dist)) + + 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 + n_placed = n_placed + 1 + placed(1, n_placed) = rx + placed(2, n_placed) = ry + placed(3, n_placed) = rz + + 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 diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 1099b51328..adec47a107 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -631,11 +631,11 @@ def check_ibm(self): packing_method = self.get(f"particle_cloud({i})%packing_method", None) self.prohibit( packing_method is None, - f"particle_cloud({i})%packing_method must be specified (1 = rejection sampling, 2 = lattice)", + f"particle_cloud({i})%packing_method must be specified (1 = rejection sampling, 2 = lattice, 3 = hemisphere shell)", ) self.prohibit( - packing_method is not None and packing_method not in [1, 2], - f"particle_cloud({i})%packing_method must be 1 (rejection sampling) or 2 (lattice)", + packing_method is not None and packing_method not in [1, 2, 3], + f"particle_cloud({i})%packing_method must be 1 (rejection sampling), 2 (lattice), or 3 (hemisphere shell)", ) num_ib_airfoils_max = get_fortran_constants().get("num_ib_airfoils_max", 5) diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 038ba9449c..57270d4e97 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -987,6 +987,8 @@ def _load(): _pb_attrs["radius"] = (REAL, _pb_tags) _pb_attrs["mass"] = (REAL, _pb_tags) _pb_attrs["min_spacing"] = (REAL, _pb_tags) + _pb_attrs["shell_inner_radius"] = (REAL, _pb_tags) + _pb_attrs["shell_outer_radius"] = (REAL, _pb_tags) _pb_attrs["moving_ibm"] = (INT, _pb_tags) _pb_attrs["seed"] = (INT, _pb_tags) _pb_attrs["packing_method"] = (INT, _pb_tags) @@ -1083,11 +1085,8 @@ def _load(): _r(f"simplex_params%perturb_vel_offset({d},{j})", REAL) # lag_params (Lagrangian bubbles) - # Members present in bubbles_lagrange_parameters: solver_approach, cluster_type, - # pressure_corrector, smooth_type, heatTransfer_model, massTransfer_model, - # write_bubbles, write_bubbles_stats, write_void_evol, pressure_force, - # gravity_force, nBubs_glb, epsilonb, charwidth, valmaxvoid. T0/Thost/c0/rho0/x0 - # were removed from the Fortran type by upstream #1085/#1093 — they must NOT be + # Members present in bubbles_lagrange_parameters. T0/Thost/c0/rho0/x0 were + # removed from the Fortran type by upstream #1085/#1093 — they must NOT be # registered (namelist read would crash). for a in ["heatTransfer_model", "massTransfer_model", "pressure_corrector", "write_bubbles", "write_bubbles_stats", "pressure_force", "gravity_force", "write_void_evol", "kahan_summation"]: _r(f"lag_params%{a}", LOG, {"bubbles"})