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 
17 struct X;
19 
22 using DDimX = ddc::UniformPointSampling<X>;
24 
26 // Our second continuous dimension
27 struct Y;
28 // Its uniform discretization
29 using DDimY = ddc::UniformPointSampling<Y>;
31 
33 // Our simulated time dimension
34 struct T;
35 // Its uniform discretization
36 using DDimT = ddc::UniformPointSampling<T>;
38 
39 
41 
45 template <class ChunkType>
46 void display(double time, ChunkType temp)
47 {
48  double const mean_temp = 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 
75 int 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  size_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
127  ddc::DiscreteDomain const
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
151  ddc::DiscreteDomain const
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 = transform_reduce(
163  x_domain,
164  0.,
167  return 1.
168  / (distance_at_left(ix) * distance_at_right(ix));
169  });
170  // max(1/dy^2)
171  double const invdy2_max = transform_reduce(
172  y_domain,
173  0.,
176  return 1.
177  / (distance_at_left(iy) * distance_at_right(iy));
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 = coordinate(ddc::select<DDimX>(ixy));
222  double const y = coordinate(ddc::select<DDimY>(ixy));
223  ghosted_initial_temp(ixy)
224  = 9.999 * ((x * x + y * y) < 0.25);
225  });
227 
228  ddc::Chunk ghosted_temp(
230  DDimX,
231  DDimY>(ghosted_x_domain, ghosted_y_domain),
233 
234 
236  // display the initial data
237  ddc::deepcopy(ghosted_temp, ghosted_last_temp);
238  display(coordinate(time_domain.front()),
239  ghosted_temp[x_domain][y_domain]);
240  // time of the iteration where the last output happened
241  ddc::DiscreteElement<DDimT> last_output = time_domain.front();
243 
245  for (auto const iter :
246  time_domain.remove_first(ddc::DiscreteVector<DDimT>(1))) {
248 
250  // Periodic boundary conditions
252  ghosted_last_temp[x_pre_ghost][y_domain],
253  ghosted_last_temp[y_domain][x_domain_end]);
255  ghosted_last_temp[y_domain][x_post_ghost],
256  ghosted_last_temp[y_domain][x_domain_begin]);
258  ghosted_last_temp[x_domain][y_pre_ghost],
259  ghosted_last_temp[x_domain][y_domain_end]);
261  ghosted_last_temp[x_domain][y_post_ghost],
262  ghosted_last_temp[x_domain][y_domain_begin]);
264 
266  // a span excluding ghosts of the temperature at the time-step we
267  // will build
268  ddc::ChunkSpan const next_temp {
269  ghosted_next_temp[x_domain][y_domain]};
270  // a read-only view of the temperature at the previous time-step
271  ddc::ChunkSpan const last_temp {ghosted_last_temp.span_view()};
273 
275  // Stencil computation on the main domain
278  next_temp.domain(),
279  DDC_LAMBDA(
282  = ddc::select<DDimX>(ixy);
284  = ddc::select<DDimY>(ixy);
285  double const dx_l = distance_at_left(ix);
286  double const dx_r = distance_at_right(ix);
287  double const dx_m = 0.5 * (dx_l + dx_r);
288  double const dy_l = distance_at_left(iy);
289  double const dy_r = distance_at_right(iy);
290  double const dy_m = 0.5 * (dy_l + dy_r);
291  next_temp(ix, iy) = last_temp(ix, iy);
292  next_temp(ix, iy)
293  += kx * max_dt
294  * (dx_l * last_temp(ix + 1, iy)
295  - 2.0 * dx_m * last_temp(ix, iy)
296  + dx_r * last_temp(ix - 1, iy))
297  / (dx_l * dx_m * dx_r);
298  next_temp(ix, iy)
299  += ky * max_dt
300  * (dy_l * last_temp(ix, iy + 1)
301  - 2.0 * dy_m * last_temp(ix, iy)
302  + dy_r * last_temp(ix, iy - 1))
303  / (dy_l * dy_m * dy_r);
304  });
306 
308  if (iter - last_output >= t_output_period) {
309  last_output = iter;
310  ddc::deepcopy(ghosted_temp, ghosted_last_temp);
311  display(coordinate(iter), ghosted_temp[x_domain][y_domain]);
312  }
314 
316  // Swap our two buffers
317  std::swap(ghosted_last_temp, ghosted_next_temp);
319  }
320 
322  if (last_output < time_domain.back()) {
323  ddc::deepcopy(ghosted_temp, ghosted_last_temp);
324  display(coordinate(time_domain.back()),
325  ghosted_temp[x_domain][y_domain]);
326  }
328 }
Definition: discrete_domain.hpp:24
constexpr discrete_element_type back() const noexcept
Definition: discrete_domain.hpp:132
constexpr discrete_element_type front() const noexcept
Definition: discrete_domain.hpp:127
constexpr DiscreteDomain remove_first(mlength_type n) const
Definition: discrete_domain.hpp:147
A DiscreteElement identifies an element of the discrete dimension.
Definition: discrete_element.hpp:124
A DiscreteVector is a vector in the discrete dimension.
Definition: discrete_vector.hpp:208
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:23
constexpr parallel_device_policy parallel_device
Definition: for_each.hpp:173
constexpr serial_host_policy serial_host
Definition: for_each.hpp:171
T transform_reduce([[maybe_unused]] 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:231
ddc::Coordinate< QueryDDims... > coordinate(DiscreteDomain< DDims... > const &domain, DiscreteElement< QueryDDims... > const &icoord) noexcept
Definition: discrete_domain.hpp:261
DDC_INLINE_FUNCTION Coordinate< CDim > distance_at_left(DiscreteElement< NonUniformPointSampling< CDim >> i)
Definition: non_uniform_point_sampling.hpp:144
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:109
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:101
ddc_detail::TaggedVector< CoordinateElement, CDims... > Coordinate
A Coordinate represents a coordinate in the continuous space.
Definition: coordinate.hpp:19
Definition: chunk.hpp:15
Definition: chunk_span.hpp:26
Definition: reducer.hpp:96
Definition: reducer.hpp:10