7 Fields

Fields are distributed mesh data structures. They combine value type, dimension, mesh, layout, and execution space into a single object that supports assignment, expressions, reductions, halo exchange, and boundary conditions.

Fields constitute the Eulerian half of Particle-in-Cell methods, since they store values on a grid.

7.1 Main classes

Class	Role
`BareField`	Core data container, with a Kokkos view to store the field data, and layout information.
`Field`	User-facing field object with mesh and boundary condition information.
`FieldLayout`	Decomposition of a global domain over MPI ranks.
`SubFieldLayout`	Layout support for subdomains.
`BConds` and `BcTypes`	Boundary condition containers and boundary condition types.
`HaloCells`	Ghost cell and halo exchange support.

7.2 Field anatomy

classDiagram
  class BareField {
    +getView()
    +getHostMirror()
    +getLayout()
    +getNghost()
    +fillHalo()
    +sum()
  }
  class Field {
    +initialize(mesh, layout, nghost)
    +updateLayout(layout, nghost)
    +setFieldBC(bc)
    +get_mesh()
    +getVolumeIntegral()
    +getVolumeAverage()
    +getFieldBC()
  }
  class BConds {
    +findBCNeighbors(field)
    +apply(field)
    +assignGhostToPhysical(field)
  }
  class FieldLayout {
    +getDomain()
    +getLocalNDIndex()
    +getNeighbors()
    +getNeighborsSendRange()
    +getNeighborsRecvRange()
  }
  class UniformCartesian {
    +getMeshSpacing()
    +getCellVolume()
    +getMeshVolume()
    +getVertexPosition()
  }
  class Mesh {
    +getOrigin()
    +getGridsize()
  }


  Field --|> BareField
  Field --> BConds
  Field --> FieldLayout
  Field --> UniformCartesian
  UniformCartesian --> Mesh

7.3 Typical field workflow

Define the index domain.
Define the mesh.
Define a field layout over the global domain.
Construct one or more fields with the mesh and layout.
Apply expressions, boundary conditions, reductions, and solver calls.

7.4 Minimal field setup

The following sets up two fields: rho, a scalar field for the charge density, and efield, a vector field for the electric field.

constexpr unsigned Dim = 3;

using Mesh_t      = ippl::UniformCartesian<double, Dim>;
using Layout_t    = ippl::FieldLayout<Dim>;
using Centering_t = Mesh_t::DefaultCentering;
using ScalarField = ippl::Field<double, Dim, Mesh_t, Centering_t>;
using Vector_t    = ippl::Vector<double, Dim>;
using VectorField = ippl::Field<Vector_t, Dim, Mesh_t, Centering_t>;

// Initialize a mesh of a unit box [0,1]^3, with 4 points.
int points = 4;
double unit_box = 1.0;

ippl::Index I(pt);
ippl::NDIndex<Dim> domain(I, I, I);

double hx = unit_box / points;
Vector_t spacing = {hx, hx, hx};
Vector_t origin  = {0.0, 0.0, 0.0};
Mesh_t mesh(domain, spacing, origin);

// specify that we parallelize in all dimensions
std::array<bool, 3> isParallel;
isParallel.fill(true);
Layout_t layout(MPI_COMM_WORLD, domain, isParallel);

// Initialize fields 
ScalarField rho;
VectorField efield;

rho.initialize(mesh, layout);
efield.initialize(mesh, layout);

rho = 0.0;
efield = Vector_t(0.0);

The Field template inherits from BareField and adds a mesh pointer plus field boundary conditions. The public type aliases in Field are useful when writing generic code:

Alias	Meaning
`Mesh_t`	Mesh type attached to the field.
`Centering_t`	Centering type. The current `Field` alias resolves to cell centering.
`Layout_t`	`FieldLayout<Dim>`.
`BareField_t`	Underlying data container type.
`view_type`	Kokkos view type used for field storage.
`BConds_t`	Boundary-condition container type for the field.

The constructor form ScalarField field(mesh, layout); is also used in tests. The default ghost depth is one cell unless a different nghost value is supplied to the constructor or initialize.

7.5 Accessing data in kernels

Field data can be accessed within a Kokkos kernel. For the most common parallel pattern, Kokkos::parallel_for, there exists an ippl::parallel_for to loop over fields in a dimension independent manner. As an example, the following loops over a field and assigns values:

auto view        = field.getView();
const int nghost = field.getNghost();

using index_array_type = typename ippl::RangePolicy<Dim>::index_array_type;
ippl::parallel_for(
    "Assign field values", field.getFieldRangePolicy(),
    KOKKOS_LAMBDA(const index_array_type& args) {
        // ippl::apply accesses the view at the given indices and obtains a
        // reference; see src/Expression/IpplOperations.h
        ippl::apply(view, args) = ...;
    });

When doing operations which require the physical position of a point in the global domain (x, y, z), we need to map from the local index of the Field data in the Kokkos View (local to each MPI rank), to the global position in the full domain.

The standard way to do this is using the domain information in the layout and the mesh spacing. For example, when assigning field values according to a function, such as a sinusoidal:

auto view        = field.getView();
const int nghost = field.getNghost();

const auto& ldom  = layout.getLocalNDIndex();
const auto hr     = mesh.getMeshSpacing();
const auto origin = mesh.getOrigin();

using index_array_type = typename ippl::RangePolicy<Dim>::index_array_type;
ippl::parallel_for(
    "Assign field values", field.getFieldRangePolicy(),
    KOKKOS_LAMBDA(const index_array_type& args) {
        ippl::Vector<int, Dim> global_indices = args + ldom.first() - nghost;
        ippl::Vector<double, Dim> global_position = (global_indices + 0.5) * hr + origin;
        ippl::apply(view, args) = sin(global_position);
    });

