The DataBox

The fundamental data type in spiner is the DataBox. A DataBox packages a multi-dimensional (up to six dimensions) array with routines for interpolating on the array and for saving the data to and loading the data from file.

To use databox, simply include the relevant header:

#include <databox.hpp>

DatBox is templated on underyling data type, which defaults to the Real type provided by ports-of-call. (This is usually a double.)

Note

The default type can be set to type float if the preprocessor macro SINGLE_PRECISION_ENABLED is defined.

Any arithmetic type is supported, although the code has only been tested carefully with floating point numbers. To set DataBox to a single type, you may wish to declare a type alias such as:

using DataBox = Spiner::DataBox<double>

Spiner is also templated on how the interpolation gridding works. This template parameter is called Grid_t. The available options at this time are:

Spiner::RegularGrid1D<T>
Spiner::PiecewiseGrid1D<T>

where here T is the arithmetic type as discussed above. The default type is RegularGrid1D. You can further alias DataBox as, for example:

using DataBox = Spiner::DataBox<double, Spiner::RegularGrid1D<double>>;

More detail on the interpolation gridding is available below and in the interpolation section.

The third template parameter, GradStencil, selects the implementation of interpToRealWithGrads and defaults to Spiner::InterpolationStencils::Linear:

using DataBox = Spiner::DataBox<double, Spiner::RegularGrid1D<double>,
                               MyGradStencil>;

The built-in Spiner::InterpolationStencils::Rational policy enables rational interpolation with gradients in 1D and 2D.

The policy supplies portable static overloads interpToRealWithGrads(x1, grids, data) and interpToRealWithGrads(x2, x1, grids, data), returning Spiner::InterpolationHelpers::Result1D<T> and Result2D<T> respectively. grids and data are passed by const reference; grid index zero is the fastest axis, and data are accessed as data(i1) or data(i2, i1). Policies must support the calling execution space and must not throw, since the DataBox entry points are noexcept. DataBox retains the rank and grid validity assertions before dispatching to the policy. The policy type is available as DataBox::GradStencilType and is preserved by slices, copies, and device transfers. It does not change value-only interpToReal calls. The internal Concept template parameter follows GradStencil as the fourth parameter; code explicitly supplying the former third parameter must move it to fourth position.

Note

In C++17 and later, you can also get the default type specialization by simply omitting the template arguments.

Note

In the function signatures below, GPU/performance portability decorators have been excluded for brevity. However they are present in the actual code.

Note

In the function signatures below, we will often refer to the type Real and the type T. These are both references to the underlying templated arithmetic type.

Creating a `DataBox`

You can create a DataBox of a given shape via the constructor:

int nx1 = 2;
int nx2 = 3;
int nx3 = 4;
Spiner::DataBox<double> db(nx3, nx2, nx1);

The constructor takes any number of shape values (e.g., nx*) up to six (or Spiner::MAXRANK) values. Zero shape values initializes an empty, size-zero array.

Note

DataBox is row-major ordered. By convention, x3 is the slowest moving index and x1 is the fastest.

If GPU support is enabled, a DataBox can be allocated on either host or device, depending on the AllocationTarget. For example, to explicitly allocate one array on the host and one on the device, you might call:

// Allocates on the host (CPU)
Spiner::DataBox<double> db_host(Spiner::AllocationTarget::Host, nx2, nx1);
// Allocates on the device (GPU)
Spiner::DataBox<double> db_dev(Spiner::AllocationTarget::Device, nx2, nx1);

Note

If GPU support is not enabled, these both allocate on host.

You can also wrap a DataBox around a pointer you allocated yourself. For example:

std::vector<double> mydata(nx1*nx2);
Spiner::DataBox<double> db(mydata.data(), nx2, nx1);

You can also resize a DataBox, which you can use to modify a DataBox in-place. For example:

Spiner::DataBox<double> db; // empty
// clears old memory, resizes the underlying array,
// and resets strides
db.resize(nx3, nx2, nx1);

Just like the constructor, resize takes an optional (first) argument for the AllocationTarget.

Warning

DataBox::resize is destructive. The underlying data is not preserved.

If you want to change the stride without changing the underlying data, you can use reshape, which modifies the dimensions of the array, without modifying the underlying memory. For example:

