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