Line data Source code
1 0 : // 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"
10 : #include "DataStructures/Tensor/Tensor.hpp"
11 : #include "DataStructures/Variables.hpp"
12 : #include "Domain/Structure/Direction.hpp"
13 : #include "Domain/Structure/Side.hpp"
14 : #include "NumericalAlgorithms/Spectral/BoundaryInterpolationMatrices.hpp"
15 : #include "NumericalAlgorithms/Spectral/BoundaryLiftingTerm.hpp"
16 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
17 : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
18 : #include "Utilities/Gsl.hpp"
19 :
20 : namespace dg {
21 : namespace detail {
22 : template <typename ValueType, size_t Dim>
23 : void lift_boundary_terms_gauss_points_impl(
24 : gsl::not_null<ValueType*> volume_dt_vars, size_t num_independent_components,
25 : const Mesh<Dim>& volume_mesh, size_t dimension,
26 : const Scalar<DataVector>& volume_det_inv_jacobian, size_t num_boundary_pts,
27 : const gsl::span<const ValueType>& boundary_corrections,
28 : const DataVector& boundary_lifting_term,
29 : const Scalar<DataVector>& magnitude_of_face_normal,
30 : const Scalar<DataVector>& face_det_jacobian);
31 :
32 : template <typename ValueType, size_t Dim>
33 : void lift_boundary_terms_gauss_points_impl(
34 : gsl::not_null<ValueType*> volume_dt_vars, size_t num_independent_components,
35 : const Mesh<Dim>& volume_mesh, size_t dimension,
36 : const Scalar<DataVector>& volume_det_inv_jacobian, size_t num_boundary_pts,
37 : const gsl::span<const ValueType>& upper_boundary_corrections,
38 : const DataVector& upper_boundary_lifting_term,
39 : const Scalar<DataVector>& upper_magnitude_of_face_normal,
40 : const Scalar<DataVector>& upper_face_det_jacobian,
41 : const gsl::span<const ValueType>& lower_boundary_corrections,
42 : const DataVector& lower_boundary_lifting_term,
43 : const Scalar<DataVector>& lower_magnitude_of_face_normal,
44 : const Scalar<DataVector>& lower_face_det_jacobian);
45 :
46 : template <typename ValueType, size_t Dim>
47 : void lift_boundary_terms_gauss_points_impl(
48 : gsl::not_null<ValueType*> volume_data, size_t num_independent_components,
49 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
50 : const gsl::span<const ValueType>& boundary_corrections);
51 : } // namespace detail
52 :
53 : /*!
54 : * \brief Lift the boundary corrections to the volume time derivatives in the
55 : * specified direction.
56 : *
57 : * The general lifting term (for the \f$\xi\f$-dimension) is:
58 : *
59 : * \f{align*}{
60 : * \partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}=\cdots
61 : * -\frac{\ell_{\breve{\imath}}\left(\xi=1\right)}
62 : * {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
63 : * \left[J\sqrt{
64 : * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
65 : * \frac{\partial\xi}{\partial x^j}}
66 : * \left(G_{\alpha} + D_{\alpha}\right)
67 : * \right]_{\breve{\jmath}\breve{k}}\left(\xi=1\right),
68 : * \f}
69 : *
70 : * where \f$\breve{\imath}\f$, \f$\breve{\jmath}\f$, and \f$\breve{k}\f$ are
71 : * indices in the logical \f$\xi\f$, \f$\eta\f$, and \f$\zeta\f$ dimensions.
72 : * The \f$G+D\f$ terms correspond to the `boundary_corrections` function
73 : * argument, and the function Spectral::boundary_lifting_term() is used to
74 : * compute and cache the terms from the lifting terms
75 : * \f$\ell_{\breve{\imath}}(\xi=\pm1)/w_{\breve{\imath}}\f$.
76 : *
77 : * \note that normal vectors are pointing out of the element.
78 : *
79 : * \param dt_vars The volume time derivatives
80 : * $\partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}$ to be updated
81 : * \param volume_det_inv_jacobian The inverse Jacobian determinant in the volume
82 : * $J^{-1}_{\breve{\imath}\breve{\jmath}\breve{k}}$
83 : * \param volume_mesh The mesh of the volume
84 : * \param direction The direction of the face on which the boundary corrections
85 : * are defined
86 : * \param boundary_corrections The boundary corrections
87 : * $\left(G_{\alpha} + D_{\alpha}\right)$
88 : * \param magnitude_of_face_normal The term $\sqrt{
89 : * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
90 : * \frac{\partial\xi}{\partial x^j}}$
91 : * \param face_det_jacobian The volume Jacobian determinant $J$ evaluated on the
92 : * face (not the surface Jacobian determinant)
93 : */
94 : template <size_t Dim, typename DtTagsList, typename BoundaryCorrectionTagsList>
95 1 : void lift_boundary_terms_gauss_points(
96 : const gsl::not_null<Variables<DtTagsList>*> dt_vars,
97 : const Scalar<DataVector>& volume_det_inv_jacobian,
98 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
99 : const Variables<BoundaryCorrectionTagsList>& boundary_corrections,
100 : const Scalar<DataVector>& magnitude_of_face_normal,
101 : const Scalar<DataVector>& face_det_jacobian) {
102 : ASSERT(volume_mesh.quadrature(direction.dimension()) ==
103 : Spectral::Quadrature::Gauss,
104 : "Must use Gauss points in the face-normal direction "
105 : << direction.dimension() << " but got mesh: " << volume_mesh);
106 : const Mesh<Dim - 1> boundary_mesh =
107 : volume_mesh.slice_away(direction.dimension());
108 : const Mesh<1> volume_stripe_mesh =
109 : volume_mesh.slice_through(direction.dimension());
110 : const size_t num_boundary_grid_points = boundary_mesh.number_of_grid_points();
111 : detail::lift_boundary_terms_gauss_points_impl(
112 : make_not_null(dt_vars->data()), dt_vars->number_of_independent_components,
113 : volume_mesh, direction.dimension(), volume_det_inv_jacobian,
114 : num_boundary_grid_points,
115 : gsl::make_span(boundary_corrections.data(), boundary_corrections.size()),
116 : direction.side() == Side::Upper
117 : ? Spectral::boundary_lifting_term(volume_stripe_mesh).second
118 : : Spectral::boundary_lifting_term(volume_stripe_mesh).first,
119 : magnitude_of_face_normal, face_det_jacobian);
120 : }
121 :
122 : /*!
123 : * \brief Lift both the upper and lower (in logical coordinates) boundary
124 : * corrections to the volume time derivatives in the specified logical
125 : * dimension.
126 : *
127 : * The upper and lower boundary corrections in the logical `dimension` are
128 : * lifted together in order to reduce the amount of striding through data that
129 : * is needed and to improve cache-friendliness.
130 : *
131 : * The general lifting term (for the \f$\xi\f$-dimension) is:
132 : *
133 : * \f{align*}{
134 : * \partial_t u_{\alpha\breve{\imath}\breve{\jmath}\breve{k}}=\cdots
135 : * -\frac{\ell_{\breve{\imath}}\left(\xi=1\right)}
136 : * {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
137 : * \left[J\sqrt{
138 : * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
139 : * \frac{\partial\xi}{\partial x^j}}
140 : * \left(G_{\alpha} + D_{\alpha}\right)
141 : * \right]_{\breve{\jmath}\breve{k}}\left(\xi=1\right)
142 : * - \frac{\ell_{\breve{\imath}}\left(\xi=-1\right)}
143 : * {w_{\breve{\imath}}J_{\breve{\imath}\breve{\jmath}\breve{k}}}
144 : * \left[J\sqrt{
145 : * \frac{\partial\xi}{\partial x^i} \gamma^{ij}
146 : * \frac{\partial\xi}{\partial x^j}}
147 : * \left(G_{\alpha} + D_{\alpha}\right)
148 : * \right]_{\breve{\jmath}\breve{k}}\left(\xi=-1\right),
149 : * \f}
150 : *
151 : * where \f$\breve{\imath}\f$, \f$\breve{\jmath}\f$, and \f$\breve{k}\f$ are
152 : * indices in the logical \f$\xi\f$, \f$\eta\f$, and \f$\zeta\f$ dimensions.
153 : * The \f$G+D\f$ terms correspond to the `upper_boundary_corrections` and
154 : * `lower_boundary_corrections` function arguments, and the function
155 : * Spectral::boundary_lifting_term() is used to compute and cache the terms
156 : * from the lifting terms
157 : * \f$\ell_{\breve{\imath}}(\xi=\pm1)/w_{\breve{\imath}}\f$.
158 : *
159 : * \note that normal vectors are pointing out of the element and therefore both
160 : * terms have the same sign.
161 : */
162 : template <size_t Dim, typename DtTagsList, typename BoundaryCorrectionTagsList>
163 1 : void lift_boundary_terms_gauss_points(
164 : const gsl::not_null<Variables<DtTagsList>*> dt_vars,
165 : const Scalar<DataVector>& volume_det_inv_jacobian,
166 : const Mesh<Dim>& volume_mesh, const size_t dimension,
167 : const Variables<BoundaryCorrectionTagsList>& upper_boundary_corrections,
168 : const Scalar<DataVector>& upper_magnitude_of_face_normal,
169 : const Scalar<DataVector>& upper_face_det_jacobian,
170 : const Variables<BoundaryCorrectionTagsList>& lower_boundary_corrections,
171 : const Scalar<DataVector>& lower_magnitude_of_face_normal,
172 : const Scalar<DataVector>& lower_face_det_jacobian) {
173 : ASSERT(volume_mesh.quadrature(dimension) == Spectral::Quadrature::Gauss,
174 : "Must use Gauss points in the face-normal dimension "
175 : << dimension << " but got mesh: " << volume_mesh);
176 : const Mesh<Dim - 1> boundary_mesh = volume_mesh.slice_away(dimension);
177 : const Mesh<1> volume_stripe_mesh = volume_mesh.slice_through(dimension);
178 : const size_t num_boundary_grid_points = boundary_mesh.number_of_grid_points();
179 : detail::lift_boundary_terms_gauss_points_impl(
180 : make_not_null(dt_vars->data()), dt_vars->number_of_independent_components,
181 : volume_mesh, dimension, volume_det_inv_jacobian, num_boundary_grid_points,
182 : gsl::make_span(upper_boundary_corrections.data(),
183 : upper_boundary_corrections.size()),
184 : Spectral::boundary_lifting_term(volume_stripe_mesh).second,
185 : upper_magnitude_of_face_normal, upper_face_det_jacobian,
186 : gsl::make_span(lower_boundary_corrections.data(),
187 : lower_boundary_corrections.size()),
188 : Spectral::boundary_lifting_term(volume_stripe_mesh).first,
189 : lower_magnitude_of_face_normal, lower_face_det_jacobian);
190 : }
191 :
192 : /*!
193 : * \brief Lift the boundary corrections from the face to the volume for Gauss
194 : * points in the direction perpendicular to the face
195 : *
196 : * This function doesn't apply any Jacobians or quadrature weights. It only
197 : * applies the lifting matrix $\ell_{\breve{\imath}}(\xi=\pm1)$ to the
198 : * `boundary_corrections` and adds the result to the `volume_data`.
199 : */
200 : template <size_t Dim, typename... VolumeTags,
201 : typename... BoundaryCorrectionTags>
202 1 : void lift_boundary_terms_gauss_points(
203 : const gsl::not_null<Variables<tmpl::list<VolumeTags...>>*> volume_data,
204 : const Variables<tmpl::list<BoundaryCorrectionTags...>>&
205 : boundary_corrections,
206 : const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction) {
207 : detail::lift_boundary_terms_gauss_points_impl(
208 : make_not_null(volume_data->data()),
209 : volume_data->number_of_independent_components, volume_mesh, direction,
210 : gsl::make_span(boundary_corrections.data(), boundary_corrections.size()));
211 : }
212 : } // namespace dg
|
