DDC 0.0.0

a discrete domain computation library

Commented example: the uniform heat equation

In examples/uniform_heat_equation.cpp is a DDC example implementing a forward finite-difference solver for the heat equation over a rectangle 2D domain with periodic boundary conditions.

As usual, the file starts with a few includes that will be used in the code.

#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <ddc/ddc.hpp>

As you can see, to use DDC, we have to include <ddc/ddc.hpp>

Definition of the discretization

Before solving the equation, DDC's primary goal is to define a discrete domain of dimensions specified by the user, along with a discretization along each dimension, which also needs to be specified.

Domains

Each point in the DiscreteDomain is a DiscreteElement. These concepts will be clarified later. Let's start by constructing this DiscreteDomain necessary for solving our 2D problem for the heat equation.

Dimensions naming

We start by defining types that we later use to name our dimensions.

First, we create X: a type that is declared but never defined to act as a tag identifying our first dimension.

struct X;

Then, we create DDimX, a type that will also act as a tag, but we define it to be a uniform discretization of X. Thus DDimX is an identifier for a discrete dimension.

struct DDimX : ddc::UniformPointSampling<X>
{
};
UniformPointSampling models a uniform discretization of the provided continuous dimension.
Definition uniform_point_sampling.hpp:36

We do the same thing for the second dimension Y.

struct Y;
struct DDimY : ddc::UniformPointSampling<Y>
{
};

And once again, now for the time dimension.

struct T;
struct DDimT : ddc::UniformPointSampling<T>
{
};

Domains

Dimension X

Once the types are defined, we can start the main function where we will define our various domains. Here for each dimension, the user needs to specify the starting and ending coordinates of the dimension of the domain, as well as the number of discretization points along each of these dimensions. Additionally, we specify here the physical characteristics specific to our equation (the thermal diffusion coefficient).

