# MultiFab Tutorial

**Time to Complete**: 15 mins

**GOALS:**Create a MultiFab

Write Data to a MultiFab

Plot MultiFab Data

In this tutorial we focus on using the MultiFab data structure. We will learn the steps to create a MultiFab, load it with data, and plot it. At the end, we make some comments about writing code that works for both 2 and 3 dimensions.

A video companion to this tutorial is available on YouTube. You can access it by clicking on the thumbnail to the right.

## What is a MultiFab?

A MultiFab is a C++ class in AMReX the stores and operates on multidimensional arrays in parallel. It contains:

Grid information in the form of a BoxArray that contains one or more components (scalar values) for a single level of the mesh.

A distribution map, that allows for parallel processing of data in the MultiFab.

Ghost cells that facilitate a variety of mesh refinement, boundary conditions, and particle algorithms.

## Defining a MultiFab

To use the MultiFab type we first need to include the header file,

```
#include <AMReX_MultiFab.H>
```

To define a MultiFab,

```
amrex::MultiFab mf(ba, dm, ncomp, ngrow);
```

Defining a MultiFab in this way requires 4 inputs:

`ba`

, a BoxArray.

`dm`

, a DistributionMapping.

`ncomp`

, the number of components (scalar values) to store in the MultiFab.

`ngrow`

, the number of layers of ghost cells.

The value of `ncomp`

is straight-forward and can be specified when
the MultiFab is defined. The correct number of ghost cells, `ngrow`

, is important
for many of the algorithms in AMReX such as boundary behavior, mesh refinement, etc.
For this reason it should be carefully selected based on application requirements.
The DistributionMapping, `dm`

, directs parts of the domain to different processes
for parallel computation. In this example, the default suffices and can be easily
defined once the BoxArray is configured. Defining the BoxArray requires a few steps
that are presented in the next section.

### Number of Components and Ghost Cells

The value of `ncomp`

determines how many scalar values to store in the MultiFab.
For example, if we want to store the value of phi at every point on the grid,
thats one component for the MultiFab. Then we would set:

```
int ncomp = 1;
```

The number of layers of ghost cells around the boundary of the domain is determined
by the value of `ngrow`

and will depend on your application. In this tutorial
we will not be doing any operations that require ghost cells,
so we will set the value to 0.

```
int ngrow = 0;
```

### BoxArray: Abstract Domain Setup

The BoxArray contains a list of boxes that cover the domain. To define a BoxArray, we will therefore first need to define the domain and some of its properties. In this example, we chose a 3-dimensional domain and therefore have three inputs. The steps are:

Define the lower and upper indices of the domain:

amrex::IntVect dom_lo(0,0,0); amrex::IntVect dom_hi(n_cell-1, n_cell-1, n_cell-1);We use two IntVects to define the high and low indices of the domain. In the case of the high indices, we define it as

`n_cell-1`

.`n_cell`

represents the number of cells we want in each dimension. Typically its value is read from the inputs file.Define a single Box with the specified domain:

amrex::Box domain(dom_lo, dom_hi);Define the BoxArray object using the domain:

amrex::BoxArray ba(domain);Next we define a maximum size for grids (“chunks”) in the domain:

ba.maxSize(max_grid_size);

`max_grid_size`

will be read from the inputs file at runtime. Its value will effect how AMReX processes the data on the MultiFab in parallel.

### Distribution Mapping

Once the BoxArray is defined. We can define a DisbributionMapping. The DistributionMapping will determine how parts of the domain are divided among processes. For us, the default behavior is sufficient. Therefore we write,

```
amrex::DistributionMapping dm(ba);
```

At this point, we have defined the four necessary parts of the MultiFab and can create it with the line,

```
amrex::MultiFab mf(ba, dm, ncomp, ngrow);
```

## Adding Data to the MultiFab

Now that we have setup an abstract domain, we want to provide physical information about the boxes within it so that we can fill it with useful data. To do that, we will define a Geometry object, and use a nested loop structure that efficiently iterates over all the boxes in the domain. AMReX will automatically decide how to parallelize these operations with MPI ranks or GPU accelerators.

### Geometry: Add Physical Properties

For the geometry, we will define the size of the box, the coordinate system and the boundary conditions. It is also convenient to derive the dimensions of each cell at this time.

Define the physical dimensions of the box:

amrex::RealBox real_box ({ 0., 0., 0.} , {1., 1., 1.});In this case, we use

`[0,1]`

for all directions.Define the Geometry object with the properties specified above:

amrex::Geometry geom(domain, &real_box);*In preparation for future operations on the MultiFab, extract the physical dimensions of each cell:

amrex::GpuArray<amrex::Real,3> dx = geom.CellSizeArray();This commands creates a 1-dimensional array of

`amrex::Real`

–single or double floating point– values that correspond to the physical dimensions of each cell side, i.e.`dx[0]`

contains the length of the cell in the x-direction,`dx[1]`

, contains the length in the y-direction, and so on.

### Loop Structure: MFIter and ParallelFor

