Commented example: the non uniform heat equation

In examples/non_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 and non-uniform space discretization.

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

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <iostream>
#include <random>
#include <string>
#include <utility>
#include <vector>
 
#include <ddc/ddc.hpp>
 
#include <Kokkos_Core.hpp>

Differences with the uniform problem resolution

For non-uniform discretization, differences arise compared to uniform discretization, particularly in terms of dimension labeling, domain construction, and even in the resolution of the problem.

Dimensions naming

Just like the uniform case, we start by defining types that we later use to name our dimensions.

struct X;

However, we then define the discretization as non uniform.

struct DDimX : ddc::NonUniformPointSampling<X>
{
};

We do the same thing for the second dimension Y.

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

And for this case, we kept an uniform discretization for the time dimension.

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

Domains

Dimension X

Just like the uniform case, we define the main characteristics of the domain.

double const x_start = -1.;
double const x_end = 1.;
std::size_t const nb_x_points = 10;
double const kx = .01;

The first step to create a DiscreteDomain with non-uniform spatial discretization is to create a C++ iterator containing the actual coordinates of the points along each dimension. For tutorial purposes, we have constructed a function generate_random_vector that generates a vector with random data ranging from the lower bound to the higher one to populate our vectors.

For the X dimension, we want to build 4 domains:

• 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.

To do so, we have to create the iterator for the main domain.

std::vector<double> const x_domain_vect = generate_random_vector(nb_x_points, x_start, x_end);

For the initialization of ghost points in the non-uniform case, it was necessary to explicitly describe the position of the pre-ghost and post-ghost points. For the pre-ghost, it was necessary to start from the first coordinate along x and subtract the difference between the last and the penultimate coordinate of the x dimension to maintain periodicity. For the post-ghost point, take the last point and add the difference between the first and second points along the x dimension.

std::vector<double> const x_pre_ghost_vect = periodic_extrapolation_left(1, x_domain_vect);
std::vector<double> const x_post_ghost_vect = periodic_extrapolation_right(1, x_domain_vect);

Then we can create our 4 domains using the init_discretization function with the init_ghosted function that takes the vectors as parameters.

auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
        = ddc::init_discrete_space<DDimX>(
                DDimX::init_ghosted<DDimX>(x_domain_vect, x_pre_ghost_vect, x_post_ghost_vect));

To summarize, in this section we saw how to build a domain with non-uniform space discretization:

• Create an iterator for your domain.
• Call the init_discretization function with init_ghosted if you need ghost points or with init for a regular domain.

Dimension Y

For the Y dimension, we do the same. We start by defining our 3 vectors.

std::vector<double> const y_domain_vect = generate_random_vector(nb_y_points, y_start, y_end);
 
std::vector<double> const y_pre_ghost_vect = periodic_extrapolation_left(1, y_domain_vect);
std::vector<double> const y_post_ghost_vect = periodic_extrapolation_right(1, y_domain_vect);

And we use them to build the 4 domains in the same way as for 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>(y_domain_vect, y_pre_ghost_vect, y_post_ghost_vect));

Time dimension

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

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

The CFL conditions are more challenging to achieve for the case of non-uniform discretization. Since the spatial steps are not uniform, we first need to find the maximum of the inverse of the square of the spatial step for the X and Y dimensions. And then we obtain the minimum value for dt.

double const invdx2_max = ddc::transform_reduce(
        x_domain,
        0.,
        ddc::reducer::max<double>(),
        [](ddc::DiscreteElement<DDimX> ix) {
            return 1. / (ddc::distance_at_left(ix) * ddc::distance_at_right(ix));
        });
 
double const invdy2_max = ddc::transform_reduce(
        y_domain,
        0.,
        ddc::reducer::max<double>(),
        [](ddc::DiscreteElement<DDimY> iy) {
            return 1. / (ddc::distance_at_left(iy) * ddc::distance_at_right(iy));
        });
 
ddc::Coordinate<T> const max_dt(.5 / (kx * invdx2_max + ky * invdy2_max));

We can calculate the number of time steps and build the DiscreteDomain for the time dimension.

ddc::DiscreteVector<DDimT> const nb_time_steps(
        std::ceil((end_time - start_time) / max_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));

Time loop

Allocation and initialization are the same as for the uniform case. Let's focus on resolving the numerical scheme. The main difference in solving the numerical equation is that we need to account for the fact that the values of dx and dy on the left and right sides are different. We use the functions distance_at_left and distance_at_right to solve the equation.

ddc::parallel_for_each(
        next_temp.domain(),
        KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> const ixy) {
            ddc::DiscreteElement<DDimX> const ix(ixy);
            ddc::DiscreteElement<DDimY> const iy(ixy);
            double const dx_l = ddc::distance_at_left(ix);
            double const dx_r = ddc::distance_at_right(ix);
            double const dx_m = 0.5 * (dx_l + dx_r);
            double const dy_l = ddc::distance_at_left(iy);
            double const dy_r = ddc::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 * ddc::step<DDimT>()
                       * (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 * ddc::step<DDimT>()
                       * (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);
        });