// allocate a 1D databox
Spiner::DataBox<double> db(nx3*nx2*nx1);
// interpret it as a 3D object
db.reshape(nx3, nx2, nx1);

Warning

Make sure not to change the underlying size of the array when using reshape. This is checked with an assert statement, so you will get errors when compiling without the NDEBUG preprocessor macro.

The method

void DataBox::reset();

sets the DataBox to be empty with zero rank.

Copying a `DataBox` to device

If GPU support is enabled, you can deep-copy a DataBox, its tabulated values, and (by default) its grids from host to device with the function

DataBox getOnDeviceDataBox(const DataBox &db_host, bool include_grids = true);

which returns a new databox with the data in db_host copied to GPU. An object-oriented method

DataBox DataBox::getOnDevice(bool include_grids = true) const;

exists as well, which returns a new object with the underlying data copied to GPU.

Set include_grids to false to skip the deep copy of the grid objects. This permits manual management of grid memory when reducing the device-memory footprint is important. The result holds non-owning shallow copies of the source’s grids (see shallowCopy in the grid documentation), so finalizing them does not free the source’s grid memory. Its grids must be made device-accessible before they are used for device-side interpolation.

Note

If GPU support is not enabled, getOnDevice and friends are no-ops.

Semantics and Memory Management

DataBox has reference semantics—meaning that copying a DataBox does not copy the underlying data. In other words,

Spiner::DataBox<double> db1(size);
Spiner::DataBox<double> db2 = db1;

shallow-copies db1 into db2. Especially for Kokkos like workflows, this is very useful.

Warning

DataBox is neither reference-counted nor garbage-collected. If you create a DataBox you must clear the memory allocated just like you would for a pointer.

Two functions are provided for freeing memory in DataBox:

void free(DataBox &db);

and

DataBox::finalize();

both will do the same thing and free the values in a DataBox in a context-dependent way. I.e., no matter what the AllocationTarget was, the appropriate memory will be freed.

Some grids, such as NonUniformGrid1D, own dynamic memory of their own. That memory is managed separately from the values and is not freed by free or finalize. Free it with

void DataBox::finalizeGrids();

which finalizes every grid held by the DataBox. Grids that do not own memory, including all RegularGrid1D and PiecewiseGrid1D grids and non-owning shallow copies, are left alone, so it is always safe to call. It is independent of finalize and may be called on a DataBox that does not own its values. For example:

Spiner::DataBox<double, Spiner::NonUniformGrid1D<double>> db(N);
db.setRange(0, points); // db's grid owns a copy of points
// ...
db.finalizeGrids(); // frees the grid coordinates
free(db);           // frees the values

Grid ownership follows one rule. An operation that makes a new handle to the same DataBox copies its grids along with their ownership. An operation that builds a different DataBox from another’s metadata either gives it non-owning shallow copies of the grids or, if it is a deep copy, new grids that it owns. It never makes two DataBox objects owners of the same grid memory. Specifically:

setRange copies the grid it is given, including its ownership, so a DataBox given a memory-owning grid is responsible for it.
Copy construction and assignment of a DataBox produce handles to the same DataBox. They report the same ownership as the source, for both values and grids. Free through only one of them.
Slices, copyShape, and getOnDevice(false) build a different DataBox and give it non-owning shallow copies of the grids.
getOnDevice(), copyMetadata, and deepCopy build a different DataBox with new grids that it owns.

Warning

Do not free a DataBox if its memory is managed externally, e.g., via a std::vector. DataBox checks for this use-case via an assert statement.

You can check whether a given DataBox is empty, unmanaged, or allocated on host or device with the

DataBox::dataStatus() const;

method. It returns an enum class, Spiner::DataStatus, which can take on the values Empty, Unmanaged, AllocatedHost, or AllocatedDevice. You can also check whether or not free should be called with the method

bool DataBox::ownsAllocatedMemory();

which returns true if a given databox is managing memory and false otherwise. The method

bool DataBox::isReference();

returns false if the databox is managing memory and true otherwise.

To make an independent copy of a host-accessible DataBox, use

void DataBox::deepCopy(const DataBox &src);

which allocates new memory and copies the values, index types, and grids of src. The grids are deep copies owned by the new DataBox, so free it with both finalizeGrids and finalize. The grids of the DataBox being copied into must not own memory. Use getOnDevice to copy a DataBox to device.