To traverse the elements of the MultiFab we will use two loops. The first loop is a MultiFab iterator,
called `MFIter`

, and the second loop is a `ParallelFor`

. The `MFIter`

loop efficiently divides
and distributes work on MultiFab data among processors. It does this by telling each processor,
“iterate through all the boxes on the MultiFab, but only do work on the boxes assigned
to you.”

Inside the `MFIter`

for loop, the
`ParallelFor`

function behaves like a triple-nested loop over the `i,j,k`

coordinates in the box.
Beginning and ending values for each index are derived
from the inputs, and do not need to be explicitly stated. When a GPU backend is enabled, `ParallelFor`

will launch a thread on the GPU to process each `i,j,k`

iteration in parallel.
On CPU only, `ParallelFor`

will enable tiling to process the data in the most efficient way for the hardware available.

Below is an example of typical usage of `MFIter`

and `ParallelFor`

to fill the MultiFab with data:

```
for(amrex::MFIter mfi(mf); mfi.isValid(); ++mfi){
const amrex::Box& bx = mfi.validbox();
const amrex::Array4<amrex::Real>& mf_array = mf.array(mfi);
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k){
amrex::Real x = (i+0.5) * dx[0];
amrex::Real y = (j+0.5) * dx[1];
amrex::Real z = (k+0.5) * dx[2];
amrex::Real rsquared = ((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5) + (z-0.5)*(z-0.5))/0.01;
mf_array(i,j,k) = 1.0 + std::exp(-rsquared);
});
}
```

In the first line,

```
for(amrex::MFIter mfi(mf); mfi.isValid(); ++mfi){
```

the MFIter object `mfi`

is defined from the MultiFab `mf`

. `mfi.isValid()`

, tells the for loop
to work on only the parts of the domain designated to that processor by `mf`

including growth or ghost cells.

The following two lines,

```
const amrex::Box& bx = mfi.validbox();
const amrex::Array4<amrex::Real>& mf_array = mf.array(mfi);
```

define a new box, `bx`

that represents the current section of the grid being iterated on. To access
and store data in the components
of the MultiFab within this section of the grid, we cast `mf`

as an Array4 object called `mf_array`

.
We declare `mf_array`

by reference so that
values changed within the `ParallelFor`

loop, change values in the components of the MultiFab `mf`

.
Within `ParallelFor`

each component of `mf`

can be accessed by
calling `mf_array(i,j,k,n)`

where `i,j,k`

represents the location and `n`

represents the
\((n+1)^{\text th}\) component in the MultiFab. Moreover, the first or \(0^{\text th}\) component, can be accessed
by dropping the 4th index, i.e. `mf_array(i,j,k)`

is equivalent to `mf_array(i,j,k,0)`

.

The line with the `ParallelFor`

function,

```
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k){
```

calls a custom looping function that traverses the three dimensions of the box, `bx`

. The range of values
for `i,j,k`

are determined by the `ParallelFor`

from the box information passed to it.
When a GPU backend is enabled, code inside `ParallelFor`

will be executed in parallel on the GPU. Functions called
within this loop are often called “lambdas”. Note that, because `ParallelFor`

is a function call, and
not a simple for loop, the closing brackets are `});`

.

The remaining lines inside `ParallelFor`

are,

```
amrex::Real x = (i+0.5) * dx[0];
amrex::Real y = (j+0.5) * dx[1];
amrex::Real z = (k+0.5) * dx[2];
amrex::Real rsquared = ((x-0.5)*(x-0.5) + (y-0.5)*(y-0.5) + (z-0.5)*(z-0.5))/0.01;
mf_array(i,j,k) = 1.0 + std::exp(-rsquared);
```

The first three lines translate integer indices to their `(x,y,z)`

location in the domain. In this case,
the `+0.5`

indicates that we want the :math:0^{text th} index for `i`

, to represent the middle of
the first cell in the x-direction. Together, treating all three variables in this way
indicates cell-center data.

The last line, stores the calculated value in the appropriate location of `mf_array`

which, in turn, stores it in the corresponding location of the MultiFab `mf`

because
of the way we declared `mf_array`

.

So far, we’ve defined and filled a MultiFab with a single component with data. The filling operation will be done in parallel when compiled with a GPU backend. When compiled for CPU-only computation, AMReX will use tiling to increase performance (see AMReX User’s Guide – Tiling). Configuring these performance optimizations would normally require different lines of code, however, AMReX handles these changes automatically for portable performance.

## Plotting MultiFab Data

AMReX can plot MultiFab data with a single function call. To access the plotting functions we include the header,

```
#include <AMReX_PlotFileUtil.H>
```

at the top of the file. We can then write the line,

```
WriteSingleLevelPlotfile("plt001", mf, {"comp0"}, geom, 0., 0);
```

The first input takes the plotfile name, “plt001”.
The second, is the MultiFab that contains th e data we want to plot, `mf`

.
The parameter, `{"comp0"}`

