Base

The Base tutorials cover the core of AMReX: initialization, runtime parameters, the mesh data structures (Box, BoxArray, DistributionMapping, Geometry, MultiFab), GPU-portable kernels with ParallelFor, ghost cells and boundary conditions, reductions, GPU containers, and I/O. They need no optional AMReX components; ExternalMPI requires MPI. Build only this topic with -DAMReXTutorials_TOPICS=Base (add -DAMReX_FFT=OFF if FFTW is not installed); see Getting Started.

Start with HelloWorld, MultiFab and HeatEquation, in that order.

HelloWorld

Level:: ★
Source:: Base/HelloWorld
Key APIs:: amrex::Initialize, amrex::Finalize, amrex::Print, amrex::AllPrint

The smallest AMReX program. amrex::Initialize sets up MPI and the GPU device; amrex::Finalize releases them. AMReX objects must be destroyed before Finalize, hence the extra scope. amrex::Print writes from the I/O rank only, amrex::AllPrint from every rank.

int main (int argc, char* argv[])
{
    // Initialize MPI and the GPU device, if enabled.
    amrex::Initialize(argc, argv);
    {
        // amrex::Print only prints on the I/O process (rank 0).
        amrex::Print() << "Hello world from AMReX version "
                       << amrex::Version() << "\n";

        amrex::Print() << "  MPI ranks:      " << amrex::ParallelDescriptor::NProcs() << "\n"
                       << "  OpenMP threads: " << amrex::OpenMP::get_max_threads() << "\n"
                       << "  GPU backend:    "
#if defined(AMREX_USE_CUDA)
                       << "CUDA"
#elif defined(AMREX_USE_HIP)
                       << "HIP"
#elif defined(AMREX_USE_SYCL)
                       << "SYCL"
#else
                       << "none"
#endif
                       << "\n";

        // amrex::AllPrint prints on every process.
        amrex::AllPrint() << "  Rank " << amrex::ParallelDescriptor::MyProc()
                          << " says hello\n";
    }
    amrex::Finalize();
}

Run:

cmake --build build --target Base_HelloWorld -j 8
cd build/Base/HelloWorld
mpiexec -n 2 ./HelloWorld

ParmParse

Level:: ★
Source:: Base/ParmParse
Key APIs:: ParmParse::query, get, queryarr, queryAdd, queryWithParser, querytable, query_enum_case_insensitive, AMREX_ENUM, ParmParse::prettyPrintUnusedInputs

Reads runtime parameters from an inputs file and the command line: scalars, strings, arrays, math expressions, enums, 2D tables, TOML-like [section] headers, and unused (often misspelled) inputs. Arguments after the inputs file, e.g. nsteps=200, override the file.

query returns whether the name was found and leaves the variable unchanged otherwise; get aborts if the name is missing; queryAdd also stores the default in the database.
Expressions without spaces (dt = 1.0/64) are evaluated by query; expressions with spaces need queryWithParser. Unknown symbols are looked up in the database, here with the my_constants. prefix.
[physics] starts a table: diffusivity = 0.1 below it is physics.diffusivity.

The inputs file:

nsteps = 100
title  = "ParmParse tutorial"   # quotes keep the spaces in a string
n_cell = 64 32 16               # an array

# Math expressions without spaces are evaluated by query()/get().
max_grid_size = 2**5            # 32
dt = 1.0/64                     # 0.015625

# Unknown symbols in an expression are looked up in the database.  The
# ParmParse object is created with "my_constants" as its parser prefix, so
# `L` is found as my_constants.L.
my_constants.L  = 2.0
my_constants.pi = 3.141592653589793
dx = L/32                       # 0.0625
omega = 2 * pi / L              # with spaces: use queryWithParser

# Names of an AMREX_ENUM, matched case-insensitively here.
time_stepper = rk2

# A 2D table: the Butcher matrix of Heun's method.
butcher_a = {{ 0.0, 0.0 }
             { 1.0, 0.0 }}

# TOML-like tables: the entries below are physics.diffusivity, etc.
[physics]
diffusivity = 0.1
gravity = 0.0 -9.81
difusivity = 0.2                # typo: reported as an unused input

[output]
plot_int = 10

Reading an enum, a table, and the entries of a section:

// Enum
auto stepper = TimeStepper::RK2;
pp.query_enum_case_insensitive("time_stepper", stepper);
amrex::Print() << "time_stepper  = " << amrex::getEnumNameString(stepper) << "\n";

// 2D table
std::vector<std::vector<double>> butcher_a{{0.0, 0.0}, {1.0, 0.0}};
pp.querytable("butcher_a", butcher_a);
amrex::Print() << "butcher_a     = ";
for (auto const& row : butcher_a) { amrex::Print() << "{ " << to_string(row) << "} "; }
amrex::Print() << "\n";

// A prefix selects the entries of a [physics] table, i.e., physics.*
amrex::ParmParse pp_phys("physics");
Real diffusivity = 0.1_rt;
found = pp_phys.query("diffusivity", diffusivity);
std::vector<Real> gravity{0.0_rt, -9.81_rt};
if (std::vector<Real> g; pp_phys.queryarr("gravity", g)) { gravity = g; }