Using `DataBox` with smart pointers

Smart pointers can be used to manage a DataBox and automatically call free for you, so long as you use them with a custom deleter. Spiner provides the following deleter for use in this scenario:

template <bool FinalizeGrids = false>
struct DBDeleter {
  template <typename T>
  void operator()(T *ptr) {
    if constexpr (FinalizeGrids) {
      ptr->finalizeGrids();
    }
    ptr->finalize();
    delete ptr;
  }
};

DBDeleter<> frees only the values, like free. DBDeleter<true> also calls finalizeGrids, which you need when the DataBox owns memory-owning grids.

It can be used, for example, with a std::unique_ptr via:

// needed for smart pointers
#include <memory>

// Creates a unique pointer pointing to a DataBox
// with memory allocated on device
std::unique_ptr<DataBox, Spiner::DBDeleter<>> pdb(
  new DataBox(Spiner::AllocationTarget::Device, N));

// Before using the databox in, e.g., a GPU or Kokkos kernel, get a
// shallow copy:
auto db = *pdb;
// some kokkos code...

// when you leave scope, the data box will be freed.

Serialization and de-serialization

Shared memory models, such as MPI Windows, require allocation of memory through an external API call (e.g., MPI_Win_allocate_shared), which tabulated data must be written to. Spiner supports this model through serialization and de-serialization. The relevant methods are as follows. The function

std::size_t DataBox::serializedSizeInBytes() const;

reports how much memory a DataBox object requires to be externally allocated. The function

std::size_t DataBox::dynamicMemorySizeInBytes() const;

reports the size of the value buffer and grid-owned dynamic memory, excluding the inline DataBox object representation. The function

std::size_t DataBox::dumpDynamicMemory(std::byte *dst) const;

writes only that dynamic memory. It is used internally by serialize after the inline object bytes have been copied.

std::size_t serialize(std::byte *dst) const;

takes a std::byte* pointer, assumed to contain enough space for a DataBox, and stores all information needed for the DataBox to reconstruct itself. The return value is the amount of memory in bytes used in the array by the serialized DataBox object. This method is non-destructive; the original DataBox is unchanged. The function

std::size_t DataBox::setPointer(T *src);

with the overload

std::size_t DataBox::setPointer(std::byte *src);

sets the underlying tabulated data from the src pointer, which is assumed to be the right size and shape, and relocates the trailing grid-owned dynamic memory. This is useful for the deSerialize function (described below) and for building your own serialization/de-serialization routines in composite objects. The function

std::size_t DataBox::deSerialize(std::byte *src);

initializes a DataBox to match the serialized DataBox contained in the src pointer.

Note

Note that the de-serialized DataBox has unmanaged memory, as it is assumed that the src pointer manages its memory for it. Therefore, one cannot free the src pointer until everything you want to do with the de-serialized DataBox is over.

The serialized memory is laid out as the inline DataBox metadata, followed by the tabulated values and grid-owned dynamic memory. The inline grid objects are already part of the DataBox metadata, so their static bytes are not emitted a second time. Dynamic memory is dumped for every dimension, including indexed and named dimensions, so resource management does not depend on mutable interpolation metadata.

The source buffer must remain alive and suitably aligned for the entire lifetime of every DataBox de-serialized from it.

Warning

DataBox previously exposed the serialize and deserialize routines using the C-style char* pointer. spiner now prefers the more std::byte* pointer. The former is still supported via overloads but should be considered deprecated.

Putting this all together, an application of serialization/de-serialization probably looks like this:

// load a databox from, e.g., file
Spiner::DataBox<double> db;
db.loadHDF(filename);

// get size of databox
std::size_t allocate_size = db.serializedSizeInBytes();

// Allocate the memory for the new databox.
// In practice this would be an API call for, e.g., shared memory
std::byte *memory = (std::byte*)malloc(allocate_size);

// serialize the old databox
std::size_t write_size = db.serialize(memory);

// make a new databox and de-serialize it
Spiner::DataBox<double> db2;
std::size_t read_size = db2.deSerialize(memory);

// read_size, write_size, and allocate_size should all be the same.
assert((read_size == write_size) && (write_size == allocate_size));

Warning

