Skip to content

Add hemisphere shell particle cloud packing #1667

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
BCKim55 wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Add hemisphere shell particle cloud packing #1667

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
101 changes: 101 additions & 0 deletions examples/3D_mibm_particle_cloud_hemi_shell/case.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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,
}
)
)
4 changes: 3 additions & 1 deletion src/common/m_derived_types.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
11 changes: 8 additions & 3 deletions src/simulation/m_checker.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 2 additions & 0 deletions src/simulation/m_global_parameters.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/simulation/m_mpi_proxy.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
145 changes: 145 additions & 0 deletions src/simulation/m_particle_cloud.fpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
6 changes: 3 additions & 3 deletions toolchain/mfc/case_validator.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
9 changes: 4 additions & 5 deletions toolchain/mfc/params/definitions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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"})
Expand Down
Loading