Run:

cmake --build build --target Base_ParmParse -j 8
cd build/Base/ParmParse
./ParmParse inputs nsteps=200

The output ends with the misspelled physics.difusivity = 0.2 listed as an unused input.

MultiFab

Level:: ★
Source:: Base/MultiFab
Key APIs:: Box, BoxArray::maxSize, DistributionMapping, Geometry, RealBox, MultiFab, ParallelFor(mf, ...), MFIter, TilingIfNotGPU, Array4, FillBoundary, MultiFab::sum/norm0/norm1/norm2, WriteSingleLevelPlotfile

Builds the basic data structures: a Box for the index space of the domain, a BoxArray that chops it into grids, a DistributionMapping that assigns grids to MPI ranks, and a Geometry for the physical, periodic domain. A MultiFab holds data on the grids owned by each rank, plus ghost cells.

// Box: the index space of the domain, cells 0 to n_cell-1 in each direction.
amrex::Box domain(amrex::IntVect(0), amrex::IntVect(n_cell-1));

// BoxArray: the domain chopped into boxes of at most max_grid_size cells
// in each direction.
amrex::BoxArray ba(domain);
ba.maxSize(max_grid_size);

// DistributionMapping: which MPI rank owns which box.
amrex::DistributionMapping dm(ba);

// Geometry: the physical domain [0,1]^d, periodic in all directions.
amrex::RealBox real_box({AMREX_D_DECL(0.0_rt, 0.0_rt, 0.0_rt)},
                        {AMREX_D_DECL(1.0_rt, 1.0_rt, 1.0_rt)});
amrex::Geometry geom(domain, real_box, amrex::CoordSys::cartesian,
                     {AMREX_D_DECL(1, 1, 1)});

The same function is written into two components in two ways. ParallelFor(mf, ...) with mf.arrays() runs one kernel over all local boxes; it is the simplest, and on GPUs it is efficient for many small boxes. An MFIter loop gives one Array4 per box (per tile on CPU, with TilingIfNotGPU() and OpenMP); use it when work is organized per box, e.g. with temporary FArrayBoxes or several kernels on the same box.

auto const& ma = mf.arrays();
amrex::ParallelFor(mf, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
{
    ma[b](i,j,k,0) = phi_value(amrex::IntVect(AMREX_D_DECL(i,j,k)), prob_lo, dx);
});