DataBox serialization is intended for transient communication between processes using the same Spiner version and a compatible architecture. The binary representation is not guaranteed to be compatible across Spiner versions, architectures, compilers, or build configurations. This is due to, for example, differences in endianness, alignment, and padding.

Do not use serialization as a persistent-storage or web-interchange format. Use HDF5 for persistent, portable storage. Size calculation, serialization, and de-serialization must all use a compatible Spiner build.

Accessing Elements of a `DataBox`

Elements of a DataBox can be accessed and set via the () operator. For example:

Spiner::DataBox<double> db(nx3, nx2, nx1);
db(2,1,0) = 5.0;

The () operator accepts between one and six indexes. If you pass in more indexes than the rank of the array, the excess indices are ignored. If you pass in fewer, the unset indices are assumed to be zero. The exception is the one-dimensional operator. You can always stride through the “flattened” array by using the one-dimensional accessor. For example:

for (int i = 0; i < nx3*nx2*nx1; ++i) {
  db(i) = static_cast<double>(i);
}

fills the three-dimensional array above with the flat index of each element.

Slicing

A new DataBox containing a shallow slice of another DataBox can be constructed with the slice method:

DataBox DataBox::slice(const int dim, const int indx, const int nvar) const;

this is fairly limited functionality. It returns a new DataBox containing only elements from indx to indx + nvar - 1 in the dim direction. All other directions are unchanged. The slowest moving dimension can be sliced to a single index with

DataBox DataBox::slice(const int indx) const;

and the slowst-moving two dimensions can be sliced to a single pair of indicies with

DataBox DataBox::slice(const int i2, int i1) const;

Index Types and Interpolation Ranges

Often-times an array mixes “continuous” and “discrete” variables. In other words, some indices of an array are discretizations of a continuous quantity, and we want to interpolate in those directions, but other indices are discrete—they may index a particle species, for example. A common example is in neutrino transport, where an array of emissivities may depend on fluid density, fluid temperature, electron fraction, neutrino energy, and neutrino species. The species can only take three discrete values, but the density, temperature, and electron fraction are all continuous.

Spiner accounts for this by assigning each dimension in the array a “type,” represented as an enum class, IndexType. Currently the type can be either Interpolated or Indexed. When a new DataBox is created, all dimensions are set to IndexType::Indexed. A dimension can be set to Interpolated via the setRange method.

void DataBox::setRange(int i, Grid_t g);

where here i is the dimension and g is the gridding object for this index. In the default setup, where grids are uniformly spaced (i.e., you use a RegularGrid1D), this is:

void DataBox::setRange(int i, T min, T max, int N);

where here i is the dimension, min is the minimum value of the independent variable, max is the maximum value of the indpendent variable, and N is the number of points in the i dimension. (Here T is the underlying templated data type.)

Warning

In these routines, the dimension is indexed from zero. Moreover, the indexing possibly the opposite direction you’re expecting. The zeroth dimension is the right most index when interpToReal is called.

Note

There is a set of lower-level objects, RegularGrid1D, and PiecewiseGrid1D, which represent these interpolation ranges internally. There is a getter method range that works with the underlying Grid_t class directly. For more details, see the relevant documentation.

It’s often desirable to have multiple databoxes with the exact same shape and interpolation structure (i.e., independent variable ranges). In this case, the method

void DataBox::copyMetadata(const DataBox &src);

can assist. This method resets and re-allocates the data in a DataBox to the exact same size and shape as src. More importantly, it also copies the relevant IndexType and independent variable range for each dimension. The grids are deep copies owned by the DataBox, so free them with finalizeGrids. The grids of the DataBox being reset must not own memory.

One can also manually set the IndexType in a given dimension with

void DataBox::setIndexType(int i, IndexType t);

and retrieve the IndexType with

IndexType &DataBox::indexType(const int i);

to see if a dimension is interpolatable.

Interpolation to a real number

The family of DataBox::interpToReal methods interpolate the “entire” DataBox to a real number. Up to four-dimensional interpolation is supported:

T DataBox::interpToReal(const T x) const;

T DataBox::interpToReal(const T x2, const T x1) const;

T DataBox::interpToReal(const T x3, const T x2, const T x1) const;

T DataBox::interpToReal(const T x4, const T x3, const T x2, const T x1) const;

where x1 is the fastest moving direction, x2 is less fast, and so on. These interpolation routines are hand-tuned for performance.

Warning

