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 = UniformPointSampling<X>;
24 
26 // Our second continuous dimension
27 struct Y;
28 // Its uniform discretization
29 using DDimY = UniformPointSampling<Y>;
31 
33 // Our simulated time dimension
34 struct T;
35 // Its uniform discretization
36 using DDimT = 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  ChunkSpan temp_slice
59  = temp[get_domain<DDimY>(temp).front()
60  + get_domain<DDimY>(temp).size() / 2];
61  std::cout << " * temperature[y:"
62  << get_domain<DDimY>(temp).size() / 2 << "] = {";
63  for_each(
65  get_domain<DDimX>(temp),
66  [=](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  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  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  = init_discrete_space(DDimX::init_ghosted(
118  Coordinate<X>(x_start),
119  Coordinate<X>(x_end),
120  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  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  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  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  = init_discrete_space(DDimY::init_ghosted(
144  Coordinate<Y>(y_start),
145  Coordinate<Y>(y_end),
146  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  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  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.,
166  [](DiscreteElement<DDimX> ix) {
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.,
175  [](DiscreteElement<DDimY> iy) {
176  return 1.
177  / (distance_at_left(iy) * distance_at_right(iy));
178  });
179  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  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  DiscreteDomain<DDimT> const time_domain
189  = init_discrete_space(DDimT::
190  init(Coordinate<T>(start_time),
191  Coordinate<T>(end_time),
192  nb_time_steps + 1));
194 
196  // Maps temperature into the full domain (including ghosts) twice:
197  // - once for the last fully computed time-step
198  Chunk ghosted_last_temp(
200  DDimX,
201  DDimY>(ghosted_x_domain, ghosted_y_domain),
203 
204  // - once for time-step being computed
205  Chunk ghosted_next_temp(
207  DDimX,
208  DDimY>(ghosted_x_domain, ghosted_y_domain),
211 
213  ChunkSpan const ghosted_initial_temp = ghosted_last_temp.span_view();
214  // Initialize the temperature on the main domain
215  for_each(
217  DiscreteDomain<DDimX, DDimY>(x_domain, y_domain),
218  DDC_LAMBDA(DiscreteElement<DDimX, DDimY> const ixy) {
219  double const x = coordinate(select<DDimX>(ixy));
220  double const y = coordinate(select<DDimY>(ixy));
221  ghosted_initial_temp(ixy)
222  = 9.999 * ((x * x + y * y) < 0.25);
223  });
225 
226  Chunk ghosted_temp(
228  DDimX,
229  DDimY>(ghosted_x_domain, ghosted_y_domain),
231 
232 
234  // display the initial data
235  deepcopy(ghosted_temp, ghosted_last_temp);
236  display(coordinate(time_domain.front()),
237  ghosted_temp[x_domain][y_domain]);
238  // time of the iteration where the last output happened
239  DiscreteElement<DDimT> last_output = time_domain.front();
241 
243  for (auto const iter :
244  time_domain.remove_first(DiscreteVector<DDimT>(1))) {
246 
248  // Periodic boundary conditions
249  deepcopy(
250  ghosted_last_temp[x_pre_ghost][y_domain],
251  ghosted_last_temp[y_domain][x_domain_end]);
252  deepcopy(
253  ghosted_last_temp[y_domain][x_post_ghost],
254  ghosted_last_temp[y_domain][x_domain_begin]);
255  deepcopy(
256  ghosted_last_temp[x_domain][y_pre_ghost],
257  ghosted_last_temp[x_domain][y_domain_end]);
258  deepcopy(
259  ghosted_last_temp[x_domain][y_post_ghost],
260  ghosted_last_temp[x_domain][y_domain_begin]);
262 
264  // a span excluding ghosts of the temperature at the time-step we
265  // will build
266  ChunkSpan const next_temp {
267  ghosted_next_temp[x_domain][y_domain]};
268  // a read-only view of the temperature at the previous time-step
269  ChunkSpan const last_temp {ghosted_last_temp.span_view()};
271 
273  // Stencil computation on the main domain
274  for_each(
276  next_temp.domain(),
277  DDC_LAMBDA(DiscreteElement<DDimX, DDimY> const ixy) {
278  DiscreteElement<DDimX> const ix = select<DDimX>(ixy);
279  DiscreteElement<DDimY> const iy = select<DDimY>(ixy);
280  double const dx_l = distance_at_left(ix);
281  double const dx_r = distance_at_right(ix);
282  double const dx_m = 0.5 * (dx_l + dx_r);
283  double const dy_l = distance_at_left(iy);
284  double const dy_r = distance_at_right(iy);
285  double const dy_m = 0.5 * (dy_l + dy_r);
286  next_temp(ix, iy) = last_temp(ix, iy);
287  next_temp(ix, iy)
288  += kx * max_dt
289  * (dx_l * last_temp(ix + 1, iy)
290  - 2.0 * dx_m * last_temp(ix, iy)
291  + dx_r * last_temp(ix - 1, iy))
292  / (dx_l * dx_m * dx_r);
293  next_temp(ix, iy)
294  += ky * max_dt
295  * (dy_l * last_temp(ix, iy + 1)
296  - 2.0 * dy_m * last_temp(ix, iy)
297  + dy_r * last_temp(ix, iy - 1))
298  / (dy_l * dy_m * dy_r);
299  });
301 
303  if (iter - last_output >= t_output_period) {
304  last_output = iter;
305  deepcopy(ghosted_temp, ghosted_last_temp);
306  display(coordinate(iter), ghosted_temp[x_domain][y_domain]);
307  }
309 
311  // Swap our two buffers
312  std::swap(ghosted_last_temp, ghosted_next_temp);
314  }
315 
317  if (last_output < time_domain.back()) {
318  deepcopy(ghosted_temp, ghosted_last_temp);
319  display(coordinate(time_domain.back()),
320  ghosted_temp[x_domain][y_domain]);
321  }
323 }
Definition: discrete_domain.hpp:22
constexpr discrete_element_type front() const noexcept
Definition: discrete_domain.hpp:125
constexpr discrete_element_type back() const noexcept
Definition: discrete_domain.hpp:130
constexpr DiscreteDomain remove_first(mlength_type n) const
Definition: discrete_domain.hpp:145
A DiscreteElement identifies an element of the discrete dimension.
Definition: discrete_element.hpp:123
A DiscreteVector is a vector in the discrete dimension.
Definition: discrete_vector.hpp:207
Definition: kokkos_allocator.hpp:12
Definition: scope_guard.hpp:12
UniformPointSampling models a uniform discretization of the provided continuous dimension.
Definition: uniform_point_sampling.hpp:21
Definition: chunk.hpp:12
Definition: chunk_span.hpp:24
constexpr parallel_device_policy parallel_device
Definition: for_each.hpp:167
constexpr serial_host_policy serial_host
Definition: for_each.hpp:165
Definition: reducer.hpp:96
Definition: reducer.hpp:10