CUT-CELL / EMBEDDED BOUNDARY (EB) spatial operator: R = -div_eb F + S on a disc, in CONSERVATIVE finite volumes with face apertures alpha_f and volume fraction kappa derived from detail::cut_fraction (project T5-PR1). More...

#include <pops/core/state/state.hpp>
#include <pops/core/foundation/types.hpp>
#include <pops/mesh/layout/box_array.hpp>
#include <pops/mesh/storage/fab2d.hpp>
#include <pops/mesh/execution/for_each.hpp>
#include <pops/mesh/geometry/geometry.hpp>
#include <pops/mesh/storage/multifab.hpp>
#include <pops/numerics/elliptic/eb/cut_fraction.hpp>
#include <pops/numerics/fv/numerical_flux.hpp>
#include <pops/numerics/fv/reconstruction.hpp>
#include <pops/numerics/spatial/embedded_boundary/domain.hpp>
#include <pops/numerics/spatial_operator.hpp>
#include <utility>
#include <vector>

Include dependency graph for operator.hpp:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	pops::detail::DiscLevelSet
	Device-safe adapter wrapping a DiscDomain as the Real(Real, Real) callable expected by cut_fraction and the EB operator. More...

struct	pops::detail::EbFaceFluxXKernel< Limiter, NumericalFlux, Model, LevelSet >
	FACE FLUX kernel for x (dir 0) of the EB transport: numerical flux at the face between (i-1, j) and (i, j), WEIGHTED by the aperture alpha_x of that face. More...

struct	pops::detail::EbFaceFluxYKernel< Limiter, NumericalFlux, Model, LevelSet >
	FACE FLUX kernel for y (dir 1) of the EB transport: analogue of EbFaceFluxXKernel in j. More...

struct	pops::detail::EbAssembleRhsKernel< Model, LevelSet >
	Kernel assembling the EB residual at cell (i, j): INACTIVE cell -> residual 0 (not advanced, like T2); ACTIVE cell -> R = S - (1/kappa_eff) [ (fx_{i+1} - fx_i)/dx + (fy_{j+1} - fy_i)/dy ] - wall_term. More...

Namespaces
namespace	pops

namespace	pops::detail

Functions
POPS_HD DiscLevelSet	pops::detail::disc_level_set (const DiscDomain &d)
	Builds the disc level set callable from a DiscDomain (sugar: disc_level_set(d)).

template<class LevelSet >
POPS_HD bool	pops::detail::eb_cell_active (const LevelSet &ls, Real xc, Real yc)
	Activity indicator (center in the disc, ls < 0) from a callable level set. POPS_HD.

POPS_HD Real	pops::detail::eb_face_aperture (Real lc, Real ln, Real h)

template<class Limiter = NoSlope, class NumericalFlux = RusanovFlux, class Model , class LevelSet >
void	pops::assemble_rhs_eb (const Model &model, const MultiFab &U, const MultiFab &aux, const LevelSet &ls, const Geometry &geom, MultiFab &R, bool recon_prim=false, Real kappa_min=detail::kEbKappaMin, Real pos_floor=Real(0))
	assemble_rhs_eb<Limiter, NumericalFlux>: residual R = -div_eb F + S on a DISC in cut-cell / EB, with face apertures alpha_f in [0, 1] and volume fraction kappa derived from detail::cut_fraction (T5-PR1).

Variables
constexpr Real	pops::detail::kEbFaceOpenEps = Real(1e-6)
	Default aperture below which a face is treated as CLOSED (immersed wall).

constexpr Real	pops::detail::kEbKappaMin = Real(1e-2)

Detailed Description

CUT-CELL / EMBEDDED BOUNDARY (EB) spatial operator: R = -div_eb F + S on a disc, in CONSERVATIVE finite volumes with face apertures alpha_f and volume fraction kappa derived from detail::cut_fraction (project T5-PR1).

Aperture of ONE face between the active cell (xc, yc) and its neighbor at (xn, yn), step h.

Volume fraction floor (small-cell clamp).

CONTEXT (the "Cartesian-ring-edge lock"; cf. docs/HOFFART_FIDELITY.md). The T2 path (spatial_operator.hpp: assemble_rhs_masked) approximates the disc with a STAIRCASE MASK: an active/inactive face is a 0/1 GATE (normal flux set to zero), the boundary is crenellated. This operator GENERALIZES the 0/1 gate to an APERTURE alpha_f in [0, 1] (the linear fraction of the face inside the disc, EXACTLY the cut_distance/h of the elliptic wall -> bit-for-bit consistency with Poisson) AND divides the residual by the cell's VOLUME FRACTION kappa (volume of the cut cell). This is a 2nd-order EB scheme for smooth transport inside the disc, which recovers the physical interior growth rate without the structural l-dependent over-rate of the staircase (Cartesian) baseline.

NOTE ON alpha_f AND THE ACTIVITY MODEL. Cell activity follows the center-in-the-disc criterion (ls(center) < 0), like GeometricMG and the T2 mask. With this criterion, the geometric aperture cut_distance(lc, ln, h)/h equals EXACTLY 1 on a face between two ACTIVE cells (lc < 0 and ln < 0 -> "interior neighbor" branch of cut_distance -> distance = h) and FRACTIONAL in (0, 1) on the leg from an active cell toward an INACTIVE neighbor (ln >= 0: linear crossing). The immersed boundary is therefore carried by TWO continuous quantities derived from cut_fraction: (i) the fractional leg (which closes the face toward the inactive side, cf. no-penetration wall below) and ABOVE ALL (ii) the VOLUME FRACTION kappa in (0, 1] of the cut cell, which corrects the divergence by the true immersed volume. It is kappa (and not the aperture of internal faces) that de-crenellates the boundary and restores 2nd-order interior transport, which the T2 mask (implicit volume = 1) does not.

CONSERVATIVE EB FORM (cell (i, j), volume kappa dx dy): kappa dx dy d_t U = - [ alpha_xp Fx_{i+1} - alpha_xm Fx_i ] dy

[ alpha_yp Fy_{j+1} - alpha_ym Fy_i ] dx
alpha_wall |wall| F_wall (immersed WALL term)
kappa dx dy S that is, after dividing by kappa dx dy (with kappa CLAMPED, cf. small-cell stability): R = S - (1/kappa) [ (alpha_xp Fx_{i+1} - alpha_xm Fx_i) / dx

(alpha_yp Fy_{j+1} - alpha_ym Fy_i) / dy ]

(1/kappa) (alpha_wall |wall| / (dx dy)) F_wall The immersed WALL flux is a NO-PENETRATION flux (zero-normal-flux): F_wall = 0. It is the FV counterpart of the conducting wall (the elliptic side applies Dirichlet; the transport applies a solid wall at the SAME geometric boundary). The wall term is therefore IDENTICALLY ZERO; it is written explicitly (and kept at zero) so the contract stays readable and so that a future nonzero wall flux (injection, slip) has its single attachment point.

MASS CONSERVATION. On two neighboring active cells, the aperture of the SHARED FACE is the SAME on both sides (geometric alpha_f, a function of the level set only: alpha_xp(i) == alpha_xm(i+1), cf. cut_fraction symmetric per face) -> the flux of that face telescopes EXACTLY in the sum Sum_cells kappa dx dy R. A face whose neighbor cell is INACTIVE is CLOSED (alpha_f forced to 0) AND the wall flux is zero -> no mass crosses the immersed boundary. The total mass over the active cells Sum n_ij kappa_ij dx dy is therefore conserved TO MACHINE PRECISION. This is the EB analogue (continuous apertures) of the T2 conservation (0/1 gates).

SMALL-CELL STABILITY (small-cell problem). The 1/kappa factor amplifies the residual when kappa becomes small on the r0/r1 shear layer; at a FIXED time step dt calibrated on full cells, an unbounded amplification blows up (Inf/NaN) the explicit step on a strongly cut cell. Two STACKED complementary guards:

The PR1 primitive floor: cut_distance clamps each half-face at theta >= 1e-3 (anti-division guard INHERITED from the elliptic wall). kappa = product of half-face averages can therefore never equal EXACTLY 0 -> no strict division by zero. BUT the induced lower bound is ~ (1e-3)^2 / 4 ~ 2.5e-7: 1/kappa can reach ~4e6, which is enough to overflow the fixed step.
RETAINED SCHEME = VOLUME CLAMP (this header): kappa_eff = max(kappa, kappa_min), kappa_min by default 1e-2 (1% of the full volume). A SCHEME-LEVEL guard, INDEPENDENT of the elliptic floor: it bounds the 1/kappa amplification to 1/kappa_min = 100, a value CALIBRATED so the fixed explicit step stays stable whatever the degree of cut. This is the simplest and most robust implicit "volume merging" for a FIXED step, at the cost of a slight LOCAL NON-conservation on the most-cut cells (effective volume > real volume). The GLOBAL mass stays conserved TO MACHINE PRECISION because the clamp acts only on the DENOMINATOR (volume), NOT on the face fluxes (numerator): the telescoping sum of fluxes is unchanged (the discrete mass consistent with the scheme uses the SAME kappa_eff, cf. conservation test). Documented alternative (outside PR2): flux redistribution (AMReX-EB's flux redistribution) spreads the excess divergence of small cells onto the full neighbors -> exact LOCAL conservation but a non-local stencil; the clamp suffices for the target (calibrated fixed step, smooth MMS, global conservation).

NAMED FUNCTORS (and not extended lambdas), like spatial_operator.hpp / _polar.hpp (#64/#97): robust device emission when the Model-template kernel is instantiated cross-TU. The RECONSTRUCTION and the numerical FLUX are REUSED verbatim from the Cartesian operator (reconstruct<>, RusanovFlux). The level set is passed BY VALUE (POPS_HD callable, e.g. captured detail::DiscDomain): device-safe.

INVARIANT: this header is PURELY ADDITIVE and OPT-IN. The Cartesian operator (assemble_rhs) and the T2 mask path (assemble_rhs_masked) stay STRICTLY UNTOUCHED; a run without an EB disc is bit-identical. assemble_rhs_eb is called only on explicit opt-in (geometry = cutcell).

kappa_eff = max(kappa, kEbKappaMin) bounds the 1/kappa amplification to 1/kEbKappaMin = 100 -> finite residual on an arbitrarily cut cell, stable fixed explicit step. See the SMALL-CELL STABILITY note of

FV EB convention:

INACTIVE neighbor (ls(xn,yn) >= 0): the face touches the immersed wall -> CLOSED (alpha = 0, no-penetration). This is the generalization of the T2 0/1 gate: the inactive side closes the face.
ACTIVE neighbor (ln < 0): the aperture reuses VERBATIM the shared primitive cut_distance (hence bit-consistent with the elliptic wall), alpha = cut_distance(lc, ln, h) / h. But for an active neighbor cut_distance takes the "interior neighbor" branch and returns h -> alpha = 1 EXACTLY, far and near the edge alike: the shared face of two active cells is always FULL. The face apertures are therefore BINARY {0, 1}; it is the VOLUME FRACTION kappa in (0, 1], and not the aperture of internal faces, that carries the cut geometry (cf. NOTE alpha_f of SYMMETRY (key to conservation): cut_distance(lc, ln, h) depends only on (lc, ln); the shared face seen from cell i (center lc, neighbor ln) and seen from cell i+1 (center ln, neighbor lc) gives the SAME aperture as soon as both are active (ln < 0: the "interior neighbor" branch returns h on both sides -> alpha = 1). The internal active/active boundary is therefore treated SYMMETRICALLY.

include/pops/numerics/spatial/embedded_boundary/operator.hpp File Reference¶

Classes

Namespaces

Functions

Variables

Detailed Description