Applications

Complete solvers that combine several topics, meant to be copied as the starting point of a new application. They may use a Source/ directory for the solver and Exec/<Problem>/ directories for problem setups. CNS needs AMReX built with AmrLevel support (AMReX_AMRLEVEL, on by default) and IncFlow needs the linear solvers (AMReX_LINEAR_SOLVERS and AMReX_LINEAR_SOLVERS_INCFLO, on by default). Both run in 2D and 3D only. Build only this topic with -DAMReXTutorials_TOPICS=Applications; see Getting Started.

CNS

Level:: Blueprint
Source:: Applications/CNS
Key APIs:: Amr, AmrLevel, LevelBld, StateDescriptor (desc_lst), DeriveList (derive_lst), AmrLevel::RK, FluxRegister, average_down, AMRErrorTag, GpuBndryFuncFab, ParmParse::queryAdd, ParallelFor, ParReduce

A compressible Euler solver built on amrex::Amr and amrex::AmrLevel, with adaptive mesh refinement, subcycling in time, refluxing and GPU support. It uses a second-order Godunov method with characteristic slope limiting and an approximate Riemann solver, advanced in time with Runge-Kutta (cns.rk_order = 2 or 3). It comes with two problems: a Sod shock tube (Exec/Sod) and a Rayleigh-Taylor instability with gravity (Exec/RT). At the end of a run, the code checks that the total mass changed only by round-off, if no mass can enter or leave the domain.

Compared with AdvectionAmrCore, which manages level data itself on top of AmrCore, here AmrLevel provides the state data with old and new time levels, FillPatch, checkpoint/restart and plotfile I/O.

Code structure

File	Contents
`Source/main.cpp`	Initialize, build `Amr`, time loop, mass check
`Source/CNS.H`, `Source/CNS.cpp`	The `AmrLevel` subclass, one object per level
`Source/CNS_setup.cpp`	State data, boundary conditions, derived variables, `LevelBld`
`Source/CNS_advance.cpp`	`advance()` and the flux computation
`Source/CNS_hydro_K.H`	GPU kernels: primitive variables, slopes, Riemann solver
`Source/CNS_K.H`	GPU kernels: time step, temperature
`Source/CNS_derive.H`, `Source/CNS_derive.cpp`	Derived plot variables (pressure, velocity)
`Source/CNS_index.H`, `Source/CNS_parm.H`, `Source/CNS_parm.cpp`	Component indices; physical parameters (EOS, gravity)
`Exec/<Problem>/cns_prob.H`	`ProbParm` and the initial condition `cns_initdata`
`Exec/<Problem>/cns_prob.cpp`	Default parameters (`cns_prob_defaults`) and `prob.*` parameters

Amr owns the grid hierarchy and calls the virtual functions of CNS:

Function	Called
`variableSetUp` (static)	Once, to define the state data and derived variables
`initData`	For each level at start-up
`computeInitialDt`, `computeNewDt`	On level 0, to choose `dt` for all levels
`advance`	For each level and each (sub)step
`post_timestep`	After a level has caught up with the next finer level: reflux and average down
`errorEst`	When regridding, to tag cells for refinement
`init(old)`, `init()`	When regridding, to fill new grids from old data or from the coarser level
`post_init`, `postCoarseTimeStep`	After initialization and after each coarse step
`restart`	When restarting from a checkpoint

Driver and setup

The driver builds Amr from a LevelBld and takes coarse time steps; Amr recursively advances the finer levels and writes plotfiles and checkpoints according to the amr.* parameters:

// Amr manages the hierarchy of CNS levels: it builds the initial
// grids, regrids, subcycles in time, and writes plotfiles and
// checkpoints according to the amr.* parameters.
Amr amr(getLevelBld());
amr.init(strt_time, stop_time);

const Real mass0 = dynamic_cast<CNS&>(amr.getLevel(0)).conservedTotals()[0];

