Add hemisphere shell particle cloud packing

[BCKim55]

Description

Adds a hemisphere-shell particle cloud packing option for immersed-boundary particle clouds.

This introduces particle_cloud(i)%packing_method = 3, which randomly places spherical/circular IBM particles inside a hemisphere-shell region while enforcing:

• inner and outer shell-radius clearance,
• hemisphere plane clearance,
• bounding-box clearance,
• particle-particle non-overlap using the existing spatial hash approach.

This also adds shell_inner_radius and shell_outer_radius particle cloud parameters. The local verification example was removed from this PR to avoid introducing a new golden test in the same change.

Type of change

• New feature

Testing

• ./mfc.sh format
• ./mfc.sh validate examples/*/case.py
• ./mfc.sh validate examples/3D_mibm_particle_cloud_hemi_shell/case.py
• ./mfc.sh run examples/3D_mibm_particle_cloud_hemi_shell/case.py --clean --no-debug

Additional local checks:

• Verified generated ib_state_0.dat particle positions.
• Confirmed no shell-boundary violations.
• Confirmed no particle overlap violations.
• Confirmed VF 0.2 succeeds with min_spacing=0.02.
• Confirmed VF 0.3 succeeds when min_spacing=0.0.

Checklist

• I added or updated tests for new behavior
• I updated documentation if user-facing behavior changed

GPU changes (expand if you modified src/simulation/)

• GPU results match CPU results
• Tested on NVIDIA GPU or AMD GPU

[wilfonba]

I don't think this is implemented quite right. The example you added highlights this fact. It defines a hemispherical particle cloud using a cubic domain. I think instead you should add a particle_cloud(i)%geometry variable, and then treat the particle cloud more like the fluid patches are treated. Also, a hemispherical shell is not a packing method in the same way that rejection sampling and lattice placement are packing methods, so using the existing acking_method variable for this isn't the right thing to do. You can place particles in a hemispherical shell using either rejection sampling, lattice placement, or some other packing method that's not implemented yet.

Edit: Your s_particle_cloud_random_hemi_shell routine appears to do rejection sampling, which adds to the confusion of what packing_method means.

[BCKim55]

Thanks Ben, that makes sense. I agree that the hemisphere shell should be treated as a particle-cloud geometry rather than a packing method.

I’ll refactor this so that particle_cloud(i)%packing_method remains responsible for the placement algorithm, e.g. rejection sampling or lattice, and add a separate particle_cloud(i)%geometry option for the cloud shape. Then the current hemisphere-shell implementation will become the rejection-sampling path for geometry = hemisphere shell, while unsupported combinations such as hemisphere-shell lattice packing can error out explicitly for now.

I’ll also remove or adjust the example so it does not imply that a cubic domain defines the hemispherical cloud.

[danieljvickers]

num_dims < 3 peferred in general. This likely does not affect things, but num_dims is a case-optimization parameter, and therefore the compiler can optimize this away if you make it a check on num_dims.

[BCKim55]

Updated this to use num_dims < 3 instead of checking p == 0.

[danieljvickers]

This seems like a lot fo compute for your QOI. You do not use the [xyz]dir parameter at all after computing it. Again, either use it or do not computed.

[danieljvickers]

That said, you using a probability distribution function to weight your random seed is correct here. Just you are computing redundant quantities here that seem pointless.

[BCKim55]

1. Thanks for catching that! I completely missed that [xyz]dir wasn't being used later on. I've removed the redundant calculation to save compute.

2. Thanks for confirming the PDF logic. I've cleaned up all the redundant quantities you mentioned to optimize the compute.

[BCKim55]

Thanks for the feedback. I updated the hemisphere-shell path to remove the redundant box-bounds logic and rely on the shell radii for the placement region. I also changed the dimensionality checks to use num_dims < 3 and removed the extra direction variables in the 3D sampling path.

[BCKim55]

Updated this to use num_dims < 3 instead of checking p == 0.

[BCKim55]

1. Thanks for catching that! I completely missed that [xyz]dir wasn't being used later on. I've removed the redundant calculation to save compute.

2. Thanks for confirming the PDF logic. I've cleaned up all the redundant quantities you mentioned to optimize the compute.

[sbryngelson]

Request changes

1. Missing regression coverage. This PR adds the hemisphere-shell execution path but removes the local example and adds no test fixture. The existing particle-cloud case is box-only, so CI never executes geometry = 2. Please add a small deterministic 2D or 3D case that checks generated particle positions for shell/plane clearance and non-overlap.

2. Documentation was not updated. docs/documentation/case.md does not document geometry, shell_inner_radius, or shell_outer_radius, and its length_[x,y,z] description is inaccurate for geometry = 2, where those extents are ignored. Please update the Particle Clouds section accordingly.

[BCKim55]

Hi Spencer, I addressed the requested validation, documentation, and regression coverage changes. CI is now passing, and the PR is ready for another look when you have time.

[codecov]

Codecov Report

❌ Patch coverage is 68.51852% with 34 lines in your changes missing coverage. Please review.
✅ Project coverage is 61.67%. Comparing base (b7d7883) to head (60e4afa).

Files with missing lines	Patch %	Lines
src/simulation/m_particle_cloud.fpp	70.00%	18 Missing and 12 partials ⚠️
src/simulation/m_checker.fpp	0.00%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #1667      +/-   ##
==========================================
+ Coverage   61.63%   61.67%   +0.03%     
==========================================
  Files          84       84              
  Lines       21511    21600      +89     
  Branches     3174     3185      +11     
==========================================
+ Hits        13258    13321      +63     
- Misses       6074     6090      +16     
- Partials     2179     2189      +10     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[sbryngelson]

Thanks for the revisions. Documentation and validation are genuinely addressed:

• case.md now documents geometry, shell_inner_radius, shell_outer_radius, and qualifies length_x[y,z] as box-only.
• The Python validator and m_checker.fpp now enforce the same condition, and shell_outer_radius > shell_inner_radius + 2*radius is exactly the sampler's feasibility requirement (r_outer > r_inner after the +/- radius insets), so the s_mpi_abort path is no longer reachable from valid input.
• geometry and both new reals are broadcast in m_mpi_proxy.fpp.
• The sampling math is right: the 2D sqrt CDF and the 3D cube-root CDF with uniform cos(polar) both give uniform density in the shell.

The test coverage point is not yet addressed, though, and that is the blocking item. Details inline.

[sbryngelson]

This test does not exercise the code the PR adds. _pack_hemi_shell_3d is a Python reimplementation of s_particle_cloud_random_hemi_shell, and the assertions below check that the Python implementation satisfies the shell and overlap constraints. Nothing here reads from m_particle_cloud.fpp, so no regression in the Fortran sampler can make this test fail.

What I asked for was verification of the generated particle positions. Please parse ib_state_0.dat from an actual run and assert the shell clearance, plane clearance, and pairwise non-overlap on those coordinates. A mirrored implementation cannot serve that purpose no matter how faithful it is.

[BCKim55]

Addressed. Instead of validating a Python mirror implementation, I now run the actual MFC particle-cloud cases and validate the Fortran-generated ib_state_0.dat output. The test checks that generated particles satisfy the box/hemi-shell bounds and the minimum non-overlap distance.

[sbryngelson]

Separately from the point above: this model is not actually equivalent to the Fortran. _is_overlapping keys cells on exact bin tuples, whereas f_cloud_particle_overlaps folds bins through f_bin_hash into hash_size buckets. Collisions cause the Fortran to reject candidates this model accepts, so the two will diverge on particle placement for the same seed.

So even taken as documentation of the algorithm, it is misleading. Another reason to test the real output instead.

[BCKim55]

Agreed. The Python mirror was not equivalent to the Fortran spatial hash implementation because it keyed cells directly rather than folding them through f_bin_hash, so it could accept candidates that the Fortran code rejects after hash collisions.

I removed the mirror implementation entirely. The tests now run the actual MFC/Fortran particle-cloud generation and validate the generated restart_data/ib_state_0.dat output directly for geometry bounds, hemi-shell clearance, and particle-particle non-overlap.

[sbryngelson]

Two coverage gaps here.

1. 3D is not covered. This block sits inside if len(dimInfo[0]) == 2 and not viscous:, so CI only ever runs the 2D branch — which is a half-annulus, not a hemisphere. The 3D path (phi/zdir/rho with the cube-root radial CDF) is the one this feature is named for and it never executes. Please add a 3D golden.

2. The box path this PR rewrote has no coverage at all. grep particle_cloud toolchain/mfc/test/cases.py returns only this new case — it is the only particle-cloud golden in the suite. Meanwhile s_particle_cloud_random_box was rewritten here (p == 0 -> num_dims < 3, the inlined overlap loop replaced by f_cloud_particle_overlaps, bin coordinates now computed after acceptance rather than before). That refactor is unverified. Please add a box golden as well, or demonstrate bit-identical output before and after the change on an existing example.

Note also that the goldens hold cons.* fields rather than ib_state. Placement perturbs the flow, so this is a real regression lock, but it does not directly assert the geometric constraints.

[BCKim55]

Thanks, this was a good catch.

1. I added 3D particle-cloud regression coverage for the actual hemisphere-shell path. The test matrix now includes 3D hemisphere-shell particle-cloud cases, so the phi/zdir/rho sampling and cube-root radial CDF branch is exercised in CI.

2. I also added box particle-cloud regression coverage for the refactored random-box path. The test matrix now includes box particle-cloud cases, so the shared overlap-check refactor and num_dims-based 2D/3D logic are covered as well.

I also enabled ib_state_wrt for these particle-cloud cases and added a verifier that reads restart_data/ib_state_0.dat directly. The regular goldens still lock the resulting flow-field regression, while the new verifier directly checks the particle geometry constraints: box bounds, shell radial clearance, shell flat-plane clearance, and minimum particle-particle spacing.

[sbryngelson]

geometry is overloaded in a way that will bite users. Here particle_cloud(i)%geometry == 2 means hemisphere shell, while two lines down geom = 2 means circle in the patch_ib vocabulary (and geom = 8 means sphere). The same subroutine uses two different meanings of geometry = 2 within six lines.

Please rename the cloud field to something like region_shape or cloud_geometry, or reuse the existing patch geometry codes. As written, a user who reads particle_cloud(1)%geometry = 2 will reasonably expect a circle.

[BCKim55]

Addressed. I renamed particle_cloud%geometry to particle_cloud%cloud_geometry throughout the Fortran code, parameter registry, validator, documentation, and tests to avoid confusion with patch_ib%geometry.

[BCKim55]

Follow-up: after the later rename, this is now the cloud_geometry = 2 path rather than particle_cloud%geometry = 2.

[sbryngelson]

Two things about the region definition.

Orientation is hard-coded and unconfigurable. The flat face is at y_centroid in 2D and z_centroid in 3D, always opening toward +y / +z. That is a reasonable initial restriction, but it needs to be stated in case.md — right now nothing tells a user which way the hemisphere faces.

Nothing bounds the cloud to the domain. Dropping the [xyz]min/max logic was the right call per @danieljvickers, but nothing replaced it: x_centroid +/- shell_outer_radius can extend past x_domain, and s_add_cloud_particle performs no bounds check, so particles land off-grid silently. The validator already has x[y,z]_centroid, shell_outer_radius, and the domain extents, so this is a cheap check to add there.

[BCKim55]

Addressed. I expanded the documentation to define the hemi-shell orientation and domain assumptions explicitly. In 2D, this is treated as a half-annulus with the flat face at y_centroid and opening in +y. In 3D, this is a hemisphere shell with the flat face at z_centroid and opening in +z.

I also added validator checks to require the full shell extent to remain inside the computational domain.

[sbryngelson]

Minor: the extraction turned integer exponents into real ones. The original was (rx - placed(1, j))**2; this is now **2._wp, and min_dist_sq = min_dist**2._wp above is the same. gfortran folds x**2.0 to x*x (I checked at both -O0 and -O2), but that is a compiler courtesy rather than a guarantee, and this is the innermost loop of the rejection sampler. Please use **2.

Also minor, in the same area: s_get_cloud_bin runs twice per accepted particle — once inside f_cloud_particle_overlaps and again at the insertion site on line 227 — and it tests num_dims == 3 while everything else in the module tests num_dims < 3.

[BCKim55]

Addressed. I changed the distance calculations back to integer exponents, e.g. min_dist**2 and coordinate differences **2.

I also refactored the overlap check so it returns the candidate hash-bin indices along with the overlap flag. That way, accepted particles reuse the same bx/by/bz values at insertion instead of calling s_get_cloud_bin a second time. I also updated the dimensionality check to use num_dims < 3 for the 2D path.

[sbryngelson]

This is a clear improvement, but two facts a user needs are still missing:

• Which way the hemisphere opens. The flat face lies at y_centroid in 2D and z_centroid in 3D, and the filled half is +y / +z. This is not configurable and is not stated anywhere.
• That the 2D region is a half-annulus, not a hemisphere. Calling both cases "hemisphere shell" is going to confuse people setting up 2D cases.

[BCKim55]

This overlaps with the geometry/orientation documentation issue discussed above, so I addressed it as part of that same update.

The documentation now explicitly distinguishes the 2D case as a half-annulus and the 3D case as a hemisphere shell. It also states the fixed orientation: the 2D half-annulus opens in +y from y_centroid, and the 3D hemisphere shell opens in +z from z_centroid.

[BCKim55]

Hi Spencer,

I’ve pushed the latest fixes and the PR is ready for re-review.

The blocking regression coverage item is now addressed for particle clouds. The new checks validate generated ib_state_0.dat particle count, shell/plane clearance, box bounds, and non-overlap. I also fixed the single-precision test reader issue that caused the Ubuntu GNU no-debug CI failure.

All GitHub checks are now passing: 60 success, 5 skipped, 0 failures.

Thanks again!

[sbryngelson]

Thanks for the revisions. The design is now right, and everything from my previous two rounds landed:

• The geometry to cloud_geometry rename resolves the collision with geom = 2 (circle) and geom = 8 (sphere) in the same routine.
• Integer exponents are back in the innermost loop, and s_get_cloud_bin now runs once per candidate.
• The Python mirror implementation is gone in favour of parsing real ib_state_0.dat output.
• 2D, 3D, and box goldens all exist. All six golden dirs are new, so the squared-distance refactor did not perturb the box path.
• Sampling math verified independently: the 2D sqrt radial CDF and the 3D cube-root CDF with uniform cos(polar) both give uniform density, and the plane-clearance cut preserves uniformity on the remaining region.

./mfc.sh precheck passes all 7 checks, the 4 new pytest cases pass, and ./mfc.sh test -l lists exactly the 6 new cases matching the 6 new golden dirs.

What is left is one test-fixture problem that undermines the 3D coverage I asked for, a validation-surface gap, and two latent bugs. Details inline.

Blocking

1. The 3D particles are sub-cell (radius = 0.6*dx), so the 3D golden may not be observing the feature it was added to cover. The 2D goldens do carry real signal; only the 3D case needs fixing.
2. The box branch of the new ib_state assertion is only valid for packing_method = 1. Separately, this exposed what looks like a genuine pre-existing bug in s_particle_cloud_lattice: its outer loops have no ymax/zmax bound, so lattice clouds can place particles outside the declared region and potentially outside the domain. Not yours to fix here, but it should get an issue.
3. The m_checker shell-thickness guard is vacuous when radius is unset, since dflt_real = -1.e6_wp swamps the right-hand side. More broadly, nothing requires radius, mass, or num_particles to be set for any cloud.
4. Domain containment is enforced only in case_validator.py. Please mirror it into m_checker, same as you did for the radii checks.
5. cloud_geometry = 2 is unguarded in 1D.

Design, and the one I would push hardest on

s_particle_cloud_random_hemi_shell is still a near-verbatim copy of s_particle_cloud_random_box. Extracting the overlap check was the right first step, but roughly 60 of its 110 lines are still duplicated. Collapsing both into a single rejection packer with a geometry select case inside the loop removes that and flattens the dispatch in s_generate_particle_clouds. The cloud_geometry enum invites a third geometry, and a third geometry triples the copy.

Housekeeping

The PR description is stale. It still describes particle_cloud(i)%packing_method = 3, which is the design @wilfonba asked you to change, and it lists examples/3D_mibm_particle_cloud_hemi_shell/case.py under Testing though no such example is in the PR. Release notes get assembled from PR bodies, so please fix it before merge.

[sbryngelson]

The 3D particles are sub-cell, so this golden may not be observing the feature.

This case sets m=n=p=29 on [0,1]^3, so dx = 1/30 = 0.0333 and radius = 0.02 = 0.6*dx. A sphere smaller than one cell contains at most one cell centre, sometimes zero depending on where the seed lands it relative to a cell centre. Two problems follow:

• If a particle marks no cells, the golden is blind to the IB it was added to cover. The 3D path executes but nothing observes it.
• Where a particle does mark a cell, the marked set is a discontinuous function of position, so it can flip across the compiler and precision matrix this golden runs on.

For contrast, the two 2D goldens (tests/E085CC5A vs tests/5A22B45F) differ by up to 0.37 at t_step 50, so the 2D cases do carry real signal at radius = 1*dx. The 3D case is the one that needs attention.

Please raise radius (or refine the grid) so particles span 3-4 cells, adjusting shell_inner_radius/shell_outer_radius to keep the packing feasible. A quick way to confirm the fix: check that the 3D golden actually differs from the same case with num_particle_clouds = 0.

[sbryngelson]

This box assertion is only valid for packing_method = 1, and the lattice packer looks like it has a real bug.

This branch asserts every particle centre lies within centroid +/- length/2. That invariant holds for rejection sampling but not for lattice packing. s_particle_cloud_lattice has no upper bound on its outer loops:

• 2D: do while (n_placed < n_target) increments row with py = ymin + real(row, wp)*row_dy and there is no py <= ymax guard.
• 3D: kz increments with no zmax guard.

spacing is derived so that the area/volume per particle yields exactly n_target sites in the region, but edge effects mean the sites that actually fit are generally fewer than n_target. So the loop runs extra rows or layers past ymax/zmax. That is the common case, not a corner case.

Two separate things to do:

1. In this PR: gate the box-bounds assertion on packing_method == 1. As written, adding ib_state_wrt: "T" to any lattice box-cloud case fails this assertion even though the harness is fine.
2. Separately: s_particle_cloud_lattice silently places particles outside the requested region, and potentially outside the domain. That is pre-existing and not this PR's to fix, but it is worth its own issue. This assertion is what exposed it.

[sbryngelson]

Correcting myself on the mechanism above, after actually simulating both loops rather than reasoning about them.

2D does not escape. I swept the triangular-lattice branch over several aspect ratios at n from 4 to 200 and py stayed inside ymax every time. The unbounded row loop looks wrong, but the in-region site count exceeds n_target in practice (the floor(...)+1 edge effects work in favour), so it terminates before escaping. My "that is the common case" claim was wrong for 2D.

3D does escape, but not via the unbounded kz loop. The cause is ncx = max(1, ceiling((xmax - xmin)/cell)) rounding the cell grid up past xmax, combined with the FCC basis +0.5*cell offsets landing sites in that overhang. On a unit region, num_particles = 10 puts particles at 1.105 against a bound of 1.0, and num_particles = 50 at 1.077. Anisotropic regions escape in y as well.

Filed as #1730 with a standalone repro.

None of this changes what I am asking for in this PR: please gate the box-bounds assertion on packing_method == 1. That is still correct, since lattice packing can violate the invariant this branch asserts. The rest is not yours to fix here.

[sbryngelson]

This guard is vacuous when radius is unset.

Nothing in s_check_inputs_particle_clouds requires radius to be specified, so an unset radius is dflt_real = -1.e6_wp. The right-hand side becomes roughly shell_inner_radius - 2e6, and any positive shell_outer_radius passes. The run then dies inside s_particle_cloud_random_hemi_shell on the generic "Invalid hemisphere-shell radii for particle cloud packing" instead of naming the parameter the user actually forgot.

The Python validator catches this (radius is None), so it only bites the direct-binary path, which is the same asymmetry as the domain-containment comment below.

The broader version of this: neither m_checker nor case_validator requires radius, mass, or num_particles to be set for any cloud, box or shell. An unset radius also makes min_dist negative in the box path, which corrupts both int(floor(px/min_dist)) and the dist_sq < min_dist_sq test. Please add unconditional positivity PROHIBITs for all three, independent of cloud_geometry.

[sbryngelson]

Domain containment is enforced only in Python.

case_validator.py gained the x/y/z extent checks, but nothing equivalent landed here. MFC supports running the binaries against a hand-written input file, and that path skips the toolchain entirely, so a shell that extends past the domain is accepted with no diagnostic.

This is the same gap I raised for the radii feasibility checks, which you correctly mirrored into m_checker. Please mirror the containment checks too.

[sbryngelson]

cloud_geometry = 2 is unguarded in 1D.

With n == 0, cloud_geometry = 2 falls into the num_dims < 3 branch of the sampler and generates a half-annulus in a 1D run. Nothing rejects it here, and the Python validator's y-extent check is explicitly skipped when n == 0, so the one check that would have caught it is skipped precisely where it is needed.

Please add @:PROHIBIT(particle_cloud(i)%cloud_geometry == 2 .and. n == 0, ...).

[sbryngelson]

The ib_state record layout is hardcoded here with no comment tying it to the writer.

The 20 and the record[16:19] position slice below come from m_data_output.fpp: NFIELDS_PER_IB = 20, and ib_buf(17:19) holds x/y/z_centroid. Nothing in this file says so. If ib_patch_parameters gains a field or ib_buf is reordered, this reads the wrong columns and either passes vacuously or fails with a message that points nowhere useful.

Please name the constants (NFIELDS_PER_IB, a CENTROID_SLICE) and add a comment citing s_write_serial_ib_state as the source of truth.

Related, and worth a guard: this also assumes the non-file_per_process layout. Under file_per_process = T the writer emits restart_data/lustre_0/ib_state_0_0000000.dat, with a leading num_local_ibs integer and a gbl_patch_id integer before each 20-real record. A particle-cloud case with that flag set would abort on "Expected IB state file does not exist" rather than skipping the check. Add file_per_process to the early-return guard in _assert_particle_cloud_ib_state.

[sbryngelson]

The flat face sits exactly on the domain wall.

y_centroid = 0.0 with y_domain%beg = 0.0 (set at cases.py:222). The plane-clearance rule then permits a particle tangent to the boundary: centre at y = radius = 0.02, surface at y = 0. Tangency is a degenerate configuration for IB ghost-point reconstruction and makes this golden more sensitive than it needs to be.

The 3D case sensibly uses z_centroid = 0.1. Please give the 2D case the same standoff.

Note also that the validator's y_centroid < y_beg check passes at exact equality, which is how this got through. Requiring at least one particle radius of clearance there would be more useful than requiring none.

[sbryngelson]

At 4 particles, hash_size = max(16, 4*4) = 16 and every bucket chain has length at most 1. The shared s_check_cloud_particle_overlap these cases were added to cover is therefore only exercised in its trivial regime: the do while (j > 0) chain walk never iterates, and no hash collision ever occurs.

One case at roughly 50 particles would exercise the chain walk and the collision handling, which is where the refactor could actually regress. It would also give the rejection loop something to reject, which nothing currently does.

[sbryngelson]

Minor: this comparison and the shell_outer_radius <= shell_inner_radius + 2 * radius one below are evaluated without the self._is_numeric(...) guard the rest of this file applies consistently to real-valued params. The domain-extent block ten lines down does guard.

None is short-circuited correctly, but any non-numeric value raises an uncaught TypeError out of validate() instead of producing a clean CaseConstraintError.

[sbryngelson]

Unrelated to this PR. This reflows a lag_params comment block that has nothing to do with particle clouds. Please drop it so the diff stays scoped.