, labels the
first component as “comp0”. For multiple components, its necessary to pass additional
variable names, such as `{"comp0","comp1"}`

. The third parameter, `geom`

, is the
Geometry we defined for the MultiFab above. The
last parameter specifies the level of the MultiFab. In this tutorial, we only had a
single MultiFab to represent the zeroth level of the mesh, therefore we pass `0`

for the level.

### Visualizing the Plotfile

The call above to `WriteSingleLevelPlotfile`

will produce a plotfile in the form
of a directory that contains a Header file, and several subdirectories. The data can be
visualized by passing this directory to one of several visualization software packages.
Information on how to do this is available in the AMReX User’s Guide Visualization
section.

## Conclusion

In this tutorial we described the steps to define the domain, and physical properties
of our simulation space. We then demonstrated how to initialize and store scalar values
on a grid across
this domain in the form of a 3-dimensional MultiFab data structure with a single component.
Initializing the data involved the `MFIter`

, `ParallelFor`

nested-“loop” structure that
required casting the MultiFab as an Array4 for `i,j,k`

access to the active section of the
grid. Finally, we showed the commands to write out the MultiFab data to a plotfile that
can be visualized with several software packages.

The complete code for this tutorial is available here.

## Additional Comments

AMReX allows for the selection of 2- or 3- dimensional simulation at compile time. For simplicity, the example above is presented as 3-dimensional only. The example below shows how to write code to assign data to a MultiFab in a way that can be automatically adapted for 2 or 3 dimensions:

```
for (MFIter mfi(phi_old); mfi.isValid(); ++mfi)
{
const Box& bx = mfi.validbox();
const Array4<Real>& phiOld = phi_old.array(mfi);
// set phi = 1 + e^(-(r-0.5)^2)
amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
{
Real x = (i+0.5) * dx[0];
Real y = (j+0.5) * dx[1];
#if (AMREX_SPACEDIM == 2)
Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5))/0.01;
#elif (AMREX_SPACEDIM == 3)
Real z= (k+0.5) * dx[2];
Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/0.01;
#endif
phiOld(i,j,k) = 1. + std::exp(-rsquared);
});
}
```

The interesting features to point out are:

the preprocessor directives,

`#if`

,`#elif`

,`AMREX_SPACEDIM`

, etc. are evaluated at compile time, and will only process the code within the appropriate`AMREX_SPACEDIM`

section. Moreover, the value of`AMREX_SPACEDIM`

can be set while compiling the AMReX application.The indices of the

`ParallelFor`

function are`i,j,k`

despite being applicable to both 2- or 3- dimensional code. AMReX has the built-in ability to revert to the two dimensional`i,j`

indices, when compiled as a 2D code.Additional code differences, not shown here, are required in the steps to define the BoxArray and Geometry of the MultiFab. An example of the needed modifications can be found in

`HeatEuation_EX0_C`

.

### A ParallelFor Only Approach

The method presented above is the most common at the time of writing the guide. However,
other methods exist to achieve similar results. For example, AMReX has added to the
capabilities of `ParallelFor`

to include the functionality of both the loops in our
previous example. To use this type of `ParallelFor`

we need to include the header,

```
#include <AMReX_MFParallelFor.H>
```

The newer approach reduces the necessary syntax. However, we still need to cast
the MultiFab as a different type, so that we can access it with `i,j,k`

indices.
We also need to explicitly pass the number of grow or ghost cells in each direction.
These two things are accomplished in the first two lines.

```
const& amrex::MultiArray4 mf_arrs = mf.arrays();
const amrex::IntVect ngs(ngrow);
amrex::ParallelFor(
mf, ngs, [=] AMREX_GPU_DEVICE(int nbx, int i, int j, int k) noexcept {
amrex::Real x = (i+0.5) * dx[0];
amrex::Real y = (j+0.5) * dx[1];
amrex::Real z = (k+0.5) * dx[2];
amrex::Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/0.01;
mf_arrs[nbx](i,j,k) = 1. + std::exp(-rsquared);
}
});
```

The other change worth mentioning is `int nbx`

, the first iterative variable in the
inputs to the `ParallelFor`

. In comparison to the first example, iterating through
this variable mimics the functionality of the `MFIter`

for loop in the first
approach shown in this tutorial.

Finally one last comment, the line,

```
const& amrex::MultiArray4 mf_arrs = mf.arrays();
```

is often written to take advantage of the compiler’s ability to determine the correct type. To do this, we replace the above line with

```
auto const& mf_arrs = mf.arrays();
```

## What’s Next?

This tutorial provided an introduction to the MultiFab data structure. The Tutorial: Heat Equation - Simple
and Example: HeatEquation_EX1_C tutorials both have source code that demonstrate MultiFab usage in a loop that evolves
the data over time. At the cost of additional complexity, the Example: HeatEquation_EX1_C tutorial makes use
of the `using namespace`

and preprocessor directives to make writing coding easier and add
additional functionality. On the other hand, the source code in Tutorial: Heat Equation - Simple more
closely resembles the code used in this example.