Do not call interpToReal with a DataBox that is the wrong shape or try to interpolate on indices that are not interpolatable. This is checked with an assert statement.

Interpolation with gradients

interpToRealWithGrads returns the interpolated value and its derivatives with respect to the continuous grid coordinates. It supports fully interpolated rank-1 and rank-2 databoxes:

DataBox::InterpResult1D DataBox::interpToRealWithGrads(const T x1) const noexcept;

DataBox::InterpResult2D DataBox::interpToRealWithGrads(const T x2, const T x1) const noexcept;

The nested result aliases name plain aggregates with members of type T:

struct InterpResult1D { T value, d_dx1; };
struct InterpResult2D { T value, d_dx2, d_dx1; };

Use named members or structured bindings. As with all databox types, arguments and returned derivatives follow slowest-to-fastest order: x2 then x1. For example, this host-side table represents f(y, x) = 1 + 2*x - 3*y + 4*x*y:

Spiner::DataBox<double> db(3, 4); // (ny, nx)
db.setRange(0, 0.0, 3.0, 4);    // x, fastest dimension
db.setRange(1, 0.0, 4.0, 3);    // y, slowest dimension
for (int j = 0; j < 3; ++j) {
  const double y = 2.0 * j;
  for (int i = 0; i < 4; ++i) {
    const double x = i;
    db(j, i) = 1.0 + 2.0*x - 3.0*y + 4.0*x*y;
  }
}
auto [value, df_dy, df_dx] = db.interpToRealWithGrads(1.0, 0.5);
// value = 1, df_dy = -1, df_dx = 6
db.finalize();

For a fully interpolated 1D databox, use:

auto [value, df_dx] = db1.interpToRealWithGrads(x);
// Equivalently:
auto result = db1.interpToRealWithGrads(x);
// result.value and result.d_dx1

The derivatives are analytic derivatives of the selected policy’s interpolant. With the default Linear policy, this uses the same cell selection and weights as interpToReal. They are not interpolated samples of a separately computed gradient field. At knots or cell boundaries, derivatives can be discontinuous; the result uses the derivative from the cell selected by the grid lookup. With Linear, out-of-range coordinates use the boundary cell’s linear extrapolation and its derivatives.

Derivatives are taken with respect to the supplied grid coordinates. If a coordinate represents a transformed variable, such as a logarithm, apply the chain rule to obtain derivatives with respect to the original variable.

All four built-in grid types support these methods. Custom grid types must provide weightsWithGrad to support this mode for linear interpolants. The databox must have the matching rank and all dimensions must be interpolatable with well-formed grids, as checked by assertions. There are no mixed indexed/interpolated or higher-dimensional overloads of interpToRealWithGrads.

Rational function interpolation

Spiner::InterpolationStencils::Rational implements the rational function interpolation described by Kerley (1977, LA-6903-MS), with analytic coordinate derivatives. Select it as the third DataBox template parameter:

using RationalBox = Spiner::DataBox<
    double, Spiner::RegularGrid1D<double>,
    Spiner::InterpolationStencils::Rational>;

// Host-side example: four samples of f(x) = x*x.
RationalBox db(4);
db.setRange(0, 0.0, 3.0, 4);
for (int i = 0; i < 4; ++i) {
  db(i) = double(i) * double(i);
}
auto [value, df_dx] = db.interpToRealWithGrads(1.5);
// value = 2.25, df_dx = 3.0
db.finalize();

The policy supports fully interpolated rank-1 and rank-2 databoxes and requires at least four points on every interpolated axis. All four built-in grid types are supported. For 2D, call db.interpToRealWithGrads(y, x); the result contains value, d_dx2 (the y derivative), and d_dx1 (the x derivative), in that order.

The one-dimensional stencil uses rational interpolation in interior intervals and quadratic formulas in the first and last intervals. Strictly outside the table, it uses linear secant extrapolation from the nearest two endpoint samples. At an exact table endpoint, it returns the interior derivative, which can differ from the exterior secant slope; derivatives can therefore jump at the endpoints. These axis rules also apply in 2D.

The policy affects only interpToRealWithGrads. Value-only interpToReal calls remain linear, even on a RationalBox. To obtain only the rationally interpolated value, use db.interpToRealWithGrads(x).value (or the two-coordinate overload).

Mixed interpolation and indexing

