DDC 0.0.0

a discrete domain computation library

examples/heat_equation.cpp

1// SPDX-License-Identifier: MIT
2
4#include <cmath>
5#include <iomanip>
6#include <iostream>
7#include <numeric>
8
9#include <ddc/ddc.hpp>
10
11#include <Kokkos_Core.hpp>
13
14
17struct X;
19
24
26// Our second continuous dimension
27struct Y;
28// Its uniform discretization
31
33// Our simulated time dimension
34struct T;
35// Its uniform discretization
38
39
41
45template <class ChunkType>
46void display(double time, ChunkType temp)
47{
48 double const mean_temp = ddc::transform_reduce(
49 temp.domain(),
50 0.,
52 temp)
53 / temp.domain().size();
54 std::cout << std::fixed << std::setprecision(3);
55 std::cout << "At t = " << time << ",\n";
56 std::cout << " * mean temperature = " << mean_temp << "\n";
57 // take a slice in the middle of the box
58 ddc::ChunkSpan temp_slice
59 = temp[ddc::get_domain<DDimY>(temp).front()
60 + ddc::get_domain<DDimY>(temp).size() / 2];
61 std::cout << " * temperature[y:"
62 << ddc::get_domain<DDimY>(temp).size() / 2 << "] = {";
65 ddc::get_domain<DDimX>(temp),
66 [=](ddc::DiscreteElement<DDimX> const ix) {
67 std::cout << std::setw(6) << temp_slice(ix);
68 });
69 std::cout << " }" << std::endl;
70}
72
73
75int main(int argc, char** argv)
76{
77 ddc::ScopeGuard scope(argc, argv);
78
79 // some parameters that would typically be read from some form of
80 // configuration file in a more realistic code
81
83 // Start of the domain of interest in the X dimension
84 double const x_start = -1.;
85 // End of the domain of interest in the X dimension
86 double const x_end = 1.;
87 // Number of discretization points in the X dimension
88 size_t const nb_x_points = 10;
89 // Thermal diffusion coefficient
90 double const kx = .01;
91 // Start of the domain of interest in the Y dimension
92 double const y_start = -1.;
93 // End of the domain of interest in the Y dimension
94 double const y_end = 1.;
95 // Number of discretization points in the Y dimension
96 size_t const nb_y_points = 100;
97 // Thermal diffusion coefficient
98 double const ky = .002;
99 // Simulated time at which to start simulation
100 double const start_time = 0.;
101 // Simulated time to reach as target of the simulation
102 double const end_time = 10.;
103 // Number of time-steps between outputs
104 ptrdiff_t const t_output_period = 10;
106
109 // Number of ghost points to use on each side in X
110 ddc::DiscreteVector<DDimX> static constexpr gwx {1};
112
114 // Initialization of the global domain in X with gwx ghost points on
115 // each side
116 auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
117 = ddc::init_discrete_space(DDimX::init_ghosted(
118 ddc::Coordinate<X>(x_start),
119 ddc::Coordinate<X>(x_end),
120 ddc::DiscreteVector<DDimX>(nb_x_points),
121 gwx));
123
125 // our zone at the start of the domain that will be mirrored to the
126 // ghost
128 x_domain_begin(x_domain.front(), x_post_ghost.extents());
129 // our zone at the end of the domain that will be mirrored to the
130 // ghost
131 ddc::DiscreteDomain const x_domain_end(
132 x_domain.back() - x_pre_ghost.extents() + 1,
133 x_pre_ghost.extents());
135
137 // Number of ghost points to use on each side in Y
138 ddc::DiscreteVector<DDimY> static constexpr gwy {1};
139
140 // Initialization of the global domain in Y with gwy ghost points on
141 // each side
142 auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
143 = ddc::init_discrete_space(DDimY::init_ghosted(
144 ddc::Coordinate<Y>(y_start),
145 ddc::Coordinate<Y>(y_end),
146 ddc::DiscreteVector<DDimY>(nb_y_points),
147 gwy));
148
149 // our zone at the start of the domain that will be mirrored to the
150 // ghost
152 y_domain_begin(y_domain.front(), y_post_ghost.extents());
153 // our zone at the end of the domain that will be mirrored to the
154 // ghost
155 ddc::DiscreteDomain const y_domain_end(
156 y_domain.back() - y_pre_ghost.extents() + 1,
157 y_pre_ghost.extents());
159
161 // max(1/dx^2)
162 double const invdx2_max = ddc::transform_reduce(
163 x_domain,
164 0.,
167 return 1.
170 });
171 // max(1/dy^2)
172 double const invdy2_max = ddc::transform_reduce(
173 y_domain,
174 0.,
177 return 1.
180 });
181 ddc::Coordinate<T> const max_dt {
182 .5 / (kx * invdx2_max + ky * invdy2_max)};
183
184 // number of time intervals required to reach the end time
185 ddc::DiscreteVector<DDimT> const nb_time_steps {
186 std::ceil((end_time - start_time) / max_dt) + .2};
187 // Initialization of the global domain in time:
188 // - the number of discrete time-points is equal to the number of
189 // steps + 1
190 ddc::DiscreteDomain<DDimT> const time_domain
192 DDimT::
193 init(ddc::Coordinate<T>(start_time),
194 ddc::Coordinate<T>(end_time),
195 nb_time_steps + 1));
197
199 // Maps temperature into the full domain (including ghosts) twice:
200 // - once for the last fully computed time-step
201 ddc::Chunk ghosted_last_temp(
203 DDimX,
204 DDimY>(ghosted_x_domain, ghosted_y_domain),
206
207 // - once for time-step being computed
208 ddc::Chunk ghosted_next_temp(
210 DDimX,
211 DDimY>(ghosted_x_domain, ghosted_y_domain),
214
216 ddc::ChunkSpan const ghosted_initial_temp
217 = ghosted_last_temp.span_view();
218 // Initialize the temperature on the main domain
221 ddc::DiscreteDomain<DDimX, DDimY>(x_domain, y_domain),
222 DDC_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> const ixy) {
223 double const x
224 = ddc::coordinate(ddc::select<DDimX>(ixy));
225 double const y
226 = ddc::coordinate(ddc::select<DDimY>(ixy));
227 ghosted_initial_temp(ixy)
228 = 9.999 * ((x * x + y * y) < 0.25);
229 });
231
232 ddc::Chunk ghosted_temp(
234 DDimX,
235 DDimY>(ghosted_x_domain, ghosted_y_domain),
237
238
240 // display the initial data
241 ddc::deepcopy(ghosted_temp, ghosted_last_temp);
242 display(ddc::coordinate(time_domain.front()),
243 ghosted_temp[x_domain][y_domain]);
244 // time of the iteration where the last output happened
245 ddc::DiscreteElement<DDimT> last_output = time_domain.front();
247
249 for (auto const iter :
252
254 // Periodic boundary conditions
256 ghosted_last_temp[x_pre_ghost][y_domain],
257 ghosted_last_temp[y_domain][x_domain_end]);
259 ghosted_last_temp[y_domain][x_post_ghost],
260 ghosted_last_temp[y_domain][x_domain_begin]);
262 ghosted_last_temp[x_domain][y_pre_ghost],
263 ghosted_last_temp[x_domain][y_domain_end]);
265 ghosted_last_temp[x_domain][y_post_ghost],
266 ghosted_last_temp[x_domain][y_domain_begin]);
268
270 // a span excluding ghosts of the temperature at the time-step we
271 // will build
272 ddc::ChunkSpan const next_temp {
273 ghosted_next_temp[x_domain][y_domain]};
274 // a read-only view of the temperature at the previous time-step
275 ddc::ChunkSpan const last_temp {ghosted_last_temp.span_view()};
277
279 // Stencil computation on the main domain
282 next_temp.domain(),
283 DDC_LAMBDA(
286 = ddc::select<DDimX>(ixy);
288 = ddc::select<DDimY>(ixy);
289 double const dx_l = ddc::distance_at_left(ix);
290 double const dx_r = ddc::distance_at_right(ix);
291 double const dx_m = 0.5 * (dx_l + dx_r);
292 double const dy_l = ddc::distance_at_left(iy);
293 double const dy_r = ddc::distance_at_right(iy);
294 double const dy_m = 0.5 * (dy_l + dy_r);
295 next_temp(ix, iy) = last_temp(ix, iy);
296 next_temp(ix, iy)
297 += kx * ddc::step<DDimT>()
298 * (dx_l * last_temp(ix + 1, iy)
299 - 2.0 * dx_m * last_temp(ix, iy)
300 + dx_r * last_temp(ix - 1, iy))
301 / (dx_l * dx_m * dx_r);
302 next_temp(ix, iy)
303 += ky * ddc::step<DDimT>()
304 * (dy_l * last_temp(ix, iy + 1)
305 - 2.0 * dy_m * last_temp(ix, iy)
306 + dy_r * last_temp(ix, iy - 1))
307 / (dy_l * dy_m * dy_r);
308 });
310
312 if (iter - last_output >= t_output_period) {
313 last_output = iter;
314 ddc::deepcopy(ghosted_temp, ghosted_last_temp);
315 display(ddc::coordinate(iter),
316 ghosted_temp[x_domain][y_domain]);
317 }
319
321 // Swap our two buffers
322 std::swap(ghosted_last_temp, ghosted_next_temp);
324 }
325
327 if (last_output < time_domain.back()) {
328 ddc::deepcopy(ghosted_temp, ghosted_last_temp);
329 display(ddc::coordinate(time_domain.back()),
330 ghosted_temp[x_domain][y_domain]);
331 }
333}
Definition: discrete_domain.hpp:24
constexpr discrete_element_type back() const noexcept
Definition: discrete_domain.hpp:120
constexpr discrete_element_type front() const noexcept
Definition: discrete_domain.hpp:115
constexpr DiscreteDomain remove_first(mlength_type n) const
Definition: discrete_domain.hpp:135
A DiscreteElement identifies an element of the discrete dimension.
Definition: discrete_element.hpp:125
A DiscreteVector is a vector in the discrete dimension.
Definition: discrete_vector.hpp:228
Definition: kokkos_allocator.hpp:14
Definition: scope_guard.hpp:14
UniformPointSampling models a uniform discretization of the provided continuous dimension.
Definition: uniform_point_sampling.hpp:24
constexpr parallel_device_policy parallel_device
Definition: for_each.hpp:198
constexpr serial_host_policy serial_host
Definition: for_each.hpp:196
DDC_INLINE_FUNCTION Coordinate< CDim > distance_at_left(DiscreteElement< NonUniformPointSampling< CDim > > i)
Definition: non_uniform_point_sampling.hpp:144
DDC_INLINE_FUNCTION Coordinate< typename DDim::continuous_dimension_type... > coordinate(DiscreteElement< DDim... > const &c)
Definition: coordinate_md.hpp:12
auto deepcopy(ChunkDst &&dst, ChunkSrc &&src)
Copy the content of a borrowed chunk into another.
Definition: deepcopy.hpp:19
void init_discrete_space(Args &&... args)
Initialize (emplace) a global singleton discrete space.
Definition: discrete_space.hpp:140
T transform_reduce(serial_host_policy policy, DiscreteDomain< DDims... > const &domain, T neutral, BinaryReductionOp &&reduce, UnaryTransformOp &&transform) noexcept
A reduction over a nD domain using the Serial execution policy.
Definition: transform_reduce.hpp:267
DDC_INLINE_FUNCTION Coordinate< CDim > distance_at_right(DiscreteElement< NonUniformPointSampling< CDim > > i)
Definition: non_uniform_point_sampling.hpp:151
void for_each(serial_host_policy, DiscreteDomain< DDims... > const &domain, Functor &&f) noexcept
iterates over a nD domain using the serial execution policy
Definition: for_each.hpp:127
detail::TaggedVector< CoordinateElement, CDims... > Coordinate
A Coordinate represents a coordinate in the continuous space.
Definition: coordinate.hpp:21
Definition: chunk.hpp:15
Definition: chunk_span.hpp:26
Definition: reducer.hpp:96
Definition: reducer.hpp:10