while ( amr.okToContinue() &&
        (amr.levelSteps(0) < max_step || max_step < 0) &&
        (amr.cumTime() < stop_time || stop_time < Real(0.0)) )
{
    // One coarse time step, which recursively advances finer levels.
    amr.coarseTimeStep(stop_time);
}

The state is registered once in CNS::variableSetUp: one cell-centered StateData with NUM_GROW ghost cells and conservative interpolation, the component names with their BCRec boundary types and fill function, and derived quantities computed on demand, e.g., for amr.derive_plot_vars:

// One cell-centered StateData with NCONS components and NUM_GROW ghost
// cells, interpolated conservatively from coarse to fine.
desc_lst.addDescriptor(State_Type, IndexType::TheCellType(),
                       StateDescriptor::Point, NUM_GROW, NCONS,
                       &cell_cons_interp, /*extrap=*/ false,
                       /*store_in_checkpoint=*/ true);

Vector<std::string> names{"density", "xmom", "ymom", "zmom", "rho_E", "rho_e", "Temp"};
Vector<BCRec> bcs(NCONS);
for (int n = 0; n < NCONS; ++n) {
    bcs[n] = component_bc(n, phys_bc);
}

StateDescriptor::BndryFunc bndryfunc(GpuBndryFuncFab<CnsFillExtDir>(CnsFillExtDir{}));
bndryfunc.setRunOnGPU(true);  // the BC function launches GPU kernels itself

desc_lst.setComponent(State_Type, URHO, names, bcs, bndryfunc);

// Derived quantities, computed on demand for plotfiles
// (amr.derive_plot_vars) and tagging.
derive_lst.add("pressure", IndexType::TheCellType(), 1, cns_derpres,
               &DeriveRec::TheSameBox);
derive_lst.addComponent("pressure", desc_lst, State_Type, UEINT, 1);

derive_lst.add("velocity", IndexType::TheCellType(), AMREX_SPACEDIM,
               {AMREX_D_DECL("x_velocity","y_velocity","z_velocity")}, cns_dervel,
               &DeriveRec::TheSameBox);
derive_lst.addComponent("velocity", desc_lst, State_Type, URHO, 1+AMREX_SPACEDIM);

Time step and synchronization

CNS::advance calls AmrLevel::RK, which runs the Runge-Kutta stages and fills the ghost cells of each stage, including coarse/fine boundaries interpolated in space and time. compute_dSdt computes the fluxes with ParallelFor kernels and adds them to the flux registers:

// As the coarse side, this level adds its fluxes to the register owned by
// the next finer level.  As the fine side, it adds to its own register.
FluxRegister* fr_as_crse = nullptr;
if (do_reflux && level < parent->finestLevel()) {
    fr_as_crse = getLevel(level+1).flux_reg.get();
    fr_as_crse->setVal(Real(0.0));
}

FluxRegister* fr_as_fine = nullptr;
if (do_reflux && level > 0) {
    fr_as_fine = flux_reg.get();
}

// AmrLevel::RK fills the ghost cells of each stage, including
// coarse/fine boundaries interpolated in space and time.
RK(rk_order, State_Type, time, dt, iteration, ncycle,
   // Given state S, compute dSdt.  dtsub is this stage's weight times dt.
   [&] (int /*stage*/, MultiFab& dSdt, MultiFab const& S, Real /*t*/, Real dtsub)
   {
       compute_dSdt(S, dSdt, dtsub, fr_as_crse, fr_as_fine);
   },
   // After each stage, update the derived components.
   [&] (int /*stage*/, MultiFab& S) { computeTemp(S, 0); });

After the finer level has caught up, post_timestep corrects the coarse level with the fine fluxes and overwrites the covered cells with the average of the fine data, which together keep the scheme conservative:

if (level < parent->finestLevel()) {
    // Replace the coarse fluxes at the coarse/fine interface with the
    // sum of the fine fluxes, which keeps the scheme conservative.
    if (do_reflux) {
        MultiFab& S = get_new_data(State_Type);
        getLevel(level+1).flux_reg->Reflux(S, Real(1.0), 0, 0, NCONS, geom);
    }
    avgDown();
}