The layout (FieldLayout) contains information about the offset of the sub-domain owned by this MPI rank with respect to the global origin, given by ldom.first(). The local indexing also contains the ghost cells, hence why it needs to be subtracted to obtain a global index. The + 0.5 in the global position computation is needed due to cell centering.

7.6 Layout and neighbors

FieldLayout<Dim> stores:

the global NDIndex<Dim> domain,
one local NDIndex<Dim> per rank,
a boolean decomposition flag for each dimension,
neighbor ranks encoded by boundary component,
send and receive ranges for halo exchange.

For periodic domains, construct the layout with isAllPeriodic = true:

ippl::FieldLayout<Dim> layout(MPI_COMM_WORLD, domain, isParallel, true);

The layout is therefore more than a partition. It is the communication map used by fields and by particle layouts that are coupled to the same mesh.

7.7 Boundary conditions

Current field boundary conditions are stored in BConds<Field, Dim>, an array-like container with 2 * Dim faces. Face numbering follows pairs by dimension:

Face	Meaning
`0`	lower face in dimension 0
`1`	upper face in dimension 0
`2`	lower face in dimension 1
`3`	upper face in dimension 1
`2 * d`	lower face in dimension `d`
`2 * d + 1`	upper face in dimension `d`

The dimensions go from 0 up to Dim-1.

Supported field boundary-condition classes are:

Class	Behavior
`PeriodicFace<Field>`	Periodic face exchange.
`NoBcFace<Field>`	No boundary modification.
`ZeroFace<Field>`	Constant zero boundary.
`ConstantFace<Field>`	Constant value boundary.
`ExtrapolateFace<Field>`	Linear extrapolation boundary.

Only cell-centered field boundary conditions are currently implemented.

Example:

using BConds_t = ippl::BConds<ScalarField, Dim>;

BConds_t bc;
for (size_t face = 0; face < 2 * Dim; ++face) {
    bc[face] = std::make_shared<ippl::PeriodicFace<ScalarField>>(face);
}

field.setFieldBC(bc);

Mixed boundary conditions are also possible. test/field/TestFieldBC.cpp uses periodic faces in the x direction, no boundary condition on the lower y face, a constant value on the upper y face, zero on the lower z face, and extrapolation on the upper z face:

BConds_t bc;

bc[0] = std::make_shared<ippl::PeriodicFace<ScalarField>>(0);
bc[1] = std::make_shared<ippl::PeriodicFace<ScalarField>>(1);
bc[2] = std::make_shared<ippl::NoBcFace<ScalarField>>(2);
bc[3] = std::make_shared<ippl::ConstantFace<ScalarField>>(3, 7.0);
bc[4] = std::make_shared<ippl::ZeroFace<ScalarField>>(4);
bc[5] = std::make_shared<ippl::ExtrapolateFace<ScalarField>>(5, 0.0, 1.0);

field.setFieldBC(bc);

setFieldBC stores the container and calls findBCNeighbors on the field. Differential operators such as grad, div, laplace, curl, and hess call fillHalo() and apply the field boundary conditions before building the expression object.

7.8 Halo cells

Ghost cells provide local data for stencil operations at subdomain boundaries. A field has a ghost depth nghost, available through getNghost(). Halo exchange is invoked explicitly with:

field.fillHalo();

Field operations that require neighbor values call fillHalo internally. For custom kernels, users must decide whether the operation needs only owned cells of the MPI rank or also ghost cells. If the kernel needs neighbors across rank boundaries, e.g. for stencil operations, call fillHalo() and apply boundary conditions first.

7.9 Expressions and reductions

Fields support scalar assignment, field expressions, math functions, and reductions. The unit tests exercise:

Operation	Example
scalar assignment	`rho = 1.0;`
deep copy	`auto rho_copy = rho.deepCopy();`
arithmetic	`rho = rho - exact;`
sum	`double q = rho.sum();`
norm	`double err = norm(rho) / norm(exact);`
volume integral	`rho.getVolumeIntegral();`
volume average	`rho.getVolumeAverage();`
gradient	`efield = grad(phi);`
divergence	`rho = div(efield);`
curl	`result = curl(vectorField);`
Hessian	`matrixField = hess(rho);`
Laplacian	`lap = laplace(rho);`

The differential operators use the spacing obtained from the UniformCartesian mesh. For example, laplace has a pre-factor of 1 / h[d]^2 in each dimension d.

7.10 Laplacian example

test/field/TestLaplace.cpp is the compact example for differential field operators. It builds a 3-dimensinal periodic scalar field on [-1, 1]^3,

\[ u(x,y,z) = \sin(\pi x)\sin(\pi y)\sin(\pi z), \]

then compares laplace(field) with the exact form

\[ \nabla^2 u = -3\pi^2 \sin(\pi x)\sin(\pi y)\sin(\pi z). \]

The relevant operator call is deliberately simple:

Lap = laplace(field);

// Error computation
Lap = Lap - Lap_exact;
Lap = pow(Lap, 2);
double error = sqrt(Lap.sum()) / sqrt(Lap_exact.sum());

To run the test, we pass the number of grid-points (indicating the mesh spacing), and the number of times we want to repeat the operation:

srun ./TestLaplace <points-per-dimension> <iterations>

7.11 Sub-layouts

TO-DO

7.12 Documentation tasks

Remaining migration work from the legacy LaTeX manual: older where examples, a compact table of mathematical functions supported by expressions, and a side-by-side update of legacy constructor forms to current constructor forms.