diff --git a/first_steps.html b/first_steps.html index 1a40a2806..2a6289d1f 100644 --- a/first_steps.html +++ b/first_steps.html @@ -54,9 +54,7 @@

a discrete #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>.

+

As you can see, to use DDC, we have to include <ddc/ddc.hpp>

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.;
diff --git a/heat_equation.html b/heat_equation.html index d5806ec35..7fb55dee5 100644 --- a/heat_equation.html +++ b/heat_equation.html @@ -54,280 +54,278 @@

a discrete 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
+
11
+
12
+
15struct X;
+
17
+ +
22
+
24// Our second continuous dimension
+
25struct Y;
+
26// Its uniform discretization
+ +
29
+
31// Our simulated time dimension
+
32struct T;
+
33// Its uniform discretization
+ +
36
+
37
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}
+
43template <class ChunkType>
+
44void display(double time, ChunkType temp)
+
45{
+
46 double const mean_temp = ddc::transform_reduce(
+
47 temp.domain(),
+
48 0.,
+ +
50 temp)
+
51 / temp.domain().size();
+
52 std::cout << std::fixed << std::setprecision(3);
+
53 std::cout << "At t = " << time << ",\n";
+
54 std::cout << " * mean temperature = " << mean_temp << "\n";
+
55 // take a slice in the middle of the box
+
56 ddc::ChunkSpan temp_slice
+
57 = temp[ddc::get_domain<DDimY>(temp).front()
+
58 + ddc::get_domain<DDimY>(temp).size() / 2];
+
59 std::cout << " * temperature[y:"
+
60 << ddc::get_domain<DDimY>(temp).size() / 2 << "] = {";
+ + +
63 ddc::get_domain<DDimX>(temp),
+
64 [=](ddc::DiscreteElement<DDimX> const ix) {
+
65 std::cout << std::setw(6) << temp_slice(ix);
+
66 });
+
67 std::cout << " }" << std::endl;
+
68}
+
70
+
71
+
73int main(int argc, char** argv)
+
74{
+
75 ddc::ScopeGuard scope(argc, argv);
+
76
+
77 // some parameters that would typically be read from some form of
+
78 // configuration file in a more realistic code
+
79
+
81 // Start of the domain of interest in the X dimension
+
82 double const x_start = -1.;
+
83 // End of the domain of interest in the X dimension
+
84 double const x_end = 1.;
+
85 // Number of discretization points in the X dimension
+
86 size_t const nb_x_points = 10;
+
87 // Thermal diffusion coefficient
+
88 double const kx = .01;
+
89 // Start of the domain of interest in the Y dimension
+
90 double const y_start = -1.;
+
91 // End of the domain of interest in the Y dimension
+
92 double const y_end = 1.;
+
93 // Number of discretization points in the Y dimension
+
94 size_t const nb_y_points = 100;
+
95 // Thermal diffusion coefficient
+
96 double const ky = .002;
+
97 // Simulated time at which to start simulation
+
98 double const start_time = 0.;
+
99 // Simulated time to reach as target of the simulation
+
100 double const end_time = 10.;
+
101 // Number of time-steps between outputs
+
102 ptrdiff_t const t_output_period = 10;
+
104
+
107 // Number of ghost points to use on each side in X
+
108 ddc::DiscreteVector<DDimX> static constexpr gwx {1};
+
110
+
112 // Initialization of the global domain in X with gwx ghost points on
+
113 // each side
+
114 auto const [x_domain, ghosted_x_domain, x_pre_ghost, x_post_ghost]
+
115 = ddc::init_discrete_space(DDimX::init_ghosted(
+
116 ddc::Coordinate<X>(x_start),
+
117 ddc::Coordinate<X>(x_end),
+
118 ddc::DiscreteVector<DDimX>(nb_x_points),
+
119 gwx));
+
121
+
123 // our zone at the start of the domain that will be mirrored to the
+
124 // ghost
+ +
126 x_domain_begin(x_domain.front(), x_post_ghost.extents());
+
127 // our zone at the end of the domain that will be mirrored to the
+
128 // ghost
+
129 ddc::DiscreteDomain const x_domain_end(
+
130 x_domain.back() - x_pre_ghost.extents() + 1,
+
131 x_pre_ghost.extents());
+
133
+
135 // Number of ghost points to use on each side in Y
+
136 ddc::DiscreteVector<DDimY> static constexpr gwy {1};
+
137
+
138 // Initialization of the global domain in Y with gwy ghost points on
+
139 // each side
+
140 auto const [y_domain, ghosted_y_domain, y_pre_ghost, y_post_ghost]
+
141 = ddc::init_discrete_space(DDimY::init_ghosted(
+
142 ddc::Coordinate<Y>(y_start),
+
143 ddc::Coordinate<Y>(y_end),
+
144 ddc::DiscreteVector<DDimY>(nb_y_points),
+
145 gwy));
+
146
+
147 // our zone at the start of the domain that will be mirrored to the
+
148 // ghost
+ +
150 y_domain_begin(y_domain.front(), y_post_ghost.extents());
+
151 // our zone at the end of the domain that will be mirrored to the
+
152 // ghost
+
153 ddc::DiscreteDomain const y_domain_end(
+
154 y_domain.back() - y_pre_ghost.extents() + 1,
+
155 y_pre_ghost.extents());
+
157
+
159 // max(1/dx^2)
+
160 double const invdx2_max = ddc::transform_reduce(
+
161 x_domain,
+
162 0.,
+ + +
165 return 1.
+ + +
168 });
+
169 // max(1/dy^2)
+
170 double const invdy2_max = ddc::transform_reduce(
+
171 y_domain,
+
172 0.,
+ + +
175 return 1.
+ + +
178 });
+
179 ddc::Coordinate<T> const max_dt {
+
180 .5 / (kx * invdx2_max + ky * invdy2_max)};
+
181
+
182 // number of time intervals required to reach the end time
+
183 ddc::DiscreteVector<DDimT> const nb_time_steps {
+
184 std::ceil((end_time - start_time) / max_dt) + .2};
+
185 // Initialization of the global domain in time:
+
186 // - the number of discrete time-points is equal to the number of
+
187 // steps + 1
+
188 ddc::DiscreteDomain<DDimT> const time_domain
+ +
190 DDimT::
+
191 init(ddc::Coordinate<T>(start_time),
+
192 ddc::Coordinate<T>(end_time),
+
193 nb_time_steps + 1));
+
195
+
197 // Maps temperature into the full domain (including ghosts) twice:
+
198 // - once for the last fully computed time-step
+
199 ddc::Chunk ghosted_last_temp(
+ +
201 DDimX,
+
202 DDimY>(ghosted_x_domain, ghosted_y_domain),
+ +
204
+
205 // - once for time-step being computed
+
206 ddc::Chunk ghosted_next_temp(
+ +
208 DDimX,
+
209 DDimY>(ghosted_x_domain, ghosted_y_domain),
+ +
212
+
214 ddc::ChunkSpan const ghosted_initial_temp
+
215 = ghosted_last_temp.span_view();
+
216 // Initialize the temperature on the main domain
+ + +
219 ddc::DiscreteDomain<DDimX, DDimY>(x_domain, y_domain),
+
220 DDC_LAMBDA(ddc::DiscreteElement<DDimX, DDimY> const ixy) {
+
221 double const x
+
222 = ddc::coordinate(ddc::select<DDimX>(ixy));
+
223 double const y
+
224 = ddc::coordinate(ddc::select<DDimY>(ixy));
+
225 ghosted_initial_temp(ixy)
+
226 = 9.999 * ((x * x + y * y) < 0.25);
+
227 });
+
229
+
230 ddc::Chunk ghosted_temp(
+ +
232 DDimX,
+
233 DDimY>(ghosted_x_domain, ghosted_y_domain),
+ +
235
+
236
+
238 // display the initial data
+
239 ddc::deepcopy(ghosted_temp, ghosted_last_temp);
+
240 display(ddc::coordinate(time_domain.front()),
+
241 ghosted_temp[x_domain][y_domain]);
+
242 // time of the iteration where the last output happened
+
243 ddc::DiscreteElement<DDimT> last_output = time_domain.front();
+
245
+
247 for (auto const iter :
+ +
250
+
252 // Periodic boundary conditions
+ +
254 ghosted_last_temp[x_pre_ghost][y_domain],
+
255 ghosted_last_temp[y_domain][x_domain_end]);
+ +
257 ghosted_last_temp[y_domain][x_post_ghost],
+
258 ghosted_last_temp[y_domain][x_domain_begin]);
+ +
260 ghosted_last_temp[x_domain][y_pre_ghost],
+
261 ghosted_last_temp[x_domain][y_domain_end]);
+ +
263 ghosted_last_temp[x_domain][y_post_ghost],
+
264 ghosted_last_temp[x_domain][y_domain_begin]);
+
266
+
268 // a span excluding ghosts of the temperature at the time-step we
+
269 // will build
+
270 ddc::ChunkSpan const next_temp {
+
271 ghosted_next_temp[x_domain][y_domain]};
+
272 // a read-only view of the temperature at the previous time-step
+
273 ddc::ChunkSpan const last_temp {ghosted_last_temp.span_view()};
+
275
+
277 // Stencil computation on the main domain
+ + +
280 next_temp.domain(),
+
281 DDC_LAMBDA(
+ + +
284 = ddc::select<DDimX>(ixy);
+ +
286 = ddc::select<DDimY>(ixy);
+
287 double const dx_l = ddc::distance_at_left(ix);
+
288 double const dx_r = ddc::distance_at_right(ix);
+
289 double const dx_m = 0.5 * (dx_l + dx_r);
+
290 double const dy_l = ddc::distance_at_left(iy);
+
291 double const dy_r = ddc::distance_at_right(iy);
+
292 double const dy_m = 0.5 * (dy_l + dy_r);
+
293 next_temp(ix, iy) = last_temp(ix, iy);
+
294 next_temp(ix, iy)
+
295 += kx * ddc::step<DDimT>()
+
296 * (dx_l * last_temp(ix + 1, iy)
+
297 - 2.0 * dx_m * last_temp(ix, iy)
+
298 + dx_r * last_temp(ix - 1, iy))
+
299 / (dx_l * dx_m * dx_r);
+
300 next_temp(ix, iy)
+
301 += ky * ddc::step<DDimT>()
+
302 * (dy_l * last_temp(ix, iy + 1)
+
303 - 2.0 * dy_m * last_temp(ix, iy)
+
304 + dy_r * last_temp(ix, iy - 1))
+
305 / (dy_l * dy_m * dy_r);
+
306 });
+
308
+
310 if (iter - last_output >= t_output_period) {
+
311 last_output = iter;
+
312 ddc::deepcopy(ghosted_temp, ghosted_last_temp);
+
313 display(ddc::coordinate(iter),
+
314 ghosted_temp[x_domain][y_domain]);
+
315 }
+
317
+
319 // Swap our two buffers
+
320 std::swap(ghosted_last_temp, ghosted_next_temp);
+
322 }
+
323
+
325 if (last_output < time_domain.back()) {
+
326 ddc::deepcopy(ghosted_temp, ghosted_last_temp);
+
327 display(ddc::coordinate(time_domain.back()),
+
328 ghosted_temp[x_domain][y_domain]);
+
329 }
+
331}
Definition: discrete_domain.hpp:24
KOKKOS_FUNCTION constexpr discrete_element_type front() const noexcept
Definition: discrete_domain.hpp:117
KOKKOS_FUNCTION constexpr DiscreteDomain remove_first(mlength_type n) const
Definition: discrete_domain.hpp:137