LiftFromBoundary.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <utility>
8 
9 #include "DataStructures/DataVector.hpp"
16 #include "Utilities/Gsl.hpp"
17 
18 namespace evolution::dg {
19 namespace detail {
20 template <size_t Dim>
21 void lift_boundary_terms_gauss_points_impl(
22  gsl::not_null<double*> volume_dt_vars, size_t num_independent_components,
23  const Mesh<Dim>& volume_mesh, size_t dimension,
24  const Scalar<DataVector>& volume_det_inv_jacobian, size_t num_boundary_pts,
25  const gsl::span<const double>& boundary_corrections,
27  const Scalar<DataVector>& magnitude_of_face_normal,
28  const Scalar<DataVector>& face_det_jacobian) noexcept;
29 
30 template <size_t Dim>
31 void lift_boundary_terms_gauss_points_impl(
32  gsl::not_null<double*> volume_dt_vars, size_t num_independent_components,
33  const Mesh<Dim>& volume_mesh, size_t dimension,
34  const Scalar<DataVector>& volume_det_inv_jacobian, size_t num_boundary_pts,
35  const gsl::span<const double>& upper_boundary_corrections,
36  const DataVector& upper_boundary_lifting_term,
37  const Scalar<DataVector>& upper_magnitude_of_face_normal,
38  const Scalar<DataVector>& upper_face_det_jacobian,
39  const gsl::span<const double>& lower_boundary_corrections,
40  const DataVector& lower_boundary_lifting_term,
41  const Scalar<DataVector>& lower_magnitude_of_face_normal,
42  const Scalar<DataVector>& lower_face_det_jacobian) noexcept;
43 } // namespace detail
44 
45 /*!
46  * \brief Lift the boundary corrections to the volume time derivatives in the
47  * specified direction.
48  *
49  * The general lifting term (for the \f$\xi\f$-dimension) is:
50  *
51  * \f{align*}{
52  * \partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}=\cdots
53  * -\frac{\ell_{\breve{\imath}}\left(\xi=1\right)}
54  * {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
55  * \left[J\sqrt{
56  * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
57  * \frac{\partial\xi}{\partial x^j}}
58  * \left(G_{\alpha} + D_{\alpha}\right)
59  * \right]_{\breve{\jmath}\breve{k}}\left(\xi=1\right),
60  * \f}
61  *
62  * where \f$\breve{\imath}\f$, \f$\breve{\jmath}\f$, and \f$\breve{k}\f$ are
63  * indices in the logical \f$\xi\f$, \f$\eta\f$, and \f$\zeta\f$ dimensions.
64  * The \f$G+D\f$ terms correspond to the `boundary_corrections` function
65  * argument, and the function Spectral::boundary_lifting_term() is used to
66  * compute and cache the terms from the lifting terms
67  * \f$\ell_{\breve{\imath}}(\xi=\pm1)/w_{\breve{\imath}}\f$.
68  *
69  * \note that normal vectors are pointing out of the element.
70  */
71 template <size_t Dim, typename DtTagsList, typename BoundaryCorrectionTagsList>
73  const gsl::not_null<Variables<DtTagsList>*> dt_vars,
74  const Scalar<DataVector>& volume_det_inv_jacobian,
75  const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
76  const Variables<BoundaryCorrectionTagsList>& boundary_corrections,
77  const Scalar<DataVector>& magnitude_of_face_normal,
78  const Scalar<DataVector>& face_det_jacobian) noexcept {
79  ASSERT(std::all_of(volume_mesh.quadrature().begin(),
80  volume_mesh.quadrature().end(),
81  [](const Spectral::Quadrature quadrature) noexcept {
82  return quadrature == Spectral::Quadrature::Gauss;
83  }),
84  "Must use Gauss points in all directions but got the mesh: "
85  << volume_mesh);
86  const Mesh<Dim - 1> boundary_mesh =
87  volume_mesh.slice_away(direction.dimension());
88  const Mesh<1> volume_stripe_mesh =
89  volume_mesh.slice_through(direction.dimension());
90  const size_t num_boundary_grid_points = boundary_mesh.number_of_grid_points();
91  detail::lift_boundary_terms_gauss_points_impl(
92  make_not_null(dt_vars->data()), dt_vars->number_of_independent_components,
93  volume_mesh, direction.dimension(), volume_det_inv_jacobian,
94  num_boundary_grid_points,
95  gsl::make_span(boundary_corrections.data(), boundary_corrections.size()),
96  direction.side() == Side::Upper
97  ? Spectral::boundary_lifting_term(volume_stripe_mesh).second
98  : Spectral::boundary_lifting_term(volume_stripe_mesh).first,
99  magnitude_of_face_normal, face_det_jacobian);
100 }
101 
102 /*!
103  * \brief Lift both the upper and lower (in logical coordinates) boundary
104  * corrections to the volume time derivatives in the specified logical
105  * dimension.
106  *
107  * The upper and lower boundary corrections in the logical `dimension` are
108  * lifted together in order to reduce the amount of striding through data that
109  * is needed and to improve cache-friendliness.
110  *
111  * The general lifting term (for the \f$\xi\f$-dimension) is:
112  *
113  * \f{align*}{
114  * \partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}=\cdots
115  * -\frac{\ell_{\breve{\imath}}\left(\xi=1\right)}
116  * {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
117  * \left[J\sqrt{
118  * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
119  * \frac{\partial\xi}{\partial x^j}}
120  * \left(G_{\alpha} + D_{\alpha}\right)
121  * \right]_{\breve{\jmath}\breve{k}}\left(\xi=1\right)
122  * - \frac{\ell_{\breve{\imath}}\left(\xi=-1\right)}
123  * {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
124  * \left[J\sqrt{
125  * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
126  * \frac{\partial\xi}{\partial x^j}}
127  * \left(G_{\alpha} + D_{\alpha}\right)
128  * \right]_{\breve{\jmath}\breve{k}}\left(\xi=-1\right),
129  * \f}
130  *
131  * where \f$\breve{\imath}\f$, \f$\breve{\jmath}\f$, and \f$\breve{k}\f$ are
132  * indices in the logical \f$\xi\f$, \f$\eta\f$, and \f$\zeta\f$ dimensions.
133  * The \f$G+D\f$ terms correspond to the `upper_boundary_corrections` and
134  * `lower_boundary_corrections` function arguments, and the function
135  * Spectral::boundary_lifting_term() is used to compute and cache the terms
136  * from the lifting terms
137  * \f$\ell_{\breve{\imath}}(\xi=\pm1)/w_{\breve{\imath}}\f$.
138  *
139  * \note that normal vectors are pointing out of the element and therefore both
140  * terms have the same sign.
141  */
142 template <size_t Dim, typename DtTagsList, typename BoundaryCorrectionTagsList>
144  const gsl::not_null<Variables<DtTagsList>*> dt_vars,
145  const Scalar<DataVector>& volume_det_inv_jacobian,
146  const Mesh<Dim>& volume_mesh, const size_t dimension,
147  const Variables<BoundaryCorrectionTagsList>& upper_boundary_corrections,
148  const Scalar<DataVector>& upper_magnitude_of_face_normal,
149  const Scalar<DataVector>& upper_face_det_jacobian,
150  const Variables<BoundaryCorrectionTagsList>& lower_boundary_corrections,
151  const Scalar<DataVector>& lower_magnitude_of_face_normal,
152  const Scalar<DataVector>& lower_face_det_jacobian) noexcept {
153  ASSERT(std::all_of(volume_mesh.quadrature().begin(),
154  volume_mesh.quadrature().end(),
155  [](const Spectral::Quadrature quadrature) noexcept {
156  return quadrature == Spectral::Quadrature::Gauss;
157  }),
158  "Must use Gauss points in all directions but got the mesh: "
159  << volume_mesh);
160  const Mesh<Dim - 1> boundary_mesh = volume_mesh.slice_away(dimension);
161  const Mesh<1> volume_stripe_mesh = volume_mesh.slice_through(dimension);
162  const size_t num_boundary_grid_points = boundary_mesh.number_of_grid_points();
163  detail::lift_boundary_terms_gauss_points_impl(
164  make_not_null(dt_vars->data()), dt_vars->number_of_independent_components,
165  volume_mesh, dimension, volume_det_inv_jacobian, num_boundary_grid_points,
166  gsl::make_span(upper_boundary_corrections.data(),
167  upper_boundary_corrections.size()),
168  Spectral::boundary_lifting_term(volume_stripe_mesh).second,
169  upper_magnitude_of_face_normal, upper_face_det_jacobian,
170  gsl::make_span(lower_boundary_corrections.data(),
171  lower_boundary_corrections.size()),
172  Spectral::boundary_lifting_term(volume_stripe_mesh).first,
173  lower_magnitude_of_face_normal, lower_face_det_jacobian);
174 }
175 } // namespace evolution::dg
utility
evolution::dg
Functionality for evolving hyperbolic partial differential equations using the discontinuous Galerkin...
Definition: ConservativeDuDt.hpp:22
Side.hpp
evolution::dg::lift_boundary_terms_gauss_points
void lift_boundary_terms_gauss_points(const gsl::not_null< Variables< DtTagsList > * > dt_vars, const Scalar< DataVector > &volume_det_inv_jacobian, const Mesh< Dim > &volume_mesh, const Direction< Dim > &direction, const Variables< BoundaryCorrectionTagsList > &boundary_corrections, const Scalar< DataVector > &magnitude_of_face_normal, const Scalar< DataVector > &face_det_jacobian) noexcept
Lift the boundary corrections to the volume time derivatives in the specified direction.
Definition: LiftFromBoundary.hpp:72
Spectral.hpp
Direction< Dim >
Spectral::boundary_lifting_term
const std::pair< DataVector, DataVector > & boundary_lifting_term(const Mesh< 1 > &mesh) noexcept
Terms used during the lifting portion of a discontinuous Galerkin scheme when using Gauss points.
cstddef
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:46
Spectral::Quadrature
Quadrature
The choice of quadrature method to compute integration weights.
Definition: Spectral.hpp:90
Mesh::slice_through
Mesh< sizeof...(D)> slice_through(D... d) const noexcept
Returns a Mesh with the dimensions d, ... present (zero-indexed).
Definition: Mesh.hpp:227
ASSERT
#define ASSERT(a, m)
Assert that an expression should be true.
Definition: Assert.hpp:49
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
gsl::span
Create a span/view on a range, which is cheap to copy (one pointer).
Definition: Gsl.hpp:292
gsl::make_span
constexpr span< ElementType > make_span(ElementType *ptr, typename span< ElementType >::index_type count)
Definition: Gsl.hpp:801
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Gsl.hpp
Tensor.hpp
make_not_null
gsl::not_null< T * > make_not_null(T *ptr) noexcept
Construct a not_null from a pointer. Often this will be done as an implicit conversion,...
Definition: Gsl.hpp:880
Direction.hpp
Mesh::slice_away
Mesh< Dim - 1 > slice_away(size_t d) const noexcept
Returns a Mesh with dimension d removed (zero-indexed).
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: ReadSpecPiecewisePolynomial.hpp:13