In the case where an array has some dimensions that are discrete and some that are interpolatable, one can fuse interpolation and indexing into a single operation. These operations are still named DataBox::interpToReal, but one of the input arguments is an integer instead of a floating point number. The location of the integer in the function signature indicates which dimension in the DataBox is indexed. For example:

T DataBox::interpToReal(const T x3, const T x2, const T x1, const int idx) const;

interpolates the three slower-moving indices and indexes the fastest moving index. On the other hand,

T DataBox::interpToReal(const T x4, const T x3, const T x2, const int idx, const T x1) const;

interpolates the fastest moving index, then indexes the second-fastest, then interpolates the remaining three slower. The above fused operations are the only ones currently supported.

Interpolating into another `DataBox`

There is limited functionality for filling a DataBox with the interpolated values of another DataBox. For example, the method

void DataBox::interpFromDB(const DataBox &src, const T x);

fills the DataBox with the faster moving elements of src interpolated to x in the slowest-moving direction. Similarly,

void DataBox::interpFromDB(const DataBox &src, const T x2, const T x1);

interpolates the two slowest-moving directions of src to x2 and x1.

Warning

interpFromDB never allocates, reshapes, or frees memory. The DataBox being filled must already exist with storage and the shape of the remaining faster-moving dimensions of src. In other words, it must have rank src.rank() - 1 (or src.rank() - 2) and dim(i) == src.dim(i) for each remaining dimension. Both DataBox objects must be accessible from the execution space in which interpFromDB is called, for example both on device inside a kernel.

Only values are written. The grids and index types of the DataBox being filled are left unchanged, so call setRange on it yourself if you want to interpolate it later. For example:

// src has shape (NZ, NY, NX)
Spiner::DataBox<double> dst(NY, NX);
dst.interpFromDB(src, z);
dst.setRange(1, src.range(1));
dst.setRange(0, src.range(0));
double val = dst.interpToReal(y, x);
// ...
free(dst);

Warning

Be careful when the grids own their own memory, as NonUniformGrid1D does. setRange copies a grid’s ownership along with it, so passing src.range(i) would make both DataBox objects owners of the same grid memory. To share the grids, pass a non-owning handle instead:

Spiner::NonUniformGrid1D<double> x_grid;
x_grid.shallowCopy(src.range(0));
dst.setRange(0, x_grid);

Then only src frees the grid memory, through finalizeGrids.

File I/O

If hdf5 is enabled, Spiner can save an array to or load an array from disk. Each array so-saved is also saved with the IndexType and independent variable ranges bundled with it, so that knowledge of how to interpolate the data is automatically available.

The following methods are supported:

herr_t DataBox::saveHDF(const std::string &filename) const;

saves the DataBox to a file with filename.

herr_t DataBox::saveHDF(hid_t loc, const std::string &groupname) const;

saves the DataBox as an hdf5 group at the location loc in an hdf5 file.

DataBox::loadHDF(const std::string &filename);

fills the DataBox from information in the root of a file with filename.

DataBox::loadHDF(hid_t loc, const std::string &groupname);

fills the DataBox from information in the group with groupname based at location loc in the file.

Warning

HDF5 I/O is only supported for single- and double-precision types at this time.

Miscellany

Here we list a few convenience functions available that were not covered elsewhere.

T DataBox::min() const;

and

T DataBox::max() const;

compute and return the minimum and maximum values (respectively) in the array.

int rank() const;

returns the rank (number of dimensions) of the array.

int size() const;

returns the total number of elements in the underlying array.

int sizeBytes() const;

returns the total size of the underlying array in bytes.

int dim(int i) const;

returns the size in a given dimension/direction, indexed from one.

Generative AI was used to assist with modifications to this page.

The DataBox

Creating a DataBox

Copying a DataBox to device

Semantics and Memory Management

Using DataBox with smart pointers

Serialization and de-serialization

Accessing Elements of a DataBox

Slicing

Index Types and Interpolation Ranges

Interpolation to a real number

Interpolation with gradients

Rational function interpolation

Mixed interpolation and indexing

Interpolating into another DataBox

File I/O

Miscellany

Creating a `DataBox`

Copying a `DataBox` to device

Using `DataBox` with smart pointers

Accessing Elements of a `DataBox`

Interpolating into another `DataBox`