CartesianFluxDivergence.hpp
2 // See LICENSE.txt for details.
3
4 #pragma once
5
6 #include <cstddef>
7
8 #include "DataStructures/DataVector.hpp"
11 #include "Utilities/Gsl.hpp"
12
13 namespace evolution::dg::subcell {
14 // @{
15 /*!
16  * \brief Compute and add the 2nd-order flux divergence on a Cartesian mesh to
17  * the cell-centered time derivatives.
18  */
20  const double one_over_delta,
21  const DataVector& inv_jacobian,
22  const DataVector& boundary_correction,
23  const Index<1>& subcell_extents,
24  const size_t dimension) noexcept {
25  (void)dimension;
26  ASSERT(dimension == 0, "dimension must be 0 but is " << dimension);
27  for (size_t i = 0; i < subcell_extents[0]; ++i) {
28  (*dt_var)[i] += one_over_delta * inv_jacobian[i] *
29  (boundary_correction[i + 1] - boundary_correction[i]);
30  }
31 }
32
34  const double one_over_delta,
35  const DataVector& inv_jacobian,
36  const DataVector& boundary_correction,
37  const Index<2>& subcell_extents,
38  const size_t dimension) noexcept {
39  ASSERT(dimension == 0 or dimension == 1,
40  "dimension must be 0 or 1 but is " << dimension);
41  Index<2> subcell_face_extents = subcell_extents;
42  ++subcell_face_extents[dimension];
43  for (size_t j = 0; j < subcell_extents[1]; ++j) {
44  for (size_t i = 0; i < subcell_extents[0]; ++i) {
45  Index<2> index(i, j);
46  const size_t volume_index = collapsed_index(index, subcell_extents);
47  const size_t boundary_correction_lower_index =
48  collapsed_index(index, subcell_face_extents);
49  ++index[dimension];
50  const size_t boundary_correction_upper_index =
51  collapsed_index(index, subcell_face_extents);
52  (*dt_var)[volume_index] +=
53  one_over_delta * inv_jacobian[volume_index] *
54  (boundary_correction[boundary_correction_upper_index] -
55  boundary_correction[boundary_correction_lower_index]);
56  }
57  }
58 }
59
61  const double one_over_delta,
62  const DataVector& inv_jacobian,
63  const DataVector& boundary_correction,
64  const Index<3>& subcell_extents,
65  const size_t dimension) noexcept {
66  ASSERT(dimension == 0 or dimension == 1 or dimension == 2,
67  "dimension must be 0, 1, or 2 but is " << dimension);
68  Index<3> subcell_face_extents = subcell_extents;
69  ++subcell_face_extents[dimension];
70  for (size_t k = 0; k < subcell_extents[2]; ++k) {
71  for (size_t j = 0; j < subcell_extents[1]; ++j) {
72  for (size_t i = 0; i < subcell_extents[0]; ++i) {
73  Index<3> index(i, j, k);
74  const size_t volume_index = collapsed_index(index, subcell_extents);
75  const size_t boundary_correction_lower_index =
76  collapsed_index(index, subcell_face_extents);
77  ++index[dimension];
78  const size_t boundary_correction_upper_index =
79  collapsed_index(index, subcell_face_extents);
80  (*dt_var)[volume_index] +=
81  one_over_delta * inv_jacobian[volume_index] *
82  (boundary_correction[boundary_correction_upper_index] -
83  boundary_correction[boundary_correction_lower_index]);
84  }
85  }
86  }
87 }
88 // @}
89 } // namespace evolution::dg::subcell
Index
Definition: Index.hpp:31
collapsed_index
size_t collapsed_index(const Index< N > &index, const Index< N > &extents) noexcept
cstddef
Assert.hpp
Index.hpp
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49