ThunderEgg  1.0.0
DFTPatchSolver.h
Go to the documentation of this file.
1 /***************************************************************************
2  * ThunderEgg, a library for solvers on adaptively refined block-structured
3  * Cartesian grids.
4  *
5  * Copyright (c) 2020-2021 Scott Aiton
6  *
7  * This program is free software: you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation, either version 3 of the License, or
10  * (at your option) any later version.
11  *
12  * This program is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with this program. If not, see <https://www.gnu.org/licenses/>.
19  ***************************************************************************/
20 
21 #ifndef THUNDEREGG_POISSON_DFTPATCHSOLVER_H
22 #define THUNDEREGG_POISSON_DFTPATCHSOLVER_H
23 
28 #include <ThunderEgg/Domain.h>
29 #include <ThunderEgg/PatchArray.h>
30 #include <ThunderEgg/PatchOperator.h>
31 #include <ThunderEgg/PatchSolver.h>
32 #include <ThunderEgg/Vector.h>
33 #include <bitset>
34 #include <map>
35 #include <valarray>
36 
37 extern "C" void
38 dgemv_(char&, int&, int&, double&, double*, int&, double*, int&, double&, double*, int&);
39 
40 namespace ThunderEgg::Poisson {
46 template<int D>
47 class DFTPatchSolver : public PatchSolver<D>
48 {
49 private:
53  enum DftType
54  {
55  DCT_II,
56  DCT_III,
57  DCT_IV,
58  DST_II,
59  DST_III,
60  DST_IV
61  };
66  using CompareFunction = std::function<bool(const PatchInfo<D>&, const PatchInfo<D>&)>;
70  std::shared_ptr<const PatchOperator<D>> op;
74  std::map<PatchInfo<D>, std::array<std::shared_ptr<std::valarray<double>>, D>, CompareFunction>
75  plan1;
79  std::map<PatchInfo<D>, std::array<std::shared_ptr<std::valarray<double>>, D>, CompareFunction>
80  plan2;
84  std::map<const PatchInfo<D>, PatchArray<D>, CompareFunction> eigen_vals;
88  std::map<std::tuple<DftType, int>, std::shared_ptr<std::valarray<double>>> transform_matrixes;
89 
93  std::bitset<Side<D>::number_of> neumann;
94 
103  bool patchIsNeumannOnSide(const PatchInfo<D>& pinfo, Side<D> s)
104  {
105  return !pinfo.hasNbr(s) && neumann[s.getIndex()];
106  }
114  std::array<std::shared_ptr<std::valarray<double>>, D> plan(std::array<DftType, D> types)
115  {
116  std::array<std::shared_ptr<std::valarray<double>>, D> retval;
117  for (size_t axis = 0; axis < D; axis++) {
118  retval[axis] = getTransformArray(types[axis], this->getDomain().getNs()[axis]);
119  }
120  return retval;
121  }
122 
130  std::shared_ptr<std::valarray<double>> getTransformArray(DftType type, int n)
131  {
132  std::shared_ptr<std::valarray<double>> matrix_ptr;
133 
134  if (transform_matrixes.count(std::make_tuple(type, n)) == 0) {
135  matrix_ptr = std::make_shared<std::valarray<double>>(n * n);
136  transform_matrixes[std::make_tuple(type, n)] = matrix_ptr;
137  std::valarray<double>& matrix = *matrix_ptr;
138 
139  switch (type) {
140  case DftType::DCT_II:
141  for (int j = 0; j < n; j++) {
142  for (int i = 0; i < n; i++) {
143  matrix[i * n + j] = cos(M_PI / n * (i * (j + 0.5)));
144  }
145  }
146  break;
147  case DftType::DCT_III:
148  for (int i = 0; i < n; i++) {
149  matrix[i * n] = 0.5;
150  }
151  for (int j = 1; j < n; j++) {
152  for (int i = 0; i < n; i++) {
153  matrix[i * n + j] = cos(M_PI / n * ((i + 0.5) * j));
154  }
155  }
156  break;
157  case DftType::DCT_IV:
158  for (int j = 0; j < n; j++) {
159  for (int i = 0; i < n; i++) {
160  matrix[i * n + j] = cos(M_PI / n * ((i + 0.5) * (j + 0.5)));
161  }
162  }
163  break;
164  case DftType::DST_II:
165  for (int i = 0; i < n; i++) {
166  for (int j = 0; j < n; j++) {
167  matrix[i * n + j] = sin(M_PI / n * ((i + 1) * (j + 0.5)));
168  }
169  }
170  break;
171  case DftType::DST_III:
172  for (int i = 0; i < n; i += 2) {
173  matrix[i * n + n - 1] = 0.5;
174  }
175  for (int i = 1; i < n; i += 2) {
176  matrix[i * n + n - 1] = -0.5;
177  }
178  for (int i = 0; i < n; i++) {
179  for (int j = 0; j < n - 1; j++) {
180  matrix[i * n + j] = sin(M_PI / n * ((i + 0.5) * (j + 1)));
181  }
182  }
183  break;
184  case DftType::DST_IV:
185  for (int j = 0; j < n; j++) {
186  for (int i = 0; i < n; i++) {
187  matrix[i * n + j] = sin(M_PI / n * ((i + 0.5) * (j + 0.5)));
188  }
189  }
190  }
191  } else {
192  matrix_ptr = transform_matrixes.at(std::make_tuple(type, n));
193  }
194  return matrix_ptr;
195  }
204  void executePlan(const std::array<std::shared_ptr<std::valarray<double>>, D>& plan,
205  const PatchView<double, D>& in,
206  const PatchView<double, D>& out) const
207  {
208  PatchView<double, D> prev_result = in;
209 
210  PatchArray<D> local_tmp;
211  if (!(D % 2)) {
212  local_tmp = PatchArray<D>(op->getDomain().getNs(), 1, 0);
213  }
214 
215  for (size_t axis = 0; axis < D; axis++) {
216  int n = this->getDomain().getNs()[axis];
217  std::array<int, D + 1> start = in.getStart();
218  std::array<int, D + 1> end = in.getEnd();
219  end[axis] = 0;
220 
221  std::valarray<double>& matrix = *plan[axis];
222 
223  PatchView<double, D> new_result;
224  if (D % 2) {
225  if (axis % 2) {
226  new_result = in;
227  } else {
228  new_result = out;
229  }
230  } else {
231  if (axis == D - 1) {
232  new_result = out;
233  } else if (axis == D - 2) {
234  new_result = local_tmp.getView();
235  } else if (axis % 2) {
236  new_result = in;
237  } else {
238  new_result = out;
239  }
240  }
241 
242  int pstride = prev_result.getStrides()[axis];
243  int nstride = new_result.getStrides()[axis];
244 
245  char T = 'T';
246  double one = 1;
247  double zero = 0;
248  Loop::Nested<D + 1>(start, end, [&](std::array<int, D + 1> coord) {
249  dgemv_(T,
250  n,
251  n,
252  one,
253  &matrix[0],
254  n,
255  &prev_result[coord],
256  pstride,
257  zero,
258  &new_result[coord],
259  nstride);
260  });
261 
262  prev_result = new_result;
263  }
264  }
271  std::array<DftType, D> getTransformsForPatch(const ThunderEgg::PatchInfo<D>& pinfo)
272  {
273  // get transform types for each axis
274  std::array<DftType, D> transforms;
275  for (size_t axis = 0; axis < D; axis++) {
276  if (patchIsNeumannOnSide(pinfo, LowerSideOnAxis<D>(axis)) &&
277  patchIsNeumannOnSide(pinfo, HigherSideOnAxis<D>(axis))) {
278  transforms[axis] = DftType::DCT_II;
279  } else if (patchIsNeumannOnSide(pinfo, LowerSideOnAxis<D>(axis))) {
280  transforms[axis] = DftType::DCT_IV;
281  } else if (patchIsNeumannOnSide(pinfo, HigherSideOnAxis<D>(axis))) {
282  transforms[axis] = DftType::DST_IV;
283  } else {
284  transforms[axis] = DftType::DST_II;
285  }
286  }
287  return transforms;
288  }
295  std::array<DftType, D> getInverseTransformsForPatch(const ThunderEgg::PatchInfo<D>& pinfo)
296  {
297  // get transform types for each axis
298  std::array<DftType, D> transforms_inv;
299  for (size_t axis = 0; axis < D; axis++) {
300  if (patchIsNeumannOnSide(pinfo, LowerSideOnAxis<D>(axis)) &&
301  patchIsNeumannOnSide(pinfo, HigherSideOnAxis<D>(axis))) {
302  transforms_inv[axis] = DftType::DCT_III;
303  } else if (patchIsNeumannOnSide(pinfo, LowerSideOnAxis<D>(axis))) {
304  transforms_inv[axis] = DftType::DCT_IV;
305  } else if (patchIsNeumannOnSide(pinfo, HigherSideOnAxis<D>(axis))) {
306  transforms_inv[axis] = DftType::DST_IV;
307  } else {
308  transforms_inv[axis] = DftType::DST_III;
309  }
310  }
311  return transforms_inv;
312  }
319  PatchArray<D> getEigenValues(const PatchInfo<D>& pinfo)
320  {
321  PatchArray<D> retval(this->getDomain().getNs(), 1, 0);
322 
323  std::valarray<size_t> all_strides(D);
324  size_t curr_stride = 1;
325  for (size_t i = 0; i < D; i++) {
326  all_strides[i] = curr_stride;
327  curr_stride *= pinfo.ns[i];
328  }
329 
330  for (size_t axis = 0; axis < D; axis++) {
331  int n = pinfo.ns[axis];
332  double h = pinfo.spacings[axis];
333 
334  if (patchIsNeumannOnSide(pinfo, LowerSideOnAxis<D>(axis)) &&
335  patchIsNeumannOnSide(pinfo, HigherSideOnAxis<D>(axis))) {
336  for (int xi = 0; xi < n; xi++) {
337  double val = 4 / (h * h) * pow(sin(xi * M_PI / (2 * n)), 2);
338  View<double, D> slice = retval.getSliceOn(Side<D>(2 * axis), { xi });
339  Loop::OverInteriorIndexes<D>(
340  slice, [&](const std::array<int, D>& coord) { slice[coord] -= val; });
341  }
342  } else if (patchIsNeumannOnSide(pinfo, LowerSideOnAxis<D>(axis)) ||
343  patchIsNeumannOnSide(pinfo, HigherSideOnAxis<D>(axis))) {
344  for (int xi = 0; xi < n; xi++) {
345  double val = 4 / (h * h) * pow(sin((xi + 0.5) * M_PI / (2 * n)), 2);
346  View<double, D> slice = retval.getSliceOn(Side<D>(2 * axis), { xi });
347  Loop::OverInteriorIndexes<D>(
348  slice, [&](const std::array<int, D>& coord) { slice[coord] -= val; });
349  }
350  } else {
351  for (int xi = 0; xi < n; xi++) {
352  double val = 4 / (h * h) * pow(sin((xi + 1) * M_PI / (2 * n)), 2);
353  View<double, D> slice = retval.getSliceOn(Side<D>(2 * axis), { xi });
354  Loop::OverInteriorIndexes<D>(
355  slice, [&](const std::array<int, D>& coord) { slice[coord] -= val; });
356  }
357  }
358  }
359  return retval;
360  }
368  void addPatch(const PatchInfo<D>& pinfo)
369  {
370  if (plan1.count(pinfo) == 0) {
371  plan1[pinfo] = plan(getTransformsForPatch(pinfo));
372  plan2[pinfo] = plan(getInverseTransformsForPatch(pinfo));
373  eigen_vals.emplace(pinfo, getEigenValues(pinfo));
374  }
375  }
376 
377 public:
384  DFTPatchSolver(const PatchOperator<D>& op, std::bitset<Side<D>::number_of> neumann)
385  : PatchSolver<D>(op.getDomain(), op.getGhostFiller())
386  , op(op.clone())
387  , neumann(neumann)
388  {
389  CompareFunction compare = [&](const PatchInfo<D>& a, const PatchInfo<D>& b) {
390  std::bitset<Side<D>::number_of> a_neumann;
391  std::bitset<Side<D>::number_of> b_neumann;
392  for (Side<D> s : Side<D>::getValues()) {
393  a_neumann[s.getIndex()] = patchIsNeumannOnSide(a, s);
394  b_neumann[s.getIndex()] = patchIsNeumannOnSide(b, s);
395  }
396  return std::forward_as_tuple(a_neumann.to_ulong(), a.spacings[0]) <
397  std::forward_as_tuple(b_neumann.to_ulong(), b.spacings[0]);
398  };
399 
400  plan1 = std::map<PatchInfo<D>,
401  std::array<std::shared_ptr<std::valarray<double>>, D>,
402  CompareFunction>(compare);
403  plan2 = std::map<PatchInfo<D>,
404  std::array<std::shared_ptr<std::valarray<double>>, D>,
405  CompareFunction>(compare);
406  eigen_vals = std::map<const PatchInfo<D>, PatchArray<D>, CompareFunction>(compare);
407 
408  // process patches
409  for (auto pinfo : this->getDomain().getPatchInfoVector()) {
410  addPatch(pinfo);
411  }
412  }
418  DFTPatchSolver<D>* clone() const override { return new DFTPatchSolver<D>(*this); }
419  void solveSinglePatch(const PatchInfo<D>& pinfo,
420  const PatchView<const double, D>& f_view,
421  const PatchView<double, D>& u_view) const override
422  {
423  PatchArray<D> f_copy(op->getDomain().getNs(), 1, 0);
424  PatchArray<D> tmp(op->getDomain().getNs(), 1, 0);
425 
426  Loop::OverInteriorIndexes<D + 1>(
427  f_view, [&](std::array<int, D + 1> coord) { f_copy[coord] = f_view[coord]; });
428 
429  op->modifyRHSForInternalBoundaryConditions(pinfo, u_view, f_copy.getView());
430 
431  executePlan(plan1.at(pinfo), f_copy.getView(), tmp.getView());
432 
433  const PatchArray<D>& eigen_vals_view = eigen_vals.at(pinfo);
434  Loop::OverInteriorIndexes<D + 1>(
435  tmp, [&](std::array<int, D + 1> coord) { tmp[coord] /= eigen_vals_view[coord]; });
436 
437  if (neumann.all() && !pinfo.hasNbr()) {
438  tmp[tmp.getStart()] = 0;
439  }
440 
441  executePlan(plan2.at(pinfo), tmp.getView(), u_view);
442 
443  double scale = 1;
444  for (size_t axis = 0; axis < D; axis++) {
445  scale *= 2.0 / this->getDomain().getNs()[axis];
446  }
447  Loop::OverInteriorIndexes<D + 1>(u_view,
448  [&](std::array<int, D + 1> coord) { u_view[coord] *= scale; });
449  }
450 };
451 
452 extern template class DFTPatchSolver<2>;
453 extern template class DFTPatchSolver<3>;
454 } // namespace ThunderEgg::Poisson
455 #endif
ThunderEgg::Poisson::DFTPatchSolver::clone
DFTPatchSolver< D > * clone() const override
Clone this patch solver.
Definition: DFTPatchSolver.h:418
ThunderEgg::PatchSolver::getDomain
const Domain< D > & getDomain() const
Get the Domain object.
Definition: PatchSolver.h:84
ThunderEgg::View< double, D >
ThunderEgg::Poisson::DFTPatchSolver::solveSinglePatch
void solveSinglePatch(const PatchInfo< D > &pinfo, const PatchView< const double, D > &f_view, const PatchView< double, D > &u_view) const override
Perform a single solve over a patch.
Definition: DFTPatchSolver.h:419
Vector.h
Vector class.
ThunderEgg::View::getStrides
const std::array< int, D > & getStrides() const
Get the strides of the patch in each direction.
Definition: View.h:203
ThunderEgg::PatchSolver
Solves the problem on the patches using a specified interface value.
Definition: PatchSolver.h:42
ThunderEgg::PatchInfo::hasNbr
bool hasNbr(Face< D, M > s) const
Return whether the patch has a neighbor.
Definition: PatchInfo.h:284
ThunderEgg::PatchInfo::ns
std::array< int, D > ns
The number of cells in each direction.
Definition: PatchInfo.h:117
ThunderEgg::PatchArray::getSliceOn
View< double, M+1 > getSliceOn(Face< D, M > f, const std::array< int, D - M > &offset)
Get the slice on a given face.
Definition: PatchArray.h:116
ThunderEgg::PatchArray::getView
const PatchView< double, D > & getView()
Get the View for the array.
Definition: PatchArray.h:196
PatchOperator.h
PatchOperator class.
ThunderEgg::Poisson
Classes specific to the Poisson equation.
Definition: DFTPatchSolver.h:40
ThunderEgg::PatchView< double, D >
PatchArray.h
PatchArray class.
ThunderEgg::Poisson::DFTPatchSolver
Use DFT transforms to solve for the Poisson equation.
Definition: DFTPatchSolver.h:47
ThunderEgg::Face::getIndex
size_t getIndex() const
Get the index for this Face.
Definition: Face.h:452
ThunderEgg::View::getEnd
const std::array< int, D > & getEnd() const
Get the coordinate of the last element.
Definition: View.h:211
ThunderEgg::PatchArray
Array for acessing data of a patch. It supports variable striding.
Definition: PatchArray.h:36
ThunderEgg::PatchInfo
Contains metadata for a patch.
Definition: PatchInfo.h:51
ThunderEgg::PatchSolver::getGhostFiller
const GhostFiller< D > & getGhostFiller() const
Get the GhostFiller object.
Definition: PatchSolver.h:90
Domain.h
Domain class.
ThunderEgg::PatchOperator
This is an Operator where derived classes only have to implement the two virtual functions that opera...
Definition: PatchOperator.h:40
ThunderEgg::PatchInfo::spacings
std::array< double, D > spacings
The cell spacings in each direction.
Definition: PatchInfo.h:125
ThunderEgg::Face
Enum-style class for the faces of an n-dimensional cube.
Definition: Face.h:41
PatchSolver.h
PatchSolver class.
ThunderEgg::PatchArray::getStart
const std::array< int, D+1 > & getStart() const
Get the coordinate of the first element.
Definition: PatchArray.h:178
ThunderEgg::View::getStart
const std::array< int, D > & getStart() const
Get the coordinate of the first element.
Definition: View.h:207
ThunderEgg::Poisson::DFTPatchSolver::DFTPatchSolver
DFTPatchSolver(const PatchOperator< D > &op, std::bitset< Side< D >::number_of > neumann)
Construct a new DFTPatchSolver object.
Definition: DFTPatchSolver.h:384