DDC 0.0.0

a discrete domain computation library

Commented example: the heat equation

In examples/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>
#include <Kokkos_Core.hpp>

As you can see, DDC includes all follow the same convention: <ddc/SYMBOL> where SYMBOL is a the name of a DDC symbol. So for example, in order to use a class named Chunk, you should include <ddc/Chunk> and to use a free function template named for_each, you should include <ddc/for_each>.

Then, we define the value of some parameters that would typically be read from some form of configuration file in a more realistic code.

// Start of the domain of interest in the X dimension
double const x_start = -1.;
// End of the domain of interest in the X dimension
double const x_end = 1.;
// Number of discretization points in the X dimension
size_t const nb_x_points = 10;
// Thermal diffusion coefficient
double const kx = .01;
// Start of the domain of interest in the Y dimension
double const y_start = -1.;
// End of the domain of interest in the Y dimension
double const y_end = 1.;
// Number of discretization points in the Y dimension
size_t const nb_y_points = 100;
// Thermal diffusion coefficient
double const ky = .002;
// Simulated time at which to start simulation
double const start_time = 0.;
// Simulated time to reach as target of the simulation
double const end_time = 10.;
// Number of time-steps between outputs
size_t const t_output_period = 10;

Definition of the discretization

Given the previous parameters, we will now define a 2D point discretization.

Dimensions naming

Before starting the main function, 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.

UniformPointSampling models a uniform discretization of the provided continuous dimension.
Definition: uniform_point_sampling.hpp:21

We do the same thing for the second dimension Y.

// Our second continuous dimension
struct Y;
// Its uniform discretization

And once again, now for the time dimension.

// Our simulated time dimension
struct T;
// Its uniform discretization

Domains

Dimension X

Once the types are defined, we can start the main function where we will define our various domains.

int main(int argc, char** argv)
{
ScopeGuard scope(argc, argv);
// some parameters that would typically be read from some form of
// configuration file in a more realistic code
// Start of the domain of interest in the X dimension
double const x_start = -1.;
// End of the domain of interest in the X dimension
double const x_end = 1.;
// Number of discretization points in the X dimension
size_t const nb_x_points = 10;
// Thermal diffusion coefficient
double const kx = .01;
// Start of the domain of interest in the Y dimension
double const y_start = -1.;
// End of the domain of interest in the Y dimension
double const y_end = 1.;
// Number of discretization points in the Y dimension
size_t const nb_y_points = 100;
// Thermal diffusion coefficient
double const ky = .002;
// Simulated time at which to start simulation
double const start_time = 0.;
// Simulated time to reach as target of the simulation
double const end_time = 10.;
// Number of time-steps between outputs
size_t const t_output_period = 10;
Definition: scope_guard.hpp:12

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.

// Number of ghost points to use on each side in X
DiscreteVector<DDimX> static constexpr gwx {1};
A DiscreteVector is a vector in the discrete dimension.
Definition: discrete_vector.hpp:207

This can be made static constexpr since here this is a compile-time constant. 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.

// Initialization of the global domain in X with gwx ghost points on
// each side
auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
= init_discrete_space(DDimX::init_ghosted(
Coordinate<X>(x_start),
Coordinate<X>(x_end),
DiscreteVector<DDimX>(nb_x_points),
gwx));

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.

// our zone at the start of the domain that will be mirrored to the
// ghost
x_domain_begin(x_domain.front(), x_post_ghost.extents());
// our zone at the end of the domain that will be mirrored to the
// ghost
DiscreteDomain const x_domain_end(
x_domain.back() - x_pre_ghost.extents() + 1,
x_pre_ghost.extents());
Definition: discrete_domain.hpp:22

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.

// Number of ghost points to use on each side in Y
DiscreteVector<DDimY> static constexpr gwy {1};
// Initialization of the global domain in Y with gwy ghost points on
// each side
auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
= init_discrete_space(DDimY::init_ghosted(
Coordinate<Y>(y_start),
Coordinate<Y>(y_end),
DiscreteVector<DDimY>(nb_y_points),
gwy));
// our zone at the start of the domain that will be mirrored to the
// ghost
y_domain_begin(y_domain.front(), y_post_ghost.extents());
// our zone at the end of the domain that will be mirrored to the
// ghost
DiscreteDomain 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.

// max(1/dx^2)
double const invdx2_max = transform_reduce(
x_domain,
0.,
return 1.
/ (distance_at_left(ix) * distance_at_right(ix));
});
// max(1/dy^2)
double const invdy2_max = transform_reduce(
y_domain,
0.,
return 1.
/ (distance_at_left(iy) * distance_at_right(iy));
});
Coordinate<T> const max_dt {
.5 / (kx * invdx2_max + ky * invdy2_max)};
// number of time intervals required to reach the end time
DiscreteVector<DDimT> const nb_time_steps {
std::ceil((end_time - start_time) / max_dt) + .2};
// Initialization of the global domain in time:
// - the number of discrete time-points is equal to the number of
// steps + 1
DiscreteDomain<DDimT> const time_domain
= init_discrete_space(DDimT::
init(Coordinate<T>(start_time),
Coordinate<T>(end_time),
nb_time_steps + 1));
A DiscreteElement identifies an element of the discrete dimension.
Definition: discrete_element.hpp:123
Definition: reducer.hpp:96

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.

