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.
As you can see, to use DDC, we have to include <ddc/ddc.hpp>
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.
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.
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.
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.
We do the same thing for the second dimension Y.
And once again, now for the time dimension.
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).
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.
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.
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 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.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
Then we initialize the domain along this dimension just like we did with the X 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.
Then we use the CFL condition to determine the time step of the simulation.
Finally, we determine the number of time steps and as we did with the X and Y dimensions, we create the time domain.
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:
Note that the DeviceAllocator is responsible for allocating memory on the default memory space.
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.
Then, we iterate over each DiscreteElement of the domain to fill ghosted_initial_temp with the initial values of the simulation.
To display the data, a chunk is created on the host.
We deepcopy the data from the ghosted_last_temp chunk to ghosted_temp on the host.
And we display the initial data.
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.We then solve the equation.