AMReX_SPACE.H)AMReX_ParallelDescriptor.H)AMReX_IntVect.H)AMReX_IndexType.H)AMReX_Box.H)AMReX_RealBox.H)AMReX_Geometry.H)AMReX_BoxArray.H)AMReX_DistributionMapping.H)AMReX is a public software framework designed for building massively parallel block-structured adaptive mesh refinement applications. Here I keep track of the basic syntax or classes used in AMReX for my personal notes.
Most are copied and paste from AMReX documentaion .
Typing them down helps me go through the documentation and filter out unnecessary information for my personal preference.
AMReX_SPACE.H)For a given problem, AMReX stores the dimensionality
of the problem via macro: AMReX_SPACEDIM.
AMReX_D_EXPR(a,b,c): evaluates the expression a,b,c based on AMReX_SPACEDIM, i.e. evaluates expression a only for DIM=1, a and b for DIM=2. Used to evaluate different expressions1AMREX_D_EXPR(vect[0] *= s, vect[1] *= s, vect[2] *= s);
AMREX_D_DECL(a,b,c): expands to a comma-separated list of, 1, 2,
or 3 arguments based on AMReX_SPACEDIM. Mainly used to write portable
function calls and declarations.1return IntVect(AMREX_D_DECL(p[0] + s, p[1] + s, p[2] + s));
IntVect(p[0]+s) or IntVect(p[0]+s, p[1]+s), etc…AMREX_D_TERM(a,b,c): expands to a whitespace-separated list of, 1, 2,
or 3 arguments based on AMReX_SPACEDIM.1vol = AMREX_D_TERM(dx, *dy, *dz);
vol=dx or vol=dx*dy or vol=dx*dy*dz.AMReX_D_PICK(a,b,c): expands to a single result equal to the
1st, 2nd, or 3rd arguments of the call based on AMReX_SPACEDIM.1maxsize = AMREX_D_PICK(1024, 128, 32);
AMReX_ParallelDescriptor.H)Mainly used functions are
1// Return the rank
2int myproc = ParallelDescriptor::MyProc();
3
4// Return the number of processes
5int nprocs = ParallelDescriptor::NProcs();
6
7if (ParallelDescriptor::IOProcessor()) {
8 // Only the I/O process executes this
9 // Mainly used to do printing
10 std::cout << "Some info..." << std::endl;
11}
12
13int ioproc = ParallelDescriptor::IOProcessorNumber(); // I/O rank
14
15// Wait until all MPI ranks have reached here.
16// Ensure all ranks have finished some operation before continuing.
17ParallelDescriptor::Barrier();
18
19// Broadcast 100 ints from the I/O Processor
20Vector<int> a(100);
21ParallelDescriptor::Bcast(a.data(), a.size(),
22 ParallelDescriptor::IOProcessorNumber())
23
24// See AMReX_ParallelDescriptor.H for many other Reduce functions
25ParallelDescriptor::ReduceRealSum(x);
AMReX_Vector.H): basically std::vector but provides bound-checking
when compiled with DEBUG=TRUE.1amrex::Vector<amrex::Real> myvector {1.0_rt, 2.0_rt};
AMReX_Array.H): basically std::array. 0-based index 1D array.1amrex::Array<amrex::Real, 3> arr {1.0_rt, 2.0_rt, 3.0_rt};
AMREX_GPU_HOST_DEVICE, so that it works on both HOST and DEVICE.
But modern C++ after C++11 should make std::array compatible on both HOST and DEVICE.
But it is still safer to use this. Access via arr[2].1amrex::GpuArray<amrex::Real, 3> arr {1.0_rt, 2.0_rt, 3.0_rt};
arr(1,2).1amrex::Array2D<amrex::Real, 1, 2, 1, 3, amrex::Order::C> arr
2 {{1.0_rt, 2.0_rt, 3.0_rt,
3 4.0_rt, 5.0_rt, 6.0_rt}};
AMReX_IntVect.H)IntVect: integer vector in AMREX_SPACEDIM dimensional space.
It essentially records the index representing a given cell.
1amrex::IntVect iv(AMREX_D_DECL(i, j, k));
2iv[0] == i;
3iv[1] == j;
4iv[2] == k;
It is essentially
1int vec[AMREX_SPACEDIM];
1IntVect iv(AMREX_D_DECL(19, 0, 5));
2IntVect iv2(AMREX_D_DECL(4, 8, 0));
3iv += iv2; // iv is now (23,8,5)
4iv *= 2; // iv is now (46,16,10);
AMR codes often need to deal with refinement and coarsening. The refinement operation can be done simply using the multiplication operation or , i.e.
1IntVect iv(AMREX_D_DECL(127,127,127));
2IntVect rr(AMREX_D_DECL(2,2,2,));
3
4// Refine all component by 2
5iv *= 2;
6
7// Or use refine function?
8iv.refine(rr);
For coarsening, one should ideally use builtin coarsen function.
For example, i=-1/2 would give i=0, but what we usually want is i=-1.
1IntVect iv(AMREX_D_DECL(127,127,127));
2IntVect coarsening_ratio(AMREX_D_DECL(2,2,2));
3iv.coarsen(2); // Coarsen each component by 2
4iv.coarsen(coarsening_ratio); // Component-wise coarsening
5
6// Return an IntVect w/o modifying iv
7const auto& iv2 = amrex::coarsen(iv, 2);
8
9// iv not modified
10IntVect iv3 = amrex::coarsen(iv, coarsening_ratio);
AMReX_IndexType.H)This class defines whether an index is cell-based or node-based along each dimesion. Default is to use cell-based.
One can construct IndexType with an IntVect with:
1// Node centered in x-dir and cell-centered in y and z-dir
2IndexType xface(IntVect{AMREX_D_DECL(1,0,0)});
This class provides other useful functions to tell whether its node or cell centered.
1// True if the IndexType is cell-based in all directions.
2bool cellCentered () const;
3
4// True if the IndexType is cell-based in dir-direction.
5bool cellCentered (int dir) const;
6
7// True if the IndexType is node-based in all directions.
8bool nodeCentered () const;
9
10// True if the IndexType is node-based in dir-direction.
11bool nodeCentered (int dir) const;
AMReX_Box.H)A Box defines a rectangular indexing region of AMREX_SPACEDIM dimension.
It has 2 components:
IndexType defining whether it is cell-centered or node-centeredIntVect defining the lower and upper corners of the rectangular region.A typical way of defining the Box is the following
1IntVect lo(AMREX_D_DECL(64,64,64)); // Lower corner
2IntVect hi(AMREX_D_DECL(127,127,127)); // Upper corner
3IndexType typ({AMREX_D_DECL(1,1,1)}); // Node center in all dir
4
5// Construct a cell-centered Box by default
6Box cc(lo,hi);
7
8// Construct a nodal Box, where we explicitly pass in IndexType
9Box nd(lo,hi+1,typ);
We can convert Box from cell-centered to node-centered via
1// Construct all cell-centered box by default
2Box b0 ({64,64,64}, {127,127,127});
3
4// Convert cell-centered box to node-centered.
5Box b1 = surroundingNodes(b0); // A new Box with type (node, node, node)
6Print() << b1; // ((64,64,64) (128,128,128) (1,1,1))
7Print() << b0; // Still ((64,64,64) (127,127,127) (0,0,0))
8
9// Convert the node-centered box to cell-centered again
10Box b2 = enclosedCells(b1); // A new Box with type (cell, cell, cell)
11if (b2 == b0) { // Yes, they are identical.
12 Print() << "b0 and b2 are identical!\n";
13 }
14
15// Use convert() to not have uniform cell- or node-centered box
16Box b3 = convert(b0, {0,1,0}); // A new Box with type (cell, node, cell)
17Print() << b3; // ((64,64,64) (127,128,127) (0,1,0))
18
19b3.convert({0,0,1}); // Convert b0 to type (cell, cell, node)
20Print() << b3; // ((64,64,64) (127,127,128) (0,0,1))
21
22b3.surroundingNodes(); // ((64,64,64) (128,128,128) (1,1,1))
23b3.enclosedCells(); // ((64,64,64) (127,127,127) (0,0,0))
Sometimes given a Box, we want to know its corners. There are multiple ways
to do this
1Box bx ({64,64,64}, {127,127,127});
2
3// Get the small end (lower corner) of the Box
4IntVect lo = bx.smallEnd();
5
6 // Get the big end (upper corner) along y-dir
7bx.bigEnd(1);
8
9// Get const pointer the lower end and upper end (Legacy interface)
10const int* lo = bx.loVect();
11const int* hi = bx.hiVect();
12
13// Use lbound and ubound to get lower and upper corner
14// they return Dim3 instead of IntVect -- prefered for GPU access
15const Dim3 lo = lbound(bx);
16const Dim3 hi = ubound(bx);
Lastly, we want to do refinement and coarsen of a given Box.
Note that the behavior depends on whether the Box is cell-centered or node-centered.
A refined Box always covers the same physical domain as the original Box.
A coarsened Box does NOT always cover the same physical domain – see example below.
1// Cell-centered Box
2Box ccbx ({16,16,16}, {31,31,31});
3ccbx.refine(2);
4Print() << ccbx; // ((32,32,32) (63,63,63) (0,0,0))
5Print() << ccbx.coarsen(2); // ((16,16,16) (31,31,31) (0,0,0))
6
7// Node-centered Box
8Box ndbx ({16,16,16}, {32,32,32}, {1,1,1});
9ndbx.refine(2);
10Print() << ndbx; // ((32,32,32) (64,64,64) (1,1,1))
11Print() << ndbx.coarsen(2); // ((16,16,16) (32,32,32) (1,1,1))
12
13// Node-centered in x-dir
14Box facebx ({16,16,16}, {32,31,31}, {1,0,0});
15facebx.refine(2);
16Print() << faceb x; // ((32,32,32) (64,63,63) (1,0,0))
17Print() << facebx.coarsen(2); // ((16,16,16) (32,31,31) (1,0,0))
18
19// Uncoarsenable Box
20Box uncbx ({16,16,16}, {30,30,30});
21Print() << uncbx.coarsen(2); // ((8,8,8), (15,15,15));
22Print() << uncbx.refine(2); // ((16,16,16), (31,31,31));
23 // Different from the original!
One can also expand a given box – this is almost always used to get the ghost cells.
1// Increase a box in all directions by 2 cells,
2// basically the lower corners are subtracted by 2 in all dir
3// and upper corners are added by 2 in all dir
4const Box& obx = amrex::grow(bx, 2);
5
6// Only expand in y-dir by 4 cells.
7const Box& qbx = amrex::grow(bx, 1, 4);
8
9// Only grow the lower and upper corner
10// One can also use negative values to do 'shrink'
11Box gbx = amrex::growLo(bx, 1, 3);
12Box gbx = amrex::growHi(bx, 1, 3);
These are plain structs are holds 3 fields: x,y,z.
1struct Dim3 { int x; int y; int z; };
2struct XDim3 { Real x; Real y; Real z; };
Note that they also have 3 fields regardless of dimensions.
And we see that IntVect can be easilly used to convert to Dim3.
1IntVect iv(...);
2Dim3 d3 = iv.dim3();
Given a Box, it is often used to get upper and lower bound,
and use them to write a loop
1Box bx(...);
2Dim3 lo = lbound(bx);
3Dim3 hi = ubound(bx);
4for (int k = lo.z; k <= hi.z; ++k) {
5 for (int j = lo.y; j <= hi.y; ++j) {
6 for (int i = lo.x; i <= hi.x; ++i) {
7 }
8 }
9 }
AMReX_RealBox.H)A RealBox stores the physical location in floating-point numbers of the
lower and upper corners of a rectangular domain, rather than storing
indices compared to Box.
1// A rectangular box whose length is [-1.0, 1.0] for each dir.
2RealBox rb({AMREX_D_DECL(-1.0,-1.0,-1.0)},
3 {AMREX_D_DECL( 1.0, 1.0, 1.0)});
4
5// Lower and upper corner in physical location.
6// One can also specify dir.
7auto lo = rb.lo();
8auto hi = rb.hi();
AMReX_Geometry.H)With the Box and RealBox specified, we can now create the Geometry.
1int n_cell = 64;
2
3// This defines a Box with n_cell cells in each direction.
4Box domain(IntVect{AMREX_D_DECL( 0, 0, 0)},
5 IntVect{AMREX_D_DECL(n_cell-1, n_cell-1, n_cell-1)});
6
7// This defines the physical box, [-1,1] in each direction.
8RealBox real_box({AMREX_D_DECL(-1.0,-1.0,-1.0)},
9 {AMREX_D_DECL( 1.0, 1.0, 1.0)});
10
11// This says we are using Cartesian coordinates
12int coord = 0;
13
14// This sets the boundary conditions to be doubly or triply periodic
15Array<int,AMREX_SPACEDIM> is_periodic {AMREX_D_DECL(1,1,1)};
16
17// This defines a Geometry object
18Geometry geom(domain, real_box, coord, is_periodic);
Now we normally want a Geometry object for every AMR-level.
Hence, the Box and RealBox here represents the full-domain of that level.
Every level needs a Geometry object since stuff like CellSize, volume differs.
Here are some information you can get from Geometry
1// Lower corner of the physical domain
2// Returns GpuArray<Real,AMREX_SPACEDIM>
3const auto problo = geom.ProbLoArray();
4
5// y-direction upper corner
6Real yhi = geom.ProbHi(1);
7
8// Cell size for each direction.
9const auto dx = geom.CellSizeArray();
10
11// Index domain -- the Box containing upper lower indices
12const Box& domain = geom.Domain();
13
14// Is periodic in x-direction?
15bool is_per = geom.isPeriodic(0);
AMReX_BoxArray.H)BoxArray is a collection of Boxes.
You can create a BoxArray from a single Box, i.e. a single Box
representing the full domain. And then chop up the BoxArray
into smallers Box.
1// Create a Box representing the full domain
2Box domain(IntVect{0,0,0}, IntVect{127,127,127});
3
4// Make a new BoxArray out of a single Box
5BoxArray ba(domain);
6
7// Get the number of boxes in the BoxArray, which is 1 for now.
8Print() << "BoxArray size is " << ba.size() << "\n";
9
10// Forces each Box in BoxArray to have sides <= maxSize
11// so chop up the single Box into smaller boxes 64^3 cells
12ba.maxSize(64);
13Print() << ba; // Now shows 8 Boxes
One thing to note is that BoxArray is a global data structure.
This means that all processes knows about the whole BoxArray.
It holds all Boxes in a collection, but a single process in parallel
can only own some of its Boxes via domain decomposition – this
refers to the actual data, i.e. FArrayBox or FAB, see next few sections.
Similar to Box, BoxArray has an IndexType, and all Boxes in the BoxArray
will have same type as the BoxArray itself. One can convert the
IndexType of a BoxArray via
1// Create a cell-centered BoxArray from a single Box
2BoxArray cellba(Box(IntVect{0,0,0}, IntVect{63,127,127}));
3
4// Chop up BoxArray so that the length of each Box is < 64
5cellba.maxSize(64);
6
7// Make a copy
8BoxArray faceba = cellba;
9
10// Convert BoxArray to index type (cell, cell, node)
11faceba.convert(IntVect{0,0,1});
12
13// Return an all node BoxArray
14const BoxArray& nodeba = amrex::convert(faceba, IntVect{1,1,1});
15
16// We can access individual Box in BoxArray via []
17// Note that it returns by VALUE instead of REFERENCE
18Print() << cellba[0] << "\n"; // ((0,0,0) (63,63,63) (0,0,0))
19Print() << faceba[0] << "\n"; // ((0,0,0) (63,63,64) (0,0,1))
20Print() << nodeba[0] << "\n"; // ((0,0,0) (64,64,64) (1,1,1))
BoxArray also allows all Boxes to be refined or coarsened via
1BoxArray ba(Box(IntVect{0,0,0}, IntVect{63,127,127}));
2IntVect rr {2,2,4};
3
4// Refine each Box in BoxArray in all dir by 2
5ba.refine(2);
6
7// Different refinement ratio for each dir
8ba.refine(rr);
9
10// Coarsen each Box in BoxArray by 2
11ba.coarsen(2);
12ba.coarsen(rr);
AMReX_DistributionMapping.H)DistributionMapping uses an algorithm to describe which process
owns the data living on the domains specified by the Boxes in BoxArray.
One can construct a DistributionMapping given a BoxArray.
1DistributionMapping dm {ba};
BaseFab is a class containing multi-dimensional array data for a single Box.
This is like a step up of a Box where not only does it know about
lower and upper corner indices, it knows about the data specified in that region.
The dimensionality of the array is AMREX_SPACEDIM+1. The additional
dimension is to account for the number of components, say we have
density, momentum, and energy. Then we have 5 components.
1// First create a single Box
2Box bx(IntVect{-4,8,32}, IntVect{32,64,48});
3
4// Specify the number of components
5int numcomps = 4;
6
7// Create the array
8BaseFab<Real> fab(bx,numcomps);
Now usually we don’t work BaseFab, but classes derived from it.
FArrayBox = BaseFab<Real>IArrayBox = BaseFab<int>One can easily get the associated Box and the number of components via
1// Get the Box
2Box bx = fab.box();
3
4// Get the number of components
5int ncomp = fab.nComp();
One can also change the Box and number of components a FArrayBox is defined on
via the resize method or define method. Both are accomplishes the same thing,
but use define when first constructing the object and resize when changing
an existing one.
1// using define
2fab.define(box, ncomp);
3
4// using resize
5fab.resize(box, ncomp);
One can do simple arithmetic operations using builtin functions
1// Create the Box and specify number of components
2Box box(IntVect{0,0,0}, IntVect{63,63,63});
3int ncomp = 2;
4
5// Create the ArrayBox holding the data
6FArrayBox fab1(box, ncomp);
7FArrayBox fab2(box, ncomp);
8
9// Fill fab1 with 1.0
10fab1.setVal(1.0);
11
12// Multiply component 0 by 10.0
13fab1.mult(10.0, 0);
14
15// Fill fab2 with 2.0
16fab2.setVal(2.0);
17
18// For all components, fab2 <- a * fab1 + fab2
19Real a = 3.0;
20fab2.saxpy(a, fab1);
For more complicated operations, one needs to work the array directly.
The actual 4D array object in a BaseFab class is contained using Array4.
Here is an example
1// Given FArrayBox and IArrayBox
2FArrayBox afab(...), bfab(...);
3IArrayBox ifab(...);
4
5// Get the appropriate array containing data using fab.array()
6Array4<Real> const& a = afab.array();
7Array4<Real const> const b = bfab.const_array();
8Array4<int const> m = ifab.array();
9
10// Get the lower and upper bound indices of the Box
11// and the number of componenets
12Dim3 lo = lbound(a);
13Dim3 hi = ubound(a);
14int nc = a.nComp();
15
16// Do a 4D loop loop
17for (int n = 0; n < nc; ++n) {
18 for (int k = lo.z; k <= hi.z; ++k) {
19 for (int j = lo.y; j <= hi.y; ++j) {
20 for (int i = lo.x; i <= hi.x; ++i) {
21 if (m(i,j,k) > 0) {
22 a(i,j,k,n) *= 2.0;
23 } else {
24 a(i,j,k,n) = 2.0*a(i,j,k,n) + 0.5*(b(i-1,j,k,n)+b(i+1,j,k,n));
25 }
26 }
27 }
28 }
29 }
30
31// Using ParallelFor instead of manual loop (GPU friendly)
32amrex::ParallelFor(bx, NUM_STATE,
33[=] AMREX_GPU_DEVICE (int i, int j, int k, int n) noexcept
34{
35 if (m(i,j,k) > 0) {
36 a(i,j,k,n) *= 2.0;
37 } else {
38 a(i,j,k,n) = 2.0*a(i,j,k,n) + 0.5*(b(i-1,j,k,n)+b(i+1,j,k,n));
39 }
40}
FabArray is a collection of BaseFab, similar to BoxArray is a collection of Box.
Similar to BoxArray. Similar to BaseFab, we usually don’t work with FabArray
directly, but classes derived from it. This is usually
MultiFab = FabArray<FArrayBox>> = FabArray<BaseFab<Real>>iMultiFab = FabArray<IArrayBox>> = FabArray<BaseFab<int>>Now because FabArray is a parallel data structure where the data, i.e. FArrayBox or FAB,
are distributed among parallel processes. For each process, a FabArray contains only
the FAB objects owned by that process and only operates on that data.
For operations that require data owned by other processes, remote communications are needed,
so we FabArray needs DistributionMapping that specifies which process owns which Box.
So for each process, FabArray knows about the full BoxArray, but only the FAB that is assigned
to that process.
To create a FabArray or MultiFab,
1// ba is BoxArray
2// dm is DistributionMapping
3// 4 components
4int ncomp = 4;
5
6// 1 ghost cell, so all Boxes have grown by 1 cell.
7// If BoxArray has Box{(7,7,7) (15,15,15)}
8// then the one used for FArrayBox in MultiFab
9// is then Box{(6,6,6) (16,16,16)}.
10int ngrow = 1;
11MultiFab mf(ba, dm, ncomp, ngrow);
When creating the MultiFab, we use the individual Box in BoxArray
to create the individual FArrayBox. So a MultiFab contains both
BoxArray and a collection of FArrayBox which also knows about its
corresponding Box. Now usually these two Box are the same, except
when we create MultiFab with ghost cells.
Internally the FArrayBox is created via
1// Create the FAB but with grown Box
2FArrayBox fab(amrex::grow(ba[0], ngrow),ncomp);
So the Box used to create the FArrayBox is enlarged, but the
Box stored in MultiFab stays unchanged and this is contains
the valid zones. And notice that FArrayBox itself doesn’t have
the concept of ghost cells, it is simply allocated using a grown
Box.
For a given MultiFab, we can iterate through all the FArrayBox
via MFIter
1for (MFIter mfi(mf); mfi.isValid(); ++mfi)
2{
3 // Valid Box without ghost cells
4 const Box& valid = mfi.validbox();
5
6 // With ghost cells
7 const Box& fabbx = mfi.fabbox();
8
9 // Get the Array
10 Array4<Real const> const& a = mf.array(mfi);
11
12 Print() << valid << "\n";
13 Print() << fabbx << "\n";
14}
We can also create a MultiFab from another MultiFab.
1// mf0 is an already defined MultiFab
2// Get properties like the BoxArray and DistributionMap
3const BoxArray& ba = mf0.boxArray();
4const DistributionMapping& dm = mf0.DistributionMap();
5int ncomp = mf0.nComp();
6int ngrow = mf0.nGrow();
7
8// Create new MultiFabs with existing properties
9MultiFab mf1(ba,dm,ncomp,ngrow);
10
11// change the IndexType to face centered at x-dir for x-flux
12MultiFab xflux(amrex::convert(ba, IntVect{1,0,0}), dm, ncomp, 0);
We can also do some basic operations on a MultiFab or between
MultiFab built with the same BoxArray and DistributionMapping.
1// Minimum value in component 3 of MultiFab mf
2// no ghost cells included
3Real dmin = mf.min(3);
4
5// Maximum value in component 3 of MultiFab mf
6// including 1 ghost cell
7Real dmax = mf.max(3,1);
8
9// Set all values to zero including ghost cells
10mf.setVal(0.0);
11
12// Add mfsrc to mfdst
13MultiFab::Add(mfdst, mfsrc, sc, dc, nc, ng);
14
15// Copy from mfsrc to mfdst
16MultiFab::Copy(mfdst, mfsrc, sc, dc, nc, ng);
17
18// MultiFab mfdst: destination
19// MultiFab mfsrc: source
20// int sc : starting component index in mfsrc for this operation
21// int dc : starting component index in mfdst for this operation
22// int nc : number of components for this operation
23// int ng : number of ghost cells involved in this operation
24// mfdst and mfsrc may have more ghost cells
Now the BoxArray used to define the MultiFab is usuallly non-intersecting,
except when due to a nodal index type. But if MultiFab has ghost cells,
then FArrayBox which was constructed via a grown BoxArray can have overlaps,
even though the BoxArray is still the original one, which means that
the Box from FArrayBox is larger than the one in BoxArray.
This means that parallel communication is needed to fill ghost cells
from other FArrayBox which can be possibly owned by a different process.
To do this, use FillBoundary
1MultiFab mf(...parameters omitted...);
2Geometry geom(...parameters omitted...);
3
4// Fill ghost cells for all components
5// Periodic boundaries are not filled.
6mf.FillBoundary();
7
8// Fill ghost cells for all components
9// Periodic boundaries are filled.
10// This matters when grids are touching physical periodic boundaries.
11// Note that if physical boundary is not periodic, like outflow,
12// it just leaves them as garbage data since no grid overlap with ghost cells.
13mf.FillBoundary(geom.periodicity());
14
15// Fill 3 components starting from component 2
16mf.FillBoundary(2, 3);
17mf.FillBoundary(geom.periodicity(), 2, 3);
Another type of parallel communication is copying data from One
MultiFab to another MultiFab with a different or same BoxArray
with a different DistributionMapping.
The copy will only be performed on the regions of intersection
as one would expect.
1mfdst.ParallelCopy(mfsrc, compsrc, compdst, ncomp, ngsrc, ngdst, period, op);
We previously discussed manipulating individual FArrayBox in MultiFab
using MFIter. Let’s discuss this in more detail. Here is an example
1for (MFIter mfi(mf); mfi.isValid(); ++mfi) // Loop over grids
2{
3 // This is the valid Box of the current FArrayBox.
4 // By "valid", we mean the original ungrown Box in BoxArray.
5 const Box& box = mfi.validbox();
6
7 // A reference to the current FArrayBox in this loop iteration.
8 FArrayBox& fab = mf[mfi];
9
10 // Obtain Array4 from FArrayBox. We can also do
11 // Array4<Real> const& a = mf.array(mfi);
12 Array4<Real> const& a = fab.array();
13
14 // Call function f1 to work on the region specified by box.
15 // Note that the whole region of the Fab includes ghost
16 // cells (if there are any), and is thus larger than or
17 // equal to "box".
18 // So we don't have to worry about accessing data outside
19 // of box when doing loops.
20 f1(box, a);
21}
22
23// f1 operates on Box and can be something like
24void f1 (Box const& bx, Array4<Real> const& a) {
25 // Or use ParallelFor...
26 const auto lo = lbound(bx);
27 const auto hi = ubound(bx);
28 const int ncomp = a.nComp();
29 for (int n = 0; n < ncomp; ++n) {
30 for (int k = lo.z; k <= hi.z; ++k) {
31 for (int j = lo.y; j <= hi.y; ++j) {
32 for (int i = lo.x; i <= hi.x; ++i) {
33 a(i,j,k,n) = ...;
34 }
35 }
36 }
37 }
38}
Here note that MFIter only loops over grid owned by this process.
We can also work with multiple MultiFab for a given loop. Note that
they need to have the same DistributionMapping, but not necessary
the same BoxArray, i.e. they can differ due to index types.
1// U and F are MultiFabs
2for (MFIter mfi(F); mfi.isValid(); ++mfi) // Loop over grids
3{
4 const Box& box = mfi.validbox();
5
6 Array4<Real const> const& u = U.const_array(mfi);
7 Array4<Real > const& f = F.array(mfi);
8
9 f2(box, u, f);
10}
Let’s first dicuss what is Loop Tiling
Tiling, also known as cache blocking – a loop transformation technique for improving data locality. It transforms loops into tiles or smaller blocks first and then do the looping within tiles. The idea is that it minimizes reading data in memory. Here is an online webpage explaining loop tiling.
Before talking about loop tiling, let’s first discuss how CPU gets its data. Upon a memory request for CPU, it fetches the data with following flow:
RAM (GBs) → L3 (8-32MB) → L2 (256KB-1MB) → L1 (32-64KB) → registers → CPU
Also note that RAM is basically a 2D block of bytes where it has N number of rows and M number of columns. Here M is the cache line size, which is typically 64 bytes, i.e. each row has 64 memory addresses. For example, a single precision number uses 4 bytes and double precision uses 8 bytes. Now when a CPU fetches the desired data, it gets all data of that specific row contains the data requested, i.e. it gets 64 bytes of data every single time. Now this is why it is faster to access data adjacent in memory address, since the neighboring data are already loaded in cache.
The order at which CPU checks whether data exists is the reverse order of the above. Each level is smaller, faster and closer to the CPU. Now once the data gets filled out in those levels, it gets removed or data eviction following some policies. The most common one is Least Recently Used or LRU, which removes data that is not used in the longest time.
Now let’s consider this example
1for (int i = 0; i < n; i++) {
2 for (int j = 0; j < m; j++) {
3 c[i][j] = a[i] * b[j];
4 }
5 }
When doing this loop, CPU first gets a[0] and then subsequently loads in b[0], b[1], b[2], …, b[m]. Now if array b is large, then the earlier elements of array b gets removed since it exceeds the cache size. But if array b is small enough so that all elements can live within cache without it being evicted, then we maximize speed.
Now consider array a is small and array b is large so that we just need to do tiling on b. This means we can do the following:
1for (int jj = 0; jj < m; jj += TILE_SIZE) {
2 for (int i = 0; i < n; i++) {
3 for (int j = jj; j < MIN(jj + TILE_SIZE, m); j++) {
4 c[i][j] = a[i] * b[j];
5 }
6 }
7 }
Now the access pattern goes like:
Now notice that we’re accessing the same subset of array b, from b[0] to b[TILZE_SIZE-1] over and over again. Assume this is small enough, and so these elements never got evicted from cache. Hence we eliminated the stall of waiting to fetch these data from RAM. Now if array a is also large, then we should do the same tiling on a, which then looks like
1for (int ii = 0; ii < n; ii += TILE_SIZE_I) {
2 for (int jj = 0; jj < m; jj += TILE_SIZE_J) {
3 for (int i = ii; i < MIN(n, ii + TILE_SIZE_I); i++) {
4 for (int j = jj; j < MIN(m, jj + TILE_SIZE_J); j++) {
5 c[i][j] = a[i] * b[j];
6 }
7 }
8 }
9}
In this case, we would want the total tile size, TILE_SIZE_I * TILE_SIZE_J to be small enough so that all these subset of array a and b can live in cache. Now ideally we want the total size of the data for each tile to be smaller than the L1 cache. However, we don’t want the TILE_SIZE to be too small for two reasons:
When should we use it?
1for (int i = 0; i < n; i++) {
2 for (int j = 0; j < n; j++) {
3 // access to array a is column-major
4 // use tiling to reduce time between
5 // consecutive access to neighboring of array a
6 b[i][j] = a[j][i];
7 }
8}
From previous example we see that the major downside is that the syntax becomes ugly. But we can easily enable tiling with the following synteax when looping over boxes.
1// * true * turns on tiling
2for (MFIter mfi(mf,true); mfi.isValid(); ++mfi) // Loop over tiles
3{
4 // tilebox() instead of validbox()
5 const Box& box = mfi.tilebox();
6
7 ...
8}
We can also specify the tile size when defining MFIter,
1// No tiling in x-direction. Tile size is 16 for y and 32 for z.
2// The default tile size is IntVect(1024000,8,8)
3for (MFIter mfi(mf,IntVect(1024000,16,32)); mfi.isValid(); ++mfi) {...}
Note if the tile size is larger than the grid size, then it means that tiling is disabled in that direction. So the example above has no tiling in x-dir.
ParallelFor allows you to not manually write out all the loops when
iterating over a Box. This works on either GPU or CPU so its convenient to use.
Here is an example
1#ifdef AMREX_USE_OMP
2#pragma omp parallel if (Gpu::notInLaunchRegion())
3#endif
4 for (MFIter mfi(mfa,TilingIfNotGPU()); mfi.isValid(); ++mfi)
5 {
6 const Box& bx = mfi.tilebox();
7 Array4<Real> const& a = mfa[mfi].array();
8 Array4<Real const> const& b = mfb[mfi].const_array();
9 Array4<Real const> const& c = mfc[mfi].const_array();
10
11 // Assume a single component.
12 ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k)
13 {
14 a(i,j,k) += b(i,j,k) * c(i,j,k);
15 });
16 }
Here are some other versions
1// 1D for loop
2ParallelFor(N, [=] AMREX_GPU_DEVICE (int i) { ... });
3
4// 4D for loop
5ParallelFor(box, numcomps,
6 [=] AMREX_GPU_DEVICE (int i, int j, int k, int n) { ... });