DDC 0.1.0
Loading...
Searching...
No Matches
examples/uniform_heat_equation.cpp
1// Copyright (C) The DDC development team, see COPYRIGHT.md file
2//
3// SPDX-License-Identifier: MIT
4
6#include <cmath>
7#include <cstddef>
8#include <iomanip>
9#include <iostream>
10#include <string>
11#include <utility>
12
13#include <ddc/ddc.hpp>
14
15#include <Kokkos_Core.hpp>
16#if defined(DDC_EXAMPLE_WITH_PDI)
17#include <ddc/pdi.hpp>
18
19#include <paraconf.h>
20#include <pdi.h>
21#endif
23
25struct X;
27
29struct DDimX : ddc::UniformPointSampling<X>
30{
31};
33
35struct Y;
36struct DDimY : ddc::UniformPointSampling<Y>
37{
38};
40
42struct T;
43struct DDimT : ddc::UniformPointSampling<T>
44{
45};
47
55template <class ChunkType>
56void display(double time, ChunkType temp)
57{
58 double const mean_temp = ddc::transform_reduce(
59 temp.domain(),
60 0.,
61 ddc::reducer::sum<double>(),
62 temp)
63 / temp.domain().size();
64 std::cout << std::fixed << std::setprecision(3);
65 std::cout << "At t = " << time << ",\n";
66 std::cout << " * mean temperature = " << mean_temp << "\n";
67 ddc::ChunkSpan const temp_slice
68 = temp[ddc::get_domain<DDimY>(temp).front()
69 + ddc::get_domain<DDimY>(temp).size() / 2];
70 std::cout << " * temperature[y:"
71 << ddc::get_domain<DDimY>(temp).size() / 2 << "] = {";
72 ddc::for_each(
73 ddc::get_domain<DDimX>(temp),
74 [=](ddc::DiscreteElement<DDimX> const ix) {
75 std::cout << std::setw(6) << temp_slice(ix);
76 });
77 std::cout << " }\n" << std::flush;
78
79#if defined(DDC_EXAMPLE_WITH_PDI)
80 ddc::PdiEvent("display")
81 .with("temp", temp)
82 .and_with("mean_temp", mean_temp)
83 .and_with("temp_slice", temp_slice);
84#endif
85}
87
88int main(int argc, char** argv)
89{
90#if defined(DDC_EXAMPLE_WITH_PDI)
91 PC_tree_t pdi_conf = PC_parse_string("");
92 PDI_init(pdi_conf);
93#endif
94 Kokkos::ScopeGuard const kokkos_scope(argc, argv);
95 ddc::ScopeGuard const ddc_scope(argc, argv);
96
99 double const x_start = -1.;
100 double const x_end = 1.;
101 std::size_t const nb_x_points = 10;
102 double const kx = .01;
105 double const y_start = -1.;
106 double const y_end = 1.;
107 std::size_t const nb_y_points = 100;
108 double const ky = .002;
111 double const start_time = 0.;
112 double const end_time = 10.;
114 std::ptrdiff_t const t_output_period = 10;
116
118 ddc::DiscreteVector<DDimX> const gwx(1);
120
122 auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
123 = ddc::init_discrete_space<DDimX>(DDimX::init_ghosted<DDimX>(
124 ddc::Coordinate<X>(x_start),
125 ddc::Coordinate<X>(x_end),
126 ddc::DiscreteVector<DDimX>(nb_x_points),
127 gwx));
129
131 ddc::DiscreteDomain<DDimX> const
132 x_domain_begin(x_domain.front(), x_post_ghost.extents());
133 ddc::DiscreteDomain<DDimX> const x_domain_end(
134 x_domain.back() - x_pre_ghost.extents() + 1,
135 x_pre_ghost.extents());
137
139 ddc::DiscreteVector<DDimY> const gwy(1);
140
141 auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
142 = ddc::init_discrete_space<DDimY>(DDimY::init_ghosted<DDimY>(
143 ddc::Coordinate<Y>(y_start),
144 ddc::Coordinate<Y>(y_end),
145 ddc::DiscreteVector<DDimY>(nb_y_points),
146 gwy));
147
148 ddc::DiscreteDomain<DDimY> const
149 y_domain_begin(y_domain.front(), y_post_ghost.extents());
150
151 ddc::DiscreteDomain<DDimY> const y_domain_end(
152 y_domain.back() - y_pre_ghost.extents() + 1,
153 y_pre_ghost.extents());
155
157 double const dx = ddc::step<DDimX>();
158 double const dy = ddc::step<DDimY>();
159 double const invdx2 = 1. / (dx * dx);
160 double const invdy2 = 1. / (dy * dy);
161
162 ddc::Coordinate<T> const dt(.5 / (kx * invdx2 + ky * invdy2));
164
166 ddc::DiscreteVector<DDimT> const nb_time_steps(
167 std::ceil((end_time - start_time) / dt) + .2);
168
169 ddc::DiscreteDomain<DDimT> const time_domain
170 = ddc::init_discrete_space<DDimT>(DDimT::init<DDimT>(
171 ddc::Coordinate<T>(start_time),
172 ddc::Coordinate<T>(end_time),
173 nb_time_steps + 1));
175
177 ddc::Chunk ghosted_last_temp(
178 "ghosted_last_temp",
179 ddc::DiscreteDomain<
180 DDimX,
181 DDimY>(ghosted_x_domain, ghosted_y_domain),
182 ddc::DeviceAllocator<double>());
183
184 ddc::Chunk ghosted_next_temp(
185 "ghosted_next_temp",
186 ddc::DiscreteDomain<
187 DDimX,
188 DDimY>(ghosted_x_domain, ghosted_y_domain),
189 ddc::DeviceAllocator<double>());
191
193 ddc::ChunkSpan const ghosted_initial_temp
194 = ghosted_last_temp.span_view();
196
198 ddc::parallel_for_each(
199 ddc::DiscreteDomain<DDimX, DDimY>(x_domain, y_domain),
200 KOKKOS_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> const ixy) {
201 double const x = ddc::coordinate(
202 ddc::DiscreteElement<DDimX>(ixy));
203 double const y = ddc::coordinate(
204 ddc::DiscreteElement<DDimY>(ixy));
205 ghosted_initial_temp(ixy)
206 = 9.999 * ((x * x + y * y) < 0.25);
207 });
209
211 ddc::Chunk ghosted_temp
212 = ddc::create_mirror(ghosted_last_temp.span_cview());
214
216 ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp);
218
220 display(ddc::coordinate(time_domain.front()),
221 ghosted_temp[x_domain][y_domain]);
223
225 ddc::DiscreteElement<DDimT> last_output_iter = time_domain.front();
227
229 for (ddc::DiscreteElement<DDimT> const iter :
230 time_domain.remove_first(ddc::DiscreteVector<DDimT>(1))) {
232
234 ddc::parallel_deepcopy(
235 ghosted_last_temp[ddc::DiscreteDomain<
236 DDimX,
237 DDimY>(x_pre_ghost, y_domain)],
238 ghosted_last_temp[ddc::DiscreteDomain<
239 DDimX,
240 DDimY>(y_domain, x_domain_end)]);
241 ddc::parallel_deepcopy(
242 ghosted_last_temp[ddc::DiscreteDomain<
243 DDimX,
244 DDimY>(y_domain, x_post_ghost)],
245 ghosted_last_temp[ddc::DiscreteDomain<
246 DDimX,
247 DDimY>(y_domain, x_domain_begin)]);
248 ddc::parallel_deepcopy(
249 ghosted_last_temp[ddc::DiscreteDomain<
250 DDimX,
251 DDimY>(x_domain, y_pre_ghost)],
252 ghosted_last_temp[ddc::DiscreteDomain<
253 DDimX,
254 DDimY>(x_domain, y_domain_end)]);
255 ddc::parallel_deepcopy(
256 ghosted_last_temp[ddc::DiscreteDomain<
257 DDimX,
258 DDimY>(x_domain, y_post_ghost)],
259 ghosted_last_temp[ddc::DiscreteDomain<
260 DDimX,
261 DDimY>(x_domain, y_domain_begin)]);
263
265 ddc::ChunkSpan const next_temp(
266 ghosted_next_temp[ddc::DiscreteDomain<
267 DDimX,
268 DDimY>(x_domain, y_domain)]);
269 ddc::ChunkSpan const last_temp(ghosted_last_temp.span_cview());
271
273 ddc::parallel_for_each(
274 next_temp.domain(),
275 KOKKOS_LAMBDA(
276 ddc::DiscreteElement<DDimX, DDimY> const ixy) {
277 ddc::DiscreteElement<DDimX> const ix(ixy);
278 ddc::DiscreteElement<DDimY> const iy(ixy);
279 double const dt = ddc::step<DDimT>();
280
281 next_temp(ix, iy) = last_temp(ix, iy);
282 next_temp(ix, iy) += kx * dt
283 * (last_temp(ix + 1, iy)
284 - 2.0 * last_temp(ix, iy)
285 + last_temp(ix - 1, iy))
286 * invdx2;
287
288 next_temp(ix, iy) += ky * dt
289 * (last_temp(ix, iy + 1)
290 - 2.0 * last_temp(ix, iy)
291 + last_temp(ix, iy - 1))
292 * invdy2;
293 });
295
297 if (iter - last_output_iter >= t_output_period) {
298 last_output_iter = iter;
299 ddc::parallel_deepcopy(ghosted_temp, ghosted_next_temp);
300 display(ddc::coordinate(iter),
301 ghosted_temp[ddc::DiscreteDomain<
302 DDimX,
303 DDimY>(x_domain, y_domain)]);
304 }
306
308 std::swap(ghosted_last_temp, ghosted_next_temp);
310 }
311
313 if (last_output_iter < time_domain.back()) {
314 ddc::parallel_deepcopy(ghosted_temp, ghosted_last_temp);
315 display(ddc::coordinate(time_domain.back()),
316 ghosted_temp[ddc::DiscreteDomain<
317 DDimX,
318 DDimY>(x_domain, y_domain)]);
319 }
321
322#if defined(DDC_EXAMPLE_WITH_PDI)
323 PDI_finalize();
324 PC_tree_destroy(&pdi_conf);
325#endif
326}
ddc::DiscreteDomain::front
KOKKOS_FUNCTION constexpr discrete_element_type front() const noexcept
Definition discrete_domain.hpp:140
ddc::DiscreteDomain::back
KOKKOS_FUNCTION constexpr discrete_element_type back() const noexcept
Definition discrete_domain.hpp:145
ddc::DiscreteElement
A DiscreteElement identifies an element of the discrete dimension.
Definition discrete_element.hpp:149
ddc::DiscreteVector
A DiscreteVector is a vector in the discrete dimension.
Definition discrete_vector.hpp:257
ddc::PdiEvent::and_with
PdiEvent & and_with(std::string const &name, T &&t)
Definition pdi.hpp:109
ddc::PdiEvent::with
PdiEvent & with(std::string const &name, BorrowedChunk &&data)
Definition pdi.hpp:64
ddc::UniformPointSampling
UniformPointSampling models a uniform discretization of the provided continuous dimension.
Definition uniform_point_sampling.hpp:36
ddc
The top-level namespace of DDC.
Definition aligned_allocator.hpp:11
ddc::parallel_deepcopy
auto parallel_deepcopy(ChunkDst &&dst, ChunkSrc &&src)
Copy the content of a borrowed chunk into another.
Definition parallel_deepcopy.hpp:22
ddc::for_each
void for_each(DiscreteDomain< DDims... > const &domain, Functor &&f) noexcept
iterates over a nD domain in serial
Definition for_each.hpp:43
ddc::coordinate
KOKKOS_FUNCTION Coordinate< typename DDim::continuous_dimension_type... > coordinate(DiscreteElement< DDim... > const &c)
Definition coordinate.hpp:32
ddc::init_discrete_space
void init_discrete_space(Args &&... args)
Initialize (emplace) a global singleton discrete space.
Definition discrete_space.hpp:170
ddc::transform_reduce
T transform_reduce(DiscreteDomain< DDims... > const &domain, T neutral, BinaryReductionOp &&reduce, UnaryTransformOp &&transform) noexcept
A reduction over a nD domain in serial.
Definition transform_reduce.hpp:66
ddc::parallel_for_each
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:83
ddc::Coordinate
detail::TaggedVector< CoordinateElement, CDims... > Coordinate
A Coordinate represents a coordinate in the continuous space.
Definition coordinate.hpp:29
ddc::create_mirror
auto create_mirror(Space const &space, ChunkSpan< ElementType, Support, Layout, MemorySpace > const &src)
Definition create_mirror.hpp:20