int main(int argc, char** argv)
{
#ifdef DDC_BUILD_PDI_WRAPPER
auto pdi_conf = PC_parse_string("");
PDI_init(pdi_conf);
#endif
Kokkos::ScopeGuard const kokkos_scope(argc, argv);
ddc::ScopeGuard const ddc_scope(argc, argv);
double const x_start = -1.;
double const x_end = 1.;
std::size_t const nb_x_points = 10;
double const kx = .01;
Definition scope_guard.hpp:17

We start by defining gwx, the number of "ghost" points in the X dimension, those represented in dark grey in the figure on each side of the domain.

A DiscreteVector is a vector in the discrete dimension.
Definition discrete_vector.hpp:254

The type for this constant is DiscreteVector<DDimX> that represents a number of elements in the discretization of the X dimension.

Once done, we initialize the discretization of the X dimension. Its name has been specified before, but we now set its parameters with the init_discretization function.

auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
= ddc::init_discrete_space<DDimX>(DDimX::init_ghosted<DDimX>(
gwx));
void init_discrete_space(Args &&... args)
Initialize (emplace) a global singleton discrete space.
Definition discrete_space.hpp:151
detail::TaggedVector< CoordinateElement, CDims... > Coordinate
A Coordinate represents a coordinate in the continuous space.
Definition coordinate.hpp:26

Depending on the way the function is called, its return type can differ. Here we use it with an inner call to init_ghosted and receive four 1D domains as a result. Their type is not specified because we use C++ structured bindings, but they are all of the same type: DiscreteDomain<DDimX> that represents a set of contiguous points in the discretization of X.

init_ghosted takes as parameters the coordinate of the first and last discretized points, the number of discretized points in the domain and the number of additional points on each side of the domain. The fours DiscreteDomains returned are:

  • x_domain: the main domain from the start (included) to end (included) but excluding "ghost" points.
  • ghosted_x_domain: the domain including all "ghost" points.
  • x_pre_ghost: the "ghost" points that come before the main domain.
  • x_post_ghost: the "ghost" points that come after the main domain.

The parameters of raw C++ types like double or size_t can not be used as-is since DDC enforces strong typing. Instead, for a coordinate in the X dimension, we use Coordinate<X> and –as already mentioned– for a number of elements in the discretization of X, we use DiscreteVector<DDimX>.

Once this is done, we define two additional domains that will be mirrored to the ghost at the start and at the end of the domain.

x_domain_begin(x_domain.front(), x_post_ghost.extents());
ddc::DiscreteDomain<DDimX> const x_domain_end(
x_domain.back() - x_pre_ghost.extents() + 1,
x_pre_ghost.extents());
Definition discrete_domain.hpp:51

x_domain_begin is the sub-domain at the beginning of x_domain of the same shape as x_post_ghost that will be mirrored to it. Reciprocally for x_domain_end with x_pre_ghost.

The type of both sub-domains is DiscreteDomain<DDimX> even if only DiscreteDomain is specified, this relies on C++ class template argument deduction (CTAD). The first parameter given to the constructor is the first element of the domain, a DiscreteElement<DDimX>, the second parameter is a number of elements to include, a DiscreteVector<DDimX>.

To summarize, in this section, we have introduced the following types:

  • Coordinate<X> that represents a point in the continuous dimension X,
  • DiscreteElement<DDimX> that represents one of the elements of DDimX, the discretization of X,
  • DiscreteVector<DDimX> that represents a number of elements of DDimX,
  • DiscreteDomain<DDimX> that represents an interval in DDimX, a set of contiguous elements of the discretization.

Dimension Y

The domains in the Y dimension are handled in a way very similar to the X dimension. We first define the domain characteristics along this dimension

double const y_start = -1.;
double const y_end = 1.;
std::size_t const nb_y_points = 100;
double const ky = .002;

Then we initialize the domain along this dimension just like we did with the X dimension.

auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
= ddc::init_discrete_space<DDimY>(DDimY::init_ghosted<DDimY>(
gwy));
y_domain_begin(y_domain.front(), y_post_ghost.extents());
ddc::DiscreteDomain<DDimY> const y_domain_end(
y_domain.back() - y_pre_ghost.extents() + 1,
y_pre_ghost.extents());

Time dimension

Then we handle the domains for the simulated time dimension. We first give the simulated time at which to stard and end the simulation.

double const start_time = 0.;
double const end_time = 10.;

Then we use the CFL condition to determine the time step of the simulation.

double const dx = ddc::step<DDimX>();
double const dy = ddc::step<DDimY>();
double const invdx2 = 1. / (dx * dx);
double const invdy2 = 1. / (dy * dy);
ddc::Coordinate<T> const dt(.5 / (kx * invdx2 + ky * invdy2));

Finally, we determine the number of time steps and as we did with the X and Y dimensions, we create the time domain.

ddc::DiscreteVector<DDimT> const nb_time_steps(
std::ceil((end_time - start_time) / dt) + .2);
ddc::DiscreteDomain<DDimT> const time_domain
= ddc::init_discrete_space<DDimT>(DDimT::init<DDimT>(
ddc::Coordinate<T>(start_time),
ddc::Coordinate<T>(end_time),
nb_time_steps + 1));

Data allocation

We allocate two 2D Chunks along the X and Y dimensions which will be used to map temperature to the domains' points at t and t+dt. When constructing the Chunks one can give an optional string to label the memory allocations. This helps debugging and profiling applications using the Kokkos tools, see also Kokkos Tools.

These chunks map the temperature into the full domain (including ghosts) twice:

  • ghosted_last_temp for the last fully computed time-step.
  • ghosted_next_temp for time-step being computed.
ddc::Chunk ghosted_last_temp(
"ghosted_last_temp",
DDimX,
DDimY>(ghosted_x_domain, ghosted_y_domain),
ddc::Chunk ghosted_next_temp(
"ghosted_next_temp",
DDimX,
DDimY>(ghosted_x_domain, ghosted_y_domain),
Definition kokkos_allocator.hpp:17
Definition chunk.hpp:21

Note that the DeviceAllocator is responsible for allocating memory on the default memory space.

Initial conditions

To set the initial conditions, the ghosted_intial_temp is created and acts as a pointer to the chunk. The const qualifier makes it clear that ghosted_initial_temp always references the same chunk, ghosted_last_temp in this case.

ddc::ChunkSpan const ghosted_initial_temp
= ghosted_last_temp.span_view();
Definition chunk_span.hpp:30

Then, we iterate over each DiscreteElement of the domain to fill ghosted_initial_temp with the initial values of the simulation.

ddc::DiscreteDomain<DDimX, DDimY>(x_domain, y_domain),
KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> const ixy) {
double const x = ddc::coordinate(
double const y = ddc::coordinate(
ghosted_initial_temp(ixy)
= 9.999 * ((x * x + y * y) < 0.25);
});
A DiscreteElement identifies an element of the discrete dimension.
Definition discrete_element.hpp:146
KOKKOS_FUNCTION Coordinate< typename DDim::continuous_dimension_type... > coordinate(DiscreteElement< DDim... > const &c)
Definition coordinate.hpp:29
void parallel_for_each(std::string const &label, ExecSpace const &execution_space, DiscreteDomain< DDims... > const &domain, Functor &&f) noexcept
iterates over a nD domain using a given Kokkos execution space
Definition parallel_for_each.hpp:155

To display the data, a chunk is created on the host.

ddc::Chunk ghosted_temp(
"ghost_temp",
DDimX,
DDimY>(ghosted_x_domain, ghosted_y_domain),

We deepcopy the data from the ghosted_last_temp chunk to ghosted_temp on the host.

ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp);
auto parallel_deepcopy(ChunkDst &&dst, ChunkSrc &&src)
Copy the content of a borrowed chunk into another.
Definition parallel_deepcopy.hpp:21

And we display the initial data.

display(ddc::coordinate(time_domain.front()),
ghosted_temp[x_domain][y_domain]);
KOKKOS_FUNCTION constexpr discrete_element_type front() const noexcept
Definition discrete_domain.hpp:139
for (ddc::DiscreteElement<DDimT> const iter :
time_domain.remove_first(ddc::DiscreteVector<DDimT>(1))) {
The top-level namespace of DDC.
Definition aligned_allocator.hpp:11

To display the data, a chunk is created on the host.

ddc::Chunk ghosted_temp(
"ghost_temp",
DDimX,
DDimY>(ghosted_x_domain, ghosted_y_domain),

We deepcopy the data from the ghosted_last_temp chunk to ghosted_temp on the host.

ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_pre_ghost, y_domain)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(y_domain, x_domain_end)]);
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(y_domain, x_post_ghost)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(y_domain, x_domain_begin)]);
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_pre_ghost)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain_end)]);
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_post_ghost)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain_begin)]);
ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp);

And we display the initial data.

display(ddc::coordinate(time_domain.front()),
ghosted_temp[x_domain][y_domain]);

For the numerical scheme, two chunkspans are created:

  • next_temp a span excluding ghosts of the temperature at the time-step we will build.
  • last_temp a read-only view of the temperature at the previous time-step.Note that span_cview returns a read-only ChunkSpan.
ddc::ChunkSpan const next_temp(
ghosted_next_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain)]);
ddc::ChunkSpan const last_temp(ghosted_last_temp.span_cview());

We then solve the equation.

next_temp.domain(),
KOKKOS_LAMBDA(
double const dt = ddc::step<DDimT>();
next_temp(ix, iy) = last_temp(ix, iy);
next_temp(ix, iy) += kx * dt
* (last_temp(ix + 1, iy)
- 2.0 * last_temp(ix, iy)
+ last_temp(ix - 1, iy))
* invdx2;
next_temp(ix, iy) += ky * dt
* (last_temp(ix, iy + 1)
- 2.0 * last_temp(ix, iy)
+ last_temp(ix, iy - 1))
* invdy2;
});

Time loop

for (ddc::DiscreteElement<DDimT> const iter :
time_domain.remove_first(ddc::DiscreteVector<DDimT>(1))) {

Periodic conditions

if (iter - last_output_iter >= t_output_period) {
last_output_iter = iter;
ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp);
display(ddc::coordinate(iter),
ghosted_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain)]);
}
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_pre_ghost, y_domain)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(y_domain, x_domain_end)]);
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(y_domain, x_post_ghost)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(y_domain, x_domain_begin)]);
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_pre_ghost)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain_end)]);
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_post_ghost)],
ghosted_last_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain_begin)]);

Numerical scheme

std::swap(ghosted_last_temp, ghosted_next_temp);

For the numerical scheme, two chunkspans are created:

  • next_temp a span excluding ghosts of the temperature at the time-step we will build.
  • last_temp a read-only view of the temperature at the previous time-step.Note that span_cview returns a read-only ChunkSpan.
ddc::ChunkSpan const next_temp(
ghosted_next_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain)]);
ddc::ChunkSpan const last_temp(ghosted_last_temp.span_cview());

We then solve the equation.

if (last_output_iter < time_domain.back()) {
ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp);
display(ddc::coordinate(time_domain.back()),
ghosted_temp[ddc::DiscreteDomain<
DDimX,
DDimY>(x_domain, y_domain)]);
}
KOKKOS_FUNCTION constexpr discrete_element_type back() const noexcept
Definition discrete_domain.hpp:144
next_temp.domain(),
KOKKOS_LAMBDA(
double const dt = ddc::step<DDimT>();
next_temp(ix, iy) = last_temp(ix, iy);
next_temp(ix, iy) += kx * dt
* (last_temp(ix + 1, iy)
- 2.0 * last_temp(ix, iy)
+ last_temp(ix - 1, iy))
* invdx2;
next_temp(ix, iy) += ky * dt
* (last_temp(ix, iy + 1)
- 2.0 * last_temp(ix, iy)
+ last_temp(ix, iy - 1))
* invdy2;
});