// Maps temperature into the full domain (including ghosts) twice:
// - once for the last fully computed time-step
Chunk ghosted_last_temp(
DDimX,
DDimY>(ghosted_x_domain, ghosted_y_domain),
// - once for time-step being computed
Chunk ghosted_next_temp(
DDimX,
DDimY>(ghosted_x_domain, ghosted_y_domain),
Definition: kokkos_allocator.hpp:12
Definition: chunk.hpp:12

Initial conditions

ChunkSpan const ghosted_initial_temp = ghosted_last_temp.span_view();
// Initialize the temperature on the main domain
for_each(
DiscreteDomain<DDimX, DDimY>(x_domain, y_domain),
DDC_LAMBDA(DiscreteElement<DDimX, DDimY> const ixy) {
double const x = coordinate(select<DDimX>(ixy));
double const y = coordinate(select<DDimY>(ixy));
ghosted_initial_temp(ixy)
= 9.999 * ((x * x + y * y) < 0.25);
});
Definition: chunk_span.hpp:24
constexpr parallel_device_policy parallel_device
Definition: for_each.hpp:167
// display the initial data
deepcopy(ghosted_temp, ghosted_last_temp);
display(coordinate(time_domain.front()),
ghosted_temp[x_domain][y_domain]);
// time of the iteration where the last output happened
DiscreteElement<DDimT> last_output = time_domain.front();
constexpr discrete_element_type front() const noexcept
Definition: discrete_domain.hpp:125

Time loop

for (auto const iter :
constexpr DiscreteDomain remove_first(mlength_type n) const
Definition: discrete_domain.hpp:145

Periodic conditions

// Periodic boundary conditions
deepcopy(
ghosted_last_temp[x_pre_ghost][y_domain],
ghosted_last_temp[y_domain][x_domain_end]);
deepcopy(
ghosted_last_temp[y_domain][x_post_ghost],
ghosted_last_temp[y_domain][x_domain_begin]);
deepcopy(
ghosted_last_temp[x_domain][y_pre_ghost],
ghosted_last_temp[x_domain][y_domain_end]);
deepcopy(
ghosted_last_temp[x_domain][y_post_ghost],
ghosted_last_temp[x_domain][y_domain_begin]);

Numerical scheme

// a span excluding ghosts of the temperature at the time-step we
// will build
ChunkSpan const next_temp {
ghosted_next_temp[x_domain][y_domain]};
// a read-only view of the temperature at the previous time-step
ChunkSpan const last_temp {ghosted_last_temp.span_view()};
// Stencil computation on the main domain
for_each(
next_temp.domain(),
DDC_LAMBDA(DiscreteElement<DDimX, DDimY> const ixy) {
DiscreteElement<DDimX> const ix = select<DDimX>(ixy);
DiscreteElement<DDimY> const iy = select<DDimY>(ixy);
double const dx_l = distance_at_left(ix);
double const dx_r = distance_at_right(ix);
double const dx_m = 0.5 * (dx_l + dx_r);
double const dy_l = distance_at_left(iy);
double const dy_r = distance_at_right(iy);
double const dy_m = 0.5 * (dy_l + dy_r);
next_temp(ix, iy) = last_temp(ix, iy);
next_temp(ix, iy)
+= kx * max_dt
* (dx_l * last_temp(ix + 1, iy)
- 2.0 * dx_m * last_temp(ix, iy)
+ dx_r * last_temp(ix - 1, iy))
/ (dx_l * dx_m * dx_r);
next_temp(ix, iy)
+= ky * max_dt
* (dy_l * last_temp(ix, iy + 1)
- 2.0 * dy_m * last_temp(ix, iy)
+ dy_r * last_temp(ix, iy - 1))
/ (dy_l * dy_m * dy_r);
});

Output

if (iter - last_output >= t_output_period) {
last_output = iter;
deepcopy(ghosted_temp, ghosted_last_temp);
display(coordinate(iter), ghosted_temp[x_domain][y_domain]);
}

Final swap

// Swap our two buffers
std::swap(ghosted_last_temp, ghosted_next_temp);

Final output

if (last_output < time_domain.back()) {
deepcopy(ghosted_temp, ghosted_last_temp);
display(coordinate(time_domain.back()),
ghosted_temp[x_domain][y_domain]);
}
constexpr discrete_element_type back() const noexcept
Definition: discrete_domain.hpp:130