Making your own code

New problem: copy Exec/Sod to Exec/MyProblem. Edit ProbParm and cns_initdata in cns_prob.H, and the defaults and prob.* parameters in cns_prob.cpp. Add the problem to the list in CMakeLists.txt.
New physics: add components in CNS_index.H and their names and boundary conditions in CNS_setup.cpp, and change the kernels in CNS_hydro_K.H and compute_dSdt. Inflow boundaries need values in CnsFillExtDir and removal of the check in CNS::read_params.
New AmrLevel code: keep the structure of main.cpp, CNS_setup.cpp (LevelBld) and the AmrLevel overrides in CNS.cpp, and replace advance, the kernels and the tagging criteria. Kernels capture small parameter structs such as Parm by value; large tables belong in GPU memory, e.g., Gpu::DeviceVector.

Building and running

One executable is built per problem: targets Applications_CNS_Sod and Applications_CNS_RT, binaries CNS_Sod and CNS_RT. Both problems have an inputs file, which is read from the source tree:

cmake --build build --target Applications_CNS_Sod -j 8
cd build/Applications/CNS
mpiexec -n 2 ./CNS_Sod ../../../Applications/CNS/Exec/Sod/inputs

After copying Applications/CNS out of the repository, its CMakeLists.txt builds it on its own and fetches AMReX (or uses a local checkout given by -DFETCHCONTENT_SOURCE_DIR_AMREX):

cmake -S . -B build [-DAMReX_SPACEDIM=2] [-DAMReX_GPU_BACKEND=CUDA]
cmake --build build -j 8

With GNU Make, build in a problem directory:

cd Applications/CNS/Exec/Sod
make -j8
mpiexec -n 2 ./CNS3d.gnu.MPI.ex inputs

Without an inputs file, each problem runs to its stop_time with the defaults in cns_prob_defaults. The inputs files stop after 20 steps; remove max_step and increase amr.n_cell and amr.max_level for a real run. Boundary types (cns.lo_bc, cns.hi_bc) are 0 interior (periodic), 1 inflow (not implemented), 2 outflow, 3 symmetry, 4 slip wall and 5 no-slip wall. cns.v=2 prints the conserved totals after every coarse step.

IncFlow

Level:: Blueprint
Source:: Applications/IncFlow
Key APIs:: AmrCore, FillPatchTwoLevels, InterpFromCoarseLevel, GpuBndryFuncFab, MLPoisson, MLNodeLaplacian, MLABecLaplacian, MLMG::solve, MLMG::getFluxes, MLMG::apply, computeDivergence, average_down_faces, amrex_calc_xslope_extdir, AMRErrorTag

A small incompressible Navier-Stokes solver built on amrex::AmrCore, in the spirit of incflo but with one method for each step and no options. It solves

\[\frac{\partial u}{\partial t} + \nabla\cdot(u u) = -\nabla p + \nu \nabla^2 u, \qquad \nabla\cdot u = 0\]

with constant density, adaptive mesh refinement and periodic or no-slip wall boundaries. All levels advance together with the same time step. The projections are written directly with the MLMG linear solvers. It comes with three problems:

Exec/TaylorGreen: decaying Taylor-Green vortices on one level, checked against the exact solution.
Exec/DoubleShearLayer: two periodic shear layers that roll up into vortices (Kelvin-Helmholtz), with a fine level that follows the vorticity. The mean velocity must be conserved to round-off, also across levels.
Exec/LidDrivenCavity: no-slip walls and a moving lid, refined near the lid.

Compared with AdvectionAmrCore, which also manages its level data on top of AmrCore, the solution here is a velocity that must stay divergence free, so each step involves elliptic solves over all levels.

Code structure

File	Contents
`Source/main.cpp`	Initialize, `InitData`, `Evolve`, check the result
`Source/IncFlow.H`, `Source/IncFlow.cpp`	The `AmrCore` subclass: parameters, boundary conditions, level creation, tagging, `FillPatch`
`Source/Advance.cpp`	Time step, advection and viscous terms
`Source/Advection_K.H`	GPU kernels: slopes, face values, upwinding
`Source/Projection.cpp`	MAC and nodal projections
`Source/IO.cpp`	Plotfiles, checkpoint and restart
`Exec/<Problem>/prob.H`	The initial velocity `prob_initdata`
`Exec/<Problem>/prob.cpp`	Default parameters (`prob_add_defaults`) and the check of the result (`prob_check`)

AmrCore calls MakeNewLevelFromScratch at start-up, MakeNewLevelFromCoarse, RemakeLevel and ClearLevel when regridding, and ErrorEst to tag cells where the vorticity exceeds incflow.vort_tag. Ghost cells are filled by FillPatch from the same level, the coarser level, and the wall velocity.

Time step

Each step is a predictor-corrector method. ComputeRHS returns \(R(u) = -\nabla\cdot(u u) + \nu\nabla^2 u\) without the pressure gradient, and the nodal projection removes the gradient part of the result:

// Predictor: u* = P(u^n + dt R(u^n))
ComputeRHS(rhs_old, t_old);
for (int lev = 0; lev < nlevs; ++lev) {
    MultiFab::LinComb(vel_new[lev], 1.0_rt, vel_old[lev], 0, dt, rhs_old[lev], 0,
                      0, AMREX_SPACEDIM, 0);
}
NodalProjection(vel_new, dt);

// Corrector: u^{n+1} = P(u^n + dt/2 (R(u^n) + R(u*)))
ComputeRHS(rhs_new, t_new);
for (int lev = 0; lev < nlevs; ++lev) {
    MultiFab::LinComb(vel_new[lev], 1.0_rt, vel_old[lev], 0, 0.5_rt*dt, rhs_old[lev], 0,
                      0, AMREX_SPACEDIM, 0);
    MultiFab::Saxpy(vel_new[lev], 0.5_rt*dt, rhs_new[lev], 0, 0, AMREX_SPACEDIM, 0);
}
NodalProjection(vel_new, dt);

For the advection term, the velocity is extrapolated to both sides of each face with limited slopes, one-sided next to walls:

void face_values (amrex::IntVect const& iv, int n, int d,
                  amrex::Array4<amrex::Real const> const& q, DirInfo const& di,
                  amrex::Real& ql, amrex::Real& qr)
{
    using namespace amrex::literals;
    amrex::IntVect const im = iv - amrex::IntVect::TheDimensionVector(d);
    if (di.wall_lo && iv[d] == di.dlo) {
        ql = qr = q(im,n);
    } else if (di.wall_hi && iv[d] == di.dhi+1) {
        ql = qr = q(iv,n);
    } else {
        ql = q(im,n) + 0.5_rt * slope(im, n, d, q, di);
        qr = q(iv,n) - 0.5_rt * slope(iv, n, d, q, di);
    }
}

The upwinded normal component gives a face velocity umac. A MAC projection makes it divergence free: solve \(\nabla^2\phi = \nabla\cdot u_{mac}\) on all levels and add the fluxes \(-\nabla\phi\) from MLMG::getFluxes to umac:

if (!mac_op) {
    mac_op = std::make_unique<MLPoisson>(LevelGeom(), LevelGrids(), LevelDmap());
    mac_op->setDomainBC(LinOpBC(LinOpBCType::Neumann), LinOpBC(LinOpBCType::Neumann));
    for (int lev = 0; lev < nlevs; ++lev) {
        mac_op->setLevelBC(lev, nullptr);
    }
}

MLMG mlmg(*mac_op);
mlmg.setVerbose(mg_verbose);
mlmg.solve(GetVecOfPtrs(phi), GetVecOfConstPtrs(rhs), solver_tolerance(), 0.0_rt);

All velocity components are then upwinded with umac, and the fluxes on faces covered by a finer level are replaced by the average of the fine fluxes, so that the update is conservative. The viscous term is applied with an MLABecLaplacian, which handles the walls (Dirichlet values in the ghost cells) and coarse/fine boundaries:

if (!visc_op) {
    // Only the operator is needed, no multigrid hierarchy.
    LPInfo info;
    info.setMaxCoarseningLevel(0);
    visc_op = std::make_unique<MLABecLaplacian>(LevelGeom(), LevelGrids(), LevelDmap(),
                                                info, Vector<FabFactory<FArrayBox> const*>{},
                                                AMREX_SPACEDIM);
    visc_op->setDomainBC(LinOpBC(LinOpBCType::Dirichlet), LinOpBC(LinOpBCType::Dirichlet));
    // alpha a u - beta div(b grad(u)) with alpha = 0, beta = b = 1
    visc_op->setScalars(0.0_rt, 1.0_rt);
    for (int lev = 0; lev < nlevs; ++lev) {
        visc_op->setBCoeffs(lev, 1.0_rt);
    }
}
for (int lev = 0; lev < nlevs; ++lev) {
    // The ghost cells outside walls hold the wall velocity.
    visc_op->setLevelBC(lev, &vel[lev]);
}

MLMG mlmg(*visc_op);
mlmg.apply(GetVecOfPtrs(visc), GetVecOfPtrs(vel));

The nodal projection solves \(\nabla^2\phi = \nabla\cdot u\) with \(\phi\) on the nodes. Since the update has no pressure gradient, \(u = u_{div\,free} + \Delta t\,\nabla p\), so \(p = \phi/\Delta t\). Building a linear operator sets up its multigrid hierarchy, so all three operators are built when first needed and kept until the grids change:

if (!nodal_op) {
    // Constant sigma = 1. The empty list of factories means no embedded
    // boundaries.
    Vector<FabFactory<FArrayBox> const*> const no_eb{};
    nodal_op = std::make_unique<MLNodeLaplacian>(LevelGeom(), LevelGrids(), LevelDmap(),
                                                 LPInfo(), no_eb, 1.0_rt);
    nodal_op->setDomainBC(LinOpBC(LinOpBCType::Neumann), LinOpBC(LinOpBCType::Neumann));
}

// rhs = div(vel) on the nodes, including the coarse/fine synchronization.
// The velocity in ghost cells outside the Neumann walls is not used.
nodal_op->compRHS(GetVecOfPtrs(rhs), GetVecOfPtrs(vel), {}, {});

MLMG mlmg(*nodal_op);
mlmg.setVerbose(mg_verbose);
mlmg.solve(GetVecOfPtrs(phi), GetVecOfConstPtrs(rhs), solver_tolerance(), 0.0_rt);

Making your own code

New problem: copy Exec/TaylorGreen to Exec/MyProblem. Edit prob_initdata in prob.H, and the defaults and the check in prob.cpp. Add the problem to the list in CMakeLists.txt.
Wall velocity: xlo.velocity, yhi.velocity, … set the velocity of the walls in non-periodic directions; the normal component is zero.
More physics: a scalar can be advected like a velocity component in ComputeAdvection. Forces such as gravity are added to \(R\) in ComputeRHS. Inflow and outflow need Dirichlet values in WallVelocityFill and Dirichlet pressure boundaries in the projections.

Building and running

One executable is built per problem: targets Applications_IncFlow_<Problem>, binaries IncFlow_<Problem>:

cmake --build build --target Applications_IncFlow_TaylorGreen -j 8
cd build/Applications/IncFlow
mpiexec -n 2 ./IncFlow_TaylorGreen ../../../Applications/IncFlow/Exec/TaylorGreen/inputs

As with CNS, the directory can be copied out and built on its own with CMake, and each problem directory has a GNUmakefile. The inputs files stop after 10 steps; remove max_step for a longer run, e.g., the shear layers roll up by stop_time=1, and stop_time=20 brings the cavity close to steady state. The nodal projection needs amr.ref_ratio = 2. Restart with amr.restart=chk00005.