for (amrex::MFIter mfi(mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
{
    amrex::Box const& bx = mfi.tilebox();
    amrex::Array4<Real> const& a = mf.array(mfi);
    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
    {
        a(i,j,k,1) = phi_value(amrex::IntVect(AMREX_D_DECL(i,j,k)), prob_lo, dx);
    });
}

// Fill ghost cells from neighboring boxes, including periodic images.
mf.FillBoundary(geom.periodicity());

The MFIter loop is preceded by #pragma omp parallel if (amrex::Gpu::notInLaunchRegion()) so that tiles are distributed over OpenMP threads in CPU builds. FillBoundary fills ghost cells from neighboring boxes and periodic images; Communication covers it and the other communication functions in more detail. The program then checks valid and ghost cells against the exact function and writes the plotfile plt_multifab.

Run:

cmake --build build --target Base_MultiFab -j 8
cd build/Base/MultiFab
mpiexec -n 2 ./MultiFab inputs

HeatEquation

Level:: ★
Source:: Base/HeatEquation
Key APIs:: MultiFab, Geometry, FillBoundary, ParallelFor(mf, ...), const_arrays, WriteSingleLevelPlotfile

Solves \(\partial\phi/\partial t = D\nabla^2\phi\) with forward Euler on the periodic unit domain. The initial condition is \(\phi = 1 + \prod_d \sin(2\pi x_d)\); its sine mode decays as \(\exp(-4\pi^2 d D t)\) in \(d\) dimensions, and the program compares the result with this exact solution. The time step is \(\Delta t = \mathrm{cfl}\,\Delta x^2/(2 d D)\).

Each step fills the ghost cells and applies the standard second-order Laplacian:

// Fill ghost cells, here all from periodic images and neighbors.
phi_old.FillBoundary(geom.periodicity());

// phi_new = phi_old + dt * D * Laplacian(phi_old)
auto const& pold = phi_old.const_arrays();
auto const& pnew = phi_new.arrays();
Real const fac = dt_step * diffusion_coef;
amrex::ParallelFor(phi_new, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
{
    auto const& po = pold[b];
    pnew[b](i,j,k) = po(i,j,k) + fac * (
          AMREX_D_TERM((po(i+1,j,k) - 2.0_rt*po(i,j,k) + po(i-1,j,k)) * dxinv[0]*dxinv[0],
                     + (po(i,j+1,k) - 2.0_rt*po(i,j,k) + po(i,j-1,k)) * dxinv[1]*dxinv[1],
                     + (po(i,j,k+1) - 2.0_rt*po(i,j,k) + po(i,j,k-1)) * dxinv[2]*dxinv[2]));
});

Run:

cmake --build build --target Base_HeatEquation -j 8
cd build/Base/HeatEquation
mpiexec -n 2 ./HeatEquation inputs

Add plot_int=10 to write a plotfile every 10 steps, and e.g. n_cell=128 stop_time=0.05 for a larger run. To see where the time goes, build with TinyProfiler; see Profiling with TinyProfiler.

Communication

Level:: ★★
Source:: Base/Communication
Key APIs:: FillBoundary, FillBoundary_nowait, FillBoundary_finish, ParallelCopy, ParallelAdd, SumBoundary, Geometry::periodicity

The three basic ways to exchange MultiFab data between boxes and MPI ranks. Each valid cell holds a unique number, so every result is checked exactly.

FillBoundary copies valid data into the ghost cells of neighboring boxes, and with geom.periodicity() also across periodic boundaries. It can fill fewer ghost cells than the MultiFab has. The non-blocking version returns after starting the communication, so that work that does not need the ghost cells can be done before FillBoundary_finish:

a.FillBoundary_nowait(geom.periodicity());

// A BoxArray shifted by n_cell/4 and chopped differently, so that it
// partly lies outside the domain, with its own DistributionMapping.
amrex::BoxArray ba_dst(amrex::shift(domain, amrex::IntVect(n_cell/4)));
ba_dst.maxSize(dst_max_grid_size);
amrex::DistributionMapping dm_dst(ba_dst);
amrex::MultiFab dst(ba_dst, dm_dst, 1, 1);
dst.setVal(-1.0_rt);

a.FillBoundary_finish();

ParallelCopy copies data between MultiFabs with different BoxArrays and DistributionMappings, e.g. after regridding. The number of ghost cells of the source and the destination can be given, here 0 and 1. The destination BoxArray above is shifted and partly outside the domain; with geom.periodicity(), those cells are filled from their periodic images. ParallelAdd adds instead of copying:

dst.ParallelCopy(a, 0, 0, 1, amrex::IntVect(0), amrex::IntVect(1), geom.periodicity());
report("ParallelCopy", max_error(dst, 1, id));

// ParallelAdd adds instead of copying.
dst.ParallelAdd(a, 0, 0, 1, amrex::IntVect(0), amrex::IntVect(1), geom.periodicity());

SumBoundary is the reverse of FillBoundary: values in ghost cells are added to the valid cells they overlap, as in deposition, e.g. of particles near a box boundary. Here every valid cell deposits 1 into the \(3^d\) cells around it, itself included, some of them ghost cells, and after SumBoundary every valid cell holds \(3^d\). The ghost cells are not updated; call FillBoundary afterwards if they are needed.

Run:

cmake --build build --target Base_Communication -j 8
cd build/Base/Communication
mpiexec -n 2 ./Communication inputs

BoundaryConditions

Level:: ★★
Source:: Base/BoundaryConditions
Key APIs:: BCRec, BCType::ext_dir/foextrap/reflect_even/int_dir, PhysBCFunct, GpuBndryFuncFab, FillBoundary, MFIter, ParallelFor(box, [=] (IntVect const&) ...)

Solves the heat equation on a non-periodic domain: Dirichlet values at \(x = 0\) and \(x = 1\), zero gradient in \(y\), periodic in \(z\). The solution \(\phi = \phi_{lo} + (\phi_{hi} - \phi_{lo})x + e^{-\pi^2 D t}\sin(\pi x)\) relaxes to the linear steady state and is checked against this exact solution.

Boundary conditions are given per component as a BCRec. Periodic directions use int_dir and must match the Geometry.

        amrex::Vector<amrex::BCRec> bcs(1);
        bcs[0].setLo(0, amrex::BCType::ext_dir);
        bcs[0].setHi(0, amrex::BCType::ext_dir);
#if (AMREX_SPACEDIM >= 2)
        bcs[0].setLo(1, amrex::BCType::reflect_even);
        bcs[0].setHi(1, amrex::BCType::foextrap);
#endif
#if (AMREX_SPACEDIM == 3)
        bcs[0].setLo(2, amrex::BCType::int_dir);
        bcs[0].setHi(2, amrex::BCType::int_dir);
#endif
        amrex::Array<int,AMREX_SPACEDIM> is_periodic{};
        for (int d = 0; d < AMREX_SPACEDIM; ++d) {
            is_periodic[d] = (bcs[0].lo(d) == amrex::BCType::int_dir) ? 1 : 0;
        }

Ghost cells are filled in two steps. FillBoundary fills those that overlap other boxes, including periodic images. PhysBCFunct fills those outside the domain: AMReX handles foextrap, reflect_even, reflect_odd and hoextrap, and calls the user functor for every such cell to set the ext_dir values.

struct DirichletFill
{
    amrex::GpuArray<Real,AMREX_SPACEDIM> phi_lo;
    amrex::GpuArray<Real,AMREX_SPACEDIM> phi_hi;

    AMREX_GPU_DEVICE
    void operator() (amrex::IntVect const& iv, amrex::Array4<Real> const& dest,
                     int dcomp, int numcomp, amrex::GeometryData const& geom,
                     Real /*time*/, amrex::BCRec const* bcr, int bcomp,
                     int /*orig_comp*/) const
    {
        amrex::Box const& domain = geom.Domain();
        for (int n = 0; n < numcomp; ++n) {
            amrex::BCRec const& bc = bcr[bcomp+n];
            for (int d = 0; d < AMREX_SPACEDIM; ++d) {
                if (bc.lo(d) == amrex::BCType::ext_dir && iv[d] < domain.smallEnd(d)) {
                    dest(iv, dcomp+n) = phi_lo[d];
                } else if (bc.hi(d) == amrex::BCType::ext_dir && iv[d] > domain.bigEnd(d)) {
                    dest(iv, dcomp+n) = phi_hi[d];
                }
            }
        }
    }
};

DirichletFill dirichlet_fill{};
dirichlet_fill.phi_lo[0] = phi_lo;
dirichlet_fill.phi_hi[0] = phi_hi;
amrex::PhysBCFunct<amrex::GpuBndryFuncFab<DirichletFill>> physbc
    (geom, bcs, amrex::GpuBndryFuncFab<DirichletFill>(dirichlet_fill));

phi_old.FillBoundary(geom.periodicity());
physbc(phi_old, 0, 1, phi_old.nGrowVect(), time, 0);

For ext_dir, the ghost cell holds the value on the boundary face, so the Laplacian uses a half-cell distance next to that boundary.

Run:

cmake --build build --target Base_BoundaryConditions -j 8
cd build/Base/BoundaryConditions
mpiexec -n 2 ./BoundaryConditions inputs

Try phi_lo=2 phi_hi=-1, or stop_time=0.5 plot_int=100 to follow the approach to the steady state.

CheckpointRestart

Level:: ★★
Source:: Base/CheckpointRestart
Key APIs:: PreBuildDirectorHierarchy, MultiFabFileFullPrefix, VisMF::Write, VisMF::Read, BoxArray::writeOn/readFrom, ParallelDescriptor::ReadAndBcastFile, PlotFileData, ParallelCopy

Writes checkpoints of a heat equation run and restarts from them. A checkpoint is a directory with a Header (step, time and BoxArray), written by the I/O rank, and the MultiFab data, written by every rank with VisMF.

void WriteCheckpoint (std::string const& name, int step, Real time, amrex::MultiFab const& phi)
{
    // Create name/ and name/Level_0/ (on the I/O rank, then barrier).
    amrex::PreBuildDirectorHierarchy(name, "Level_", 1, true);

    if (amrex::ParallelDescriptor::IOProcessor()) {
        std::ofstream header(name + "/Header");
        if (!header.good()) { amrex::FileOpenFailed(name + "/Header"); }
        header.precision(17);
        header << "HeatEquation checkpoint\n" << step << "\n" << time << "\n";
        phi.boxArray().writeOn(header);
        header << "\n";
    }

    // Every rank writes the data it owns.
    amrex::VisMF::Write(phi, amrex::MultiFabFileFullPrefix(0, name, "Level_", "phi"));
}

On restart, the BoxArray comes from the checkpoint and a new DistributionMapping is made, so a run may restart on a different number of MPI ranks.

void ReadCheckpoint (std::string const& name, int& step, Real& time, amrex::MultiFab& phi,
                     int nghost)
{
    // Rank 0 reads the header and broadcasts it.
    amrex::Vector<char> buffer;
    amrex::ParallelDescriptor::ReadAndBcastFile(name + "/Header", buffer);
    std::istringstream is(buffer.dataPtr());

    std::string title;
    std::getline(is, title);
    is >> step >> time;
    is.ignore(std::numeric_limits<std::streamsize>::max(), '\n');

    amrex::BoxArray ba;
    ba.readFrom(is);

    // A new DistributionMapping; the run may use a different number of ranks.
    amrex::DistributionMapping dm(ba);
    phi.define(ba, dm, 1, nghost);
    amrex::VisMF::Read(phi, amrex::MultiFabFileFullPrefix(0, name, "Level_", "phi"));
}

With test_restart=1, the program restarts from the first checkpoint, reruns the remaining steps, and checks that the result is bitwise identical to the uninterrupted run. It also reads the final plotfile back with PlotFileData; the data come with their own DistributionMapping, so ParallelCopy moves them into the layout of another MultiFab.

amrex::PlotFileData pf(plotfile);
amrex::MultiFab phi_pf(phi.boxArray(), phi.DistributionMap(), 1, 0);
phi_pf.ParallelCopy(pf.get(0)); // all variables of level 0, here only "phi"

Run:

cmake --build build --target Base_CheckpointRestart -j 8
cd build/Base/CheckpointRestart
mpiexec -n 2 ./CheckpointRestart inputs
mpiexec -n 1 ./CheckpointRestart inputs restart=chk00010
mpiexec -n 2 ./CheckpointRestart inputs test_restart=1

Reductions

Level:: ★★
Source:: Base/Reductions
Key APIs:: ParReduce, ReduceOps, ReduceData, Reducer, ReduceOpSum/Min/Max, MultiFab::sum/min/max/norm0/norm1/norm2, Reduce::Sum, ParallelAllReduce, ParallelDescriptor::ReduceRealSum

Computes sums, minima, maxima and counts over a MultiFab and a device vector, on CPU and GPU, and combines them across MPI ranks. The data, \(\phi = i + j + k\), has simple exact answers that the program checks.

MultiFab member functions (sum, norm0, …) compute common norms and already include the MPI reduction.
ParReduce computes several custom reductions over a MultiFab in one pass.
ReduceOps/ReduceData, or the more compact Reducer, work inside MFIter loops or over arbitrary boxes and index ranges.
Reduce::Sum/Min/Max reduce a plain array, e.g. a Gpu::DeviceVector.

All except the MultiFab member functions are local to the MPI rank and are followed by ParallelAllReduce or ParallelDescriptor::Reduce*.

ParReduce with a sum, a minimum and a maximum, followed by the MPI reductions:

auto const& a = mf.const_arrays();
auto r = amrex::ParReduce(
    amrex::TypeList<amrex::ReduceOpSum, amrex::ReduceOpMin, amrex::ReduceOpMax>{},
    amrex::TypeList<Real, Real, Real>{},
    mf, amrex::IntVect(0),
    [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
        -> amrex::GpuTuple<Real, Real, Real>
    {
        Real v = a[b](i,j,k);
        return {v, v, v};
    });
Real sum = amrex::get<0>(r);
Real vmin = amrex::get<1>(r);
Real vmax = amrex::get<2>(r);
amrex::ParallelAllReduce::Sum(sum, amrex::ParallelDescriptor::Communicator());
amrex::ParallelAllReduce::Min(vmin, amrex::ParallelDescriptor::Communicator());
amrex::ParallelAllReduce::Max(vmax, amrex::ParallelDescriptor::Communicator());

Reducer with two sums of different types:

amrex::Reducer<amrex::TypeList<amrex::ReduceOpSum, amrex::ReduceOpSum>,
               amrex::TypeList<Real, Long>> reducer;
using Result_t = decltype(reducer)::Result_t;
auto const& a = mf.const_arrays();
reducer.eval(mf, amrex::IntVect(0),
[=] AMREX_GPU_DEVICE (int b, int i, int j, int k) -> Result_t
{
    Real v = a[b](i,j,k);
    Long even = (Long(v) % 2 == 0) ? 1 : 0;
    return {v*v, even};
});
Result_t r = reducer.getResult();

Run:

cmake --build build --target Base_Reductions -j 8
cd build/Base/Reductions
mpiexec -n 2 ./Reductions inputs

GpuContainers

Level:: ★★
Source:: Base/GpuContainers
Key APIs:: Gpu::DeviceVector, Gpu::PinnedVector, Gpu::copyAsync, Gpu::streamSynchronize, ParallelForRNG, amrex::Random, Scan::PrefixSum, Scan::ExclusiveSum, The_Arena, The_Async_Arena, The_Pinned_Arena

Fills a Gpu::DeviceVector with random numbers, keeps the values above a threshold with a prefix-sum based stream compaction, copies the results to pinned host memory, and checks them against the same operations on the host. The code is the same for CPU and GPU builds.

Scan::PrefixSum takes one function that produces each input and one that consumes each partial sum, so the compaction needs no temporary arrays:

// Uniform random numbers in [0,1); each GPU thread gets its own engine.
amrex::ParallelForRNG(n, [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine)
{
    pv[i] = amrex::Random(engine);
});

// Stream compaction: keep the values above threshold, in order. The
// exclusive prefix sum of the keep flags is each kept value's output
// position, and the total is the number of kept values.
amrex::Gpu::DeviceVector<Real> kept(n);
Real* pk = kept.data();
int const nkept = amrex::Scan::PrefixSum<int>(n,
    [=] AMREX_GPU_DEVICE (int i) -> int { return pv[i] > threshold; },
    [=] AMREX_GPU_DEVICE (int i, int const& pos) {
        if (pv[i] > threshold) { pk[pos] = pv[i]; }
    },
    amrex::Scan::Type::exclusive, amrex::Scan::retSum);
kept.resize(nkept);

DeviceVector memory comes from The_Arena() (device memory on GPU), and PinnedVector from The_Pinned_Arena() (page-locked host memory, needed for truly asynchronous copies). copyAsync may return before the copy is done, so synchronize the stream before using the data:

// Copy to the host. Pinned (page-locked) host memory, from
// The_Pinned_Arena(), makes asynchronous copies fast. The copy is
// asynchronous on GPU, so wait for the stream before using the data.
amrex::Gpu::PinnedVector<Real> h_values(n);
amrex::Gpu::PinnedVector<Real> h_kept(nkept);
amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, values.begin(), values.end(), h_values.begin());
amrex::Gpu::copyAsync(amrex::Gpu::deviceToHost, kept.begin(), kept.end(), h_kept.begin());
amrex::Gpu::streamSynchronize();

The_Async_Arena() is for temporaries: on GPU, memory freed there is reused only after the work already queued on the stream is done.

Run:

cmake --build build --target Base_GpuContainers -j 8
cd build/Base/GpuContainers
./GpuContainers inputs

Parser

Level:: ★★
Source:: Base/Parser
Key APIs:: Parser, IParser, setConstant, registerVariables, compile<N>, ParmParse::makeParser

Evaluates functions given as strings at runtime, here an initial condition and a tagging criterion from the inputs file, inside GPU kernels. amrex::Parser handles floating-point expressions and amrex::IParser integer ones. Expressions support + - * / ^ **, comparisons, and, or, if(cond, a, b), min, max and the usual math functions; IParser also has floor division //.

A parser is built on the host: set the constants, register the variables in the order they will be passed, and compile<N>() with N variables. ParmParse::makeParser instead resolves unknown symbols from the ParmParse database.

// 1. Parser built by hand: set the constants, then the variables in
//    the order they are passed to the compiled executor.
amrex::Parser parser(function);
parser.setConstant("pi", amrex::Math::pi<double>());
parser.setConstant("A", A);
parser.registerVariables({"x", "y", "z", "t"});
auto const f1 = parser.compile<4>(); // usable on host and device

// 2. ParmParse::makeParser looks up unknown symbols in the database,
//    here with the "my_constants" prefix.
amrex::ParmParse pp("", "my_constants");
amrex::Parser parser2 = pp.makeParser(function, {"x", "y", "z", "t"});
auto const f2 = parser2.compile<4>();

// 3. IParser for integer expressions, here of cell indices.
amrex::IParser iparser(tag);
iparser.setConstant("n", n_cell);
iparser.registerVariables({"i", "j", "k"});
auto const ftag = iparser.compile<3>();

The returned executor is a small object captured by value in device lambdas; the Parser must stay alive while it is used.

amrex::ParallelFor(phi, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
{
    amrex::IntVect const iv(AMREX_D_DECL(i,j,k));
    amrex::GpuArray<Real,3> x{0.0_rt, 0.0_rt, 0.0_rt}; // y = z = 0 in lower dimensions
    for (int d = 0; d < AMREX_SPACEDIM; ++d) {
        x[d] = prob_lo[d] + (iv[d] + 0.5_rt) * dx[d];
    }
    pa[b](i,j,k) = Real(f1(x[0], x[1], x[2], time));
    ta[b](i,j,k) = int(ftag(i, j, k));

Run:

cmake --build build --target Base_Parser -j 8
cd build/Base/Parser
./Parser inputs
./Parser inputs ic.function="x*y*exp(-t)"

CompileTimeOptions

Level:: ★★★
Source:: Base/CompileTimeOptions
Key APIs:: ParallelFor(TypeList<CompileTimeOptions<...>...>, {...}, box, f), CompileTimeOptions, if constexpr, AMREX_ENUM, BL_PROFILE

Turns runtime options into compile-time constants inside kernels. A kernel with two runtime options, the order of a Laplacian stencil and whether to add an absorption term, is written once with runtime branches and once with CompileTimeOptions. AMReX compiles a specialized kernel for every combination and launches the one that matches the runtime values. The program checks that both versions agree and times them with BL_PROFILE (see Profiling with TinyProfiler).

amrex::ParallelFor(
    amrex::TypeList<amrex::CompileTimeOptions<int(Order::second), int(Order::fourth)>,
                    amrex::CompileTimeOptions<0, 1>>{},
    {int(order), int(absorption)},
    bx, [=] AMREX_GPU_DEVICE (int i, int j, int k, auto order_ct, auto absorption_ct)
{
    amrex::IntVect const iv(AMREX_D_DECL(i,j,k));
    constexpr int ord = (order_ct.value == int(Order::second)) ? 2 : 4;
    Real r = laplacian<ord>(iv, p, dxinv);
    // A plain `if` on a compile-time constant is optimized out
    // too. (nvcc does not allow `if constexpr` blocks to be
    // the first to use a captured variable, here sigma.)
    if (absorption_ct.value) { r -= sigma * p(iv); }
    out(i,j,k) = r;
});

The kernel receives one extra argument per option, of type std::integral_constant<int, value>, so opt.value can be used in if constexpr and as a template argument. Specialized kernels can drop unused code and use fewer registers, which matters most for large GPU kernels; the cost is compile time and binary size, which grow with the number of combinations. AnyCTO does the same for other launch functions such as ParallelForRNG.

Run:

cmake --build build --target Base_CompileTimeOptions -j 8
cd build/Base/CompileTimeOptions
./CompileTimeOptions inputs

For timings, build with TinyProfiler on and use e.g. n_cell=256 nrepeat=20 max_grid_size=128 and an optimized GPU build.

TagParallelFor

Level:: ★★★
Source:: Base/TagParallelFor
Key APIs:: ParallelFor(Vector<Tag>, ncomp, f), ParallelFor(TagVector<Tag>, ncomp, f), TagVector, Array4BoxOrientationTag, Orientation, OrientationIter, adjCell, BL_PROFILE

Runs many small kernels as one launch. Filling the ghost cells outside the domain boundary touches a thin box per domain face of every box at the boundary: many small boxes of different sizes, each too small to fill a GPU. The program fills them with a reflecting boundary condition in three ways, checks that the results agree, and times them with BL_PROFILE (see Profiling with TinyProfiler): one ParallelFor per box and face inside an MFIter loop, the baseline, not recommended on GPU for many small boxes since no launch fills the device; one ParallelFor over a Vector of tags, which packs all the boxes into a single launch; and the same with a TagVector built once and reused.

A tag is a small struct with a box() member (size() for 1D work) that carries what the kernel needs, here the Array4, the thin box and the face. The tags are built on the host:

// One tag per (local box, domain face it touches): the Array4 of the box,
// the thin box of ghost cells outside the face, and the face.
amrex::Vector<amrex::Array4BoxOrientationTag<Real>>
make_boundary_tags (amrex::MultiFab& mf, amrex::Box const& domain, int nghost)
{
    amrex::Vector<amrex::Array4BoxOrientationTag<Real>> tags;
    for (amrex::MFIter mfi(mf); mfi.isValid(); ++mfi) {
        amrex::Box const& vbx = mfi.validbox();
        for (amrex::OrientationIter oit; oit; ++oit) {
            amrex::Orientation const face = oit();
            if (touches(vbx, domain, face)) {
                tags.push_back({mf.array(mfi), amrex::adjCell(vbx, face, nghost), face});
            }
        }
    }
    return tags;
}

The kernel receives (i, j, k, n, tag). ParallelFor over a Vector of tags copies the tags to the device and synchronizes the stream on every call, so for repeated use the TagVector is built once:

// 3. The TagVector is built once and reused; it only needs to be
//    rebuilt when the boxes or the MultiFab change.
amrex::TagVector<amrex::Array4BoxOrientationTag<Real>> tv(
    make_boundary_tags(phi_cached, domain, nghost));
for (int irep = 0; irep < nrepeat; ++irep) {
    BL_PROFILE("tags_cached");
    amrex::ParallelFor(tv, ncomp,
        [=] AMREX_GPU_DEVICE (int i, int j, int k, int n,
                              amrex::Array4BoxOrientationTag<Real> const& tag)
    {
        reflect(tag.fab, i, j, k, n, tag.face, dlo, dhi);
    });
    amrex::Gpu::streamSynchronize();
}

On CPU, the launch is a loop over the tags that opens its own OpenMP parallel region, so it must not be called inside one. AMReX provides Array4BoxTag, Array4BoxValTag, Array4BoxOrientationTag, Array4PairTag, Array4CopyTag and VectorTag, and any struct with box() or size() works. Use ParallelFor(mf, ...) (MultiFab) when the work is the same on every box, and tags for irregular sets of boxes: physical boundary conditions, corrections at coarse/fine faces, operations on some boxes only.

Run:

cmake --build build --target Base_TagParallelFor -j 8
cd build/Base/TagParallelFor
mpiexec -n 2 ./TagParallelFor inputs

For timings, build with TinyProfiler on and an optimized GPU build, and use e.g. n_cell=256 nrepeat=20 tiny_profiler.print_threshold=0; the difference between the methods grows with the number of boxes.

SharedMemory

Level:: ★★★
Source:: Base/SharedMemory
Key APIs:: amrex::launch<nthreads>(nblocks, stream, f), amrex::launch<nthreads>(nblocks, shared_mem_bytes, stream, f), __shared__ (CUDA, HIP), Gpu::Handler (SYCL), Gpu::Atomic::AddNoRet, HostDevice::Atomic::Add, BoxIndexer, Gpu::Device::maxBlocksPerLaunch, AMREX_USE_GPU, AMREX_USE_SYCL

A GPU kernel written by hand, with shared memory, and a separate CPU implementation. ParallelFor is portable and the right choice almost always; a raw launch is for the rare kernel whose threads must cooperate within a block. The example bins the cells by the value of one field and sums a signed integer field over each bin. Integer sums are exact, so the two methods must agree to the bit, and they need a signed 64-bit accumulator. With ParallelFor and one atomic add per cell, all threads compete for the same few addresses in global memory:

void binned_sum_atomic (amrex::MultiFab const& phi, amrex::iMultiFab const& q,
                        int nbins, Long* sums)
{
    auto const& pa = phi.const_arrays();
    auto const& qa = q.const_arrays();
    amrex::ParallelFor(phi, [=] AMREX_GPU_DEVICE (int b, int i, int j, int k)
    {
        amrex::HostDevice::Atomic::Add(sums + bin_of(pa[b](i,j,k), nbins), Long(qa[b](i,j,k)));
    });
    amrex::Gpu::streamSynchronize();
}

Atomic adds to the same address are serialized, so the fewer the bins, the slower this version. The hand-written kernel gives each thread block private int sums in shared memory, where contention is cheap, and adds them to the global Long sums once per block, so the global atomics number one per bin and block instead of one per cell. amrex::launch runs the kernel on nblocks blocks of nthreads threads; the threads of a block walk the cells of the box in a grid-stride loop, with BoxIndexer converting the linear index to (i,j,k). With CUDA and HIP the body uses threadIdx, blockIdx, __syncthreads() and a static __shared__ array:

amrex::launch<nthreads>(nblocks, amrex::Gpu::gpuStream(),
[=] AMREX_GPU_DEVICE () noexcept
{
    __shared__ int sh[max_nbins];
    int const tid = int(threadIdx.x);
    int const bdim = int(blockDim.x);
    for (int b = tid; b < nbins; b += bdim) { sh[b] = 0; }
    __syncthreads();
    auto const stride = std::uint64_t(bdim) * gridDim.x;
    for (auto icell = std::uint64_t(bdim) * blockIdx.x + tid; icell < ncells; icell += stride) {
        auto const [i, j, k] = indexer(icell);
        amrex::Gpu::Atomic::AddNoRet(sh + bin_of(pa(i,j,k), nbins), qa(i,j,k));
    }
    __syncthreads();
    for (int b = tid; b < nbins; b += bdim) {
        if (sh[b] != 0) { amrex::Gpu::Atomic::AddNoRet(sums + b, Long(sh[b])); }
    }
});

With SYCL, the shared memory is dynamic with its size passed to amrex::launch, and that overload hands the kernel a Gpu::Handler whose threadIdx(), blockDim(), gridDim(), globalIdx(), syncThreads() and sharedMemory() play these roles; atomics on the shared memory take the local_space address space; see the AMREX_USE_SYCL branch in the source. There is no CPU version of this launch. The CPU path, selected with #ifdef AMREX_USE_GPU, is an MFIter loop with a vector of private sums per OpenMP thread, merged after the parallel region; its cell loop is amrex::Loop, not ParallelFor, because it scatters into the sums:

        // Private sums per OpenMP thread, merged after the parallel region.
        int const nthreads = amrex::OpenMP::get_max_threads();
        amrex::Vector<amrex::Vector<Long>> local(nthreads, amrex::Vector<Long>(nbins, 0));
#ifdef AMREX_USE_OMP
#pragma omp parallel
#endif
        {
            auto& mine = local[amrex::OpenMP::get_thread_num()];
            for (amrex::MFIter mfi(phi, true); mfi.isValid(); ++mfi) {
                auto const& pa = phi.const_array(mfi);
                auto const& qa = q.const_array(mfi);
                amrex::Loop(mfi.tilebox(), [&] (int i, int j, int k)
                {
                    mine[bin_of(pa(i,j,k), nbins)] += qa(i,j,k);
                });
            }
        }
        for (auto const& l : local) {
            for (int b = 0; b < nbins; ++b) { sums[b] += l[b]; }
        } // end of merge

The program checks that both results agree and times them with BL_PROFILE (see Profiling with TinyProfiler).

Run:

cmake --build build --target Base_SharedMemory -j 8
cd build/Base/SharedMemory
mpiexec -n 2 ./SharedMemory inputs

For timings, build with TinyProfiler on and use a GPU build with e.g. n_cell=128 max_grid_size=128 nbins=64 nrepeat=5 tiny_profiler.print_threshold=0. The fewer the bins, the more the global atomics collide, and the larger the speedup of the shared-memory version.

ExternalMPI

Level:: ★★
Source:: Base/ExternalMPI
Key APIs:: amrex::Initialize(argc, argv, build_parm_parse, comm, func), ParmParse::add, ParmParse::addfile, ParallelDescriptor::NProcs

Uses AMReX inside an application that manages MPI itself. The application calls MPI_Init, splits MPI_COMM_WORLD into even and odd ranks, and initializes AMReX on each group’s communicator. AMReX uses a duplicate of the given communicator and does not finalize MPI; the application calls MPI_Finalize after amrex::Finalize.

With build_parm_parse = false, AMReX ignores argc/argv. The callback runs before AMReX reads its own parameters, so it is the place to set them in code.

auto add_parameters = [&] ()
{
    amrex::ParmParse pp_amrex("amrex");
    pp_amrex.add("verbose", 0);

    if (!inputs_file.empty()) {
        amrex::ParmParse::addfile(inputs_file);
    }

    amrex::ParmParse pp;
    pp.add("group", group);
};

// build_parm_parse = false: do not parse argc/argv into ParmParse.
// AMReX uses a duplicate of group_comm.
amrex::Initialize(argc, argv, false, group_comm, add_parameters);

Run (requires AMReX_MPI=ON):

cmake --build build --target Base_ExternalMPI -j 8
cd build/Base/ExternalMPI
mpiexec -n 4 ./ExternalMPI inputs