EB
These tutorials cover embedded boundaries (EB): building cut-cell
geometries from implicit functions and STL files, solving elliptic problems
with boundary conditions on the EB, and conservative cut-cell advection.
All of them need AMReX built with EB support (Poisson and
Advection also need the linear solvers) and run in 2D and 3D only. The
CMake superbuild enables both by default. Build only
this topic with -DAMReXTutorials_TOPICS=EB; see Getting Started.
Geometry
- Level:
★
- Source:
- Key APIs:
EB2::SphereIF,EB2::BoxIF,EB2::CylinderIF,EB2::PlaneIF,EB2::ParserIF,EB2::makeUnion,EB2::makeIntersection,EB2::makeDifference,EB2::rotate,EB2::translate,EB2::makeShop,EB2::Build,makeEBFabFactory,EBFArrayBoxFactory::getEBDataArrays,EBData_t,EB_WriteSingleLevelPlotfile,WriteEBSurface
Build EB geometries on \([-1,1]^D\) and query the cut-cell data. The
program builds each geometry listed in geometries:
sphere: a sphere (a circle in 2D) with the fluid inside.csg: constructive solid geometry, a box with a cylindrical hole, rotated and translated, plus a hemisphere made by intersecting a ball with a plane.parser: an implicit function given as a runtime expression,parser.function.stl: a triangulated surface fromoctahedron.stl, read through theeb2.*runtime parameters.
An implicit function is negative in the fluid and positive in the body, and the CSG operations act on the bodies:
// Body 1: a box with a cylindrical hole, rotated and translated
const Real a = 0.35_rt; // half width of the box
const Real rc = 0.15_rt; // radius of the hole
EB2::BoxIF box({AMREX_D_DECL(-a,-a,-a)}, {AMREX_D_DECL(a,a,a)}, false);
#if (AMREX_SPACEDIM == 3)
EB2::CylinderIF hole(rc, 2, {0.0_rt,0.0_rt,0.0_rt}, false); // along z
#else
EB2::SphereIF hole(rc, {0.0_rt,0.0_rt}, false); // a disk in 2D
#endif
auto block = EB2::translate(EB2::rotate(EB2::makeDifference(box, hole),
pi/6, 2),
{AMREX_D_DECL(-0.35_rt,-0.35_rt,0.0_rt)});
// Body 2: a hemisphere, the intersection of a ball and a half space
const Real rs = 0.25_rt;
EB2::SphereIF ball(rs, {AMREX_D_DECL(0.5_rt,0.5_rt,0.0_rt)}, false);
EB2::PlaneIF half_space({AMREX_D_DECL(0.5_rt,0.5_rt,0.0_rt)},
{AMREX_D_DECL(0.0_rt,1.0_rt,0.0_rt)});
auto dome = EB2::makeIntersection(ball, half_space);
auto bodies = EB2::makeUnion(block, dome);
EB2::Build(EB2::makeShop(bodies), geom, required_coarsening_level,
max_coarsening_level);
A Parser expression in \(x, y, z\) becomes an implicit function
with EB2::ParserIF. Without a GeometryShop, EB2::Build reads the
geometry from the eb2.* parameters (e.g., eb2.geom_type = stl):
Parser parser(function);
parser.registerVariables({"x","y","z"});
EB2::ParserIF pif(parser.compile<3>());
// The GeometryShop keeps the Parser alive while the EB is built.
EB2::Build(EB2::makeShop(pif, parser), geom, required_coarsening_level,
max_coarsening_level);
// Ellipse/ellipsoid with the fluid inside
const Real exact = (AMREX_SPACEDIM == 2) ? pi*0.8_rt*0.5_rt
: Real(4./3.)*pi*0.8_rt*0.5_rt*0.4_rt;
return user_function ? Real(-1.) : exact;
}
else if (name == "stl")
{
// Without a GeometryShop, EB2::Build reads the eb2.* parameters,
// e.g., eb2.geom_type = stl and eb2.stl_file.
EB2::Build(geom, required_coarsening_level, max_coarsening_level);
For each geometry, the program counts regular, cut and covered cells, and
sums the fluid volume, the EB surface area, and the fluid area of the plane
\(x = 0\). The EB data (cell flags, volume and area fractions,
centroids, …) are accessed in a kernel through the factory’s
EBDataArrays:
// The factory gives access to the EB data on the grids: cell
// flags, volume fractions, area fractions, centroids, ...
auto factory = makeEBFabFactory(geom, ba, dm, {1,1,1}, EBSupport::full);
MultiFab const& vfrac = factory->getVolFrac();
auto const& ebdata = factory->getEBDataArrays();
// Count cell types, sum fluid volume and EB area, and measure
// the fluid cross section at x = 0 with the x-area fractions.
auto r = ParReduce(TypeList<ReduceOpSum,ReduceOpSum,ReduceOpSum,
ReduceOpSum,ReduceOpSum,ReduceOpSum>{},
TypeList<Long,Long,Long,Real,Real,Real>{},
vfrac, IntVect(0),
[=] AMREX_GPU_DEVICE (int b, int i, int j, int k) noexcept
-> GpuTuple<Long,Long,Long,Real,Real,Real>
{
auto flag = ebdata.get<EBData_t::cellflag>(b,i,j,k);
Long nreg = 0, ncut = 0, ncov = 0;
Real eb_area = 0, xsec = 0;
if (flag.isRegular()) {
nreg = 1;
xsec = (i == i_mid) ? Real(1.) : Real(0.);
} else if (flag.isCovered()) {
ncov = 1;
} else {
ncut = 1;
eb_area = ebdata.get<EBData_t::bndryarea>(b,i,j,k);
xsec = (i == i_mid) ? ebdata.get<EBData_t::apx>(b,i,j,k) : Real(0.);
}
return {nreg, ncut, ncov,
ebdata.get<EBData_t::volfrac>(b,i,j,k), eb_area, xsec};
});
The volume is compared with the exact value when one is known (not for a
user-supplied parser.function or another eb2.* geometry), and for the
sphere also the EB area and the area at \(x = 0\). The program fails if a
relative error exceeds 2% at n_cell = 32 (the tolerance scales with
\(\Delta x^2\)). It writes a plotfile plt_<geometry> with the cell
type and the volume fraction, and in 3D the EB surface of
surface_geometry to eb.pvtp for ParaView.
Run (from the top of the repository, see Getting Started):
cmake --build build --target EB_Geometry -j 8
cd build/EB/Geometry
mpiexec -n 2 ./Geometry inputs
Try geometries=parser parser.function="x^2+y^2+z^2-0.25", or a finer
grid, n_cell=128.
Poisson
- Level:
★★
- Source:
- Key APIs:
EB2::SphereIF,EB2::Build,EBFArrayBoxFactory,EBData_t::centroid,EBData_t::bndrycent,MLEBABecLap,MLEBABecLap::setEBDirichlet,MLMG,EB_WriteSingleLevelPlotfile
Solve \(-\nabla^2\phi = f\) on \([0,1]^D\) outside a spherical
obstacle with the cell-centered EB solver MLEBABecLap, with Dirichlet
conditions on the domain boundary and inhomogeneous Dirichlet conditions on
the EB. The exact solution is \(\phi = \prod_d \sin(\pi x_d)\). The
program prints the max-norm error at cell centers and fails if it exceeds
\(10\,\Delta x^2\); the scheme is second order. The plotfile plt
has phi, exact, error and vfrac.
Multigrid needs the EB on coarsened grids, so EB2::Build is called with
a large max_coarsening_level. MultiFabs built with the factory carry the
EB data:
// A spherical obstacle with the fluid outside. Multigrid needs
// the EB on coarsened levels, so allow as many as possible.
EB2::SphereIF sphere(radius, {AMREX_D_DECL(0.5_rt,0.5_rt,0.5_rt)}, false);
const int required_coarsening_level = 0;
const int max_coarsening_level = 30;
EB2::Build(EB2::makeShop(sphere), geom, required_coarsening_level,
max_coarsening_level);
EBFArrayBoxFactory factory(EB2::IndexSpace::top().getLevel(geom), geom,
grids, dmap, {2,2,2}, EBSupport::full);
In cut cells, the right-hand side is evaluated at the cell centroid and the EB value at the EB centroid. Both are stored relative to the cell center in units of \(\Delta x\):
GpuArray<Real,AMREX_SPACEDIM> xc, xb;
for (int d = 0; d < AMREX_SPACEDIM; ++d) {
Real c = ebdata.get<EBData_t::centroid>(b,i,j,k,d);
Real bc = ebdata.get<EBData_t::bndrycent>(b,i,j,k,d);
xc[d] = problo[d] + (Real(iv[d])+Real(0.5)+c)*dx[d];
xb[d] = problo[d] + (Real(iv[d])+Real(0.5)+bc)*dx[d];
}
rhs_a[b](i,j,k) = fac * exact_solution(AMREX_D_DECL(xc[0],xc[1],xc[2]));
phieb_a[b](i,j,k) = exact_solution(AMREX_D_DECL(xb[0],xb[1],xb[2]));
MLEBABecLap solves \((a\alpha - b\nabla\cdot\beta\nabla)\phi =
\mathrm{rhs}\), here with \(a = 0\) and \(b = \beta = 1\). Without
setEBDirichlet or setEBHomogDirichlet, the EB condition is
homogeneous Neumann (see ABecLaplacian for the
non-EB operator):
// Operator: (a alpha - b div beta grad) phi = rhs with a = 0, b = 1,
// beta = 1, i.e., -lap(phi) = rhs.
MLEBABecLap mlebabec({geom}, {grids}, {dmap}, LPInfo(), {&factory});
// Homogeneous Dirichlet on the domain boundary. setDomainBC must be
// called first; a nullptr in setLevelBC means zero boundary values.
mlebabec.setDomainBC({AMREX_D_DECL(LinOpBCType::Dirichlet,
LinOpBCType::Dirichlet,
LinOpBCType::Dirichlet)},
{AMREX_D_DECL(LinOpBCType::Dirichlet,
LinOpBCType::Dirichlet,
LinOpBCType::Dirichlet)});
mlebabec.setLevelBC(0, nullptr);
mlebabec.setScalars(0.0_rt, 1.0_rt);
mlebabec.setACoeffs(0, 0.0_rt);
mlebabec.setBCoeffs(0, 1.0_rt);
// Inhomogeneous Dirichlet on the EB: values in phi_eb, beta = 1 there.
// (The default EB boundary condition is homogeneous Neumann.)
mlebabec.setEBDirichlet(0, phi_eb, 1.0_rt);
MLMG mlmg(mlebabec);
mlmg.setVerbose(verbose);
const Real tol_rel = std::is_same_v<Real,float> ? 1.e-5_rt : 1.e-11_rt;
const Real tol_abs = 0.0_rt;
phi.setVal(0.0_rt);
mlmg.solve({&phi}, {&rhs}, tol_rel, tol_abs);
Run (from the top of the repository, see Getting Started):
cmake --build build --target EB_Poisson -j 8
cd build/EB/Poisson
mpiexec -n 2 ./Poisson inputs
The error drops by about 4 when n_cell doubles:
MLMG: Final Iter. 12 resid, resid/resid0 = 1.220027457e-08, 4.927121892e-12
...
Max error: 0.005236049738 (tolerance 0.009765625)
PASSED
Advection
- Level:
★★★
- Source:
- Key APIs:
EBFArrayBoxFactory,EBDataArrays,EB_computeDivergence,MLEBABecLap,MLMG::getFluxes,ApplyRedistribution,EB_set_covered,EB_average_face_to_cellcenter,EB_WriteSingleLevelPlotfile
Conservative finite-volume advection of a scalar around a spherical obstacle on a periodic domain, with redistribution for small cut cells.
A uniform velocity is made discretely divergence free around the obstacle with a MAC projection solved by
MLEBABecLap. Afterwards, the sum of area fraction times velocity over the faces of every cell is zero, and nothing flows through the EB.Fluxes use second-order upwinding with limited slopes, dropping to first order next to the EB. A cut cell with volume fraction \(\kappa\) and face area fractions \(\alpha\) is updated with
\[\frac{dS}{dt} = -\frac{1}{\kappa}\sum_d \frac{(\alpha u S)_{d,\mathrm{hi}} - (\alpha u S)_{d,\mathrm{lo}}}{\Delta x_d}.\]A small cut cell would limit the time step.
ApplyRedistributionmoves part of its update to its neighbors without changing the total, so the time step is based on the full cell size. Time integration is two-stage SSP Runge-Kutta.
The projection solves \(\nabla^2\phi = \nabla\cdot\mathbf{u}\) with the
default homogeneous Neumann condition on the EB, and getFluxes returns
\(-\nabla\phi\) on face centroids:
MLEBABecLap op({geom}, {ba}, {dm}, LPInfo(), {&factory});
op.setDomainBC({AMREX_D_DECL(LinOpBCType::Periodic, LinOpBCType::Periodic,
LinOpBCType::Periodic)},
{AMREX_D_DECL(LinOpBCType::Periodic, LinOpBCType::Periodic,
LinOpBCType::Periodic)});
op.setLevelBC(0, nullptr);
op.setScalars(0.0_rt, 1.0_rt);
op.setBCoeffs(0, 1.0_rt); // homogeneous Neumann on the EB by default
MLMG mlmg(op);
mlmg.setVerbose(verbose);
const Real tol_rel = std::is_same_v<Real,float> ? 1.e-5_rt : 1.e-12_rt;
mlmg.solve({&phi}, {&divu}, tol_rel, 0.0_rt);
// The fluxes are -grad(phi) on face centroids.
Array<MultiFab,AMREX_SPACEDIM> gradphi;
for (int d = 0; d < AMREX_SPACEDIM; ++d) {
gradphi[d].define(u[d].boxArray(), dm, 1, 0, MFInfo(), factory);
}
mlmg.getFluxes({GetArrOfPtrs(gradphi)}, MLMG::Location::FaceCentroid);
for (int d = 0; d < AMREX_SPACEDIM; ++d) {
MultiFab::Add(u[d], gradphi[d], 0, 0, 1, 0);
}
The conservative divergence divides by the volume fraction:
// Conservative divergence, (sum of a u S) / (vfrac dx), with vfrac = 1
// in regular cells
{
const auto dxinv = m_geom.InvCellSizeArray();
auto const& rhs = m_dSdt.arrays();
AMREX_D_TERM(auto const& fx = m_flux[0].const_arrays();,
auto const& fy = m_flux[1].const_arrays();,
auto const& fz = m_flux[2].const_arrays();)
ParallelFor(m_dSdt, IntVect(0), ncomp,
[=] AMREX_GPU_DEVICE (int b, int i, int j, int k, int n) noexcept
{
const Real vf = ebdata.get<EBData_t::volfrac>(b,i,j,k);
if (vf > Real(0.)) {
rhs[b](i,j,k,n) = -(AMREX_D_TERM(
dxinv[0]*(fx[b](i+1,j,k,n) - fx[b](i,j,k,n)),
+ dxinv[1]*(fy[b](i,j+1,k,n) - fy[b](i,j,k,n)),
+ dxinv[2]*(fz[b](i,j,k+1,n) - fz[b](i,j,k,n)))) / vf;
} else {
rhs[b](i,j,k,n) = Real(0.);
}
});
}
m_dSdt.FillBoundary(m_geom.periodicity());
Redistribution is applied box by box where the box has cut cells:
const FabType type = flags[mfi].getType(amrex::grow(bx,nghost));
if (type == FabType::singlevalued) {
// Scratch space: weights for FluxRedist, a temporary state for
// StateRedist. It is local to the tile.
FArrayBox scratch(amrex::grow(bx,nghost), ncomp, The_Async_Arena());
scratch.setVal<RunOn::Device>(1.0_rt);
ApplyRedistribution(bx, ncomp, out, in, S.const_array(mfi),
scratch.array(), flags[mfi].const_array(),
AMREX_D_DECL(area[0]->const_array(mfi),
area[1]->const_array(mfi),
area[2]->const_array(mfi)),
vfrac.const_array(mfi),
AMREX_D_DECL(fcent[0]->const_array(mfi),
fcent[1]->const_array(mfi),
fcent[2]->const_array(mfi)),
ccent.const_array(mfi), m_bcrec.data(), m_geom,
dt, m_redistribution_type);
The state has two components, a Gaussian blob and a constant. The program checks that there are no NaNs, that the total mass (\(\sum \kappa S\,\Delta V\)) is conserved to round-off, and that the constant stays constant, which holds only if the velocity is divergence free.
Run (from the top of the repository, see Getting Started):
cmake --build build --target EB_Advection -j 8
cd build/EB/Advection
mpiexec -n 2 ./Advection inputs
Choose redistribution_type = StateRedist (default), FluxRedist
or NoRedist; in 3D, NoRedist is unstable with the default time
step. To watch the blob go around the obstacle, run with n_cell=64
stop_time=1 plot_int=10.
max |div u| before projection: 881.403378
max |div u| after projection: 1.707217148e-11
Redistribution: StateRedist, dt = 0.006040948875, initial mass = 0.007724640317
Steps: 17, time = 0.1
...
Relative change in mass: 3.817692201e-15 (tolerance 2.220446049e-13)
Max deviation of the constant: 1.914024494e-13 (tolerance 1e-08)
PASSED