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 <type_traits>
8 :
9 : #include "DataStructures/Matrix.hpp"
10 : #include "DataStructures/Tensor/Tensor.hpp"
11 : #include "DataStructures/Transpose.hpp"
12 : #include "DataStructures/Variables.hpp"
13 : #include "NumericalAlgorithms/LinearOperators/Divergence.hpp"
14 : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
15 : #include "NumericalAlgorithms/Spectral/Spectral.hpp"
16 : #include "Utilities/Blas.hpp"
17 : #include "Utilities/ContainerHelpers.hpp"
18 : #include "Utilities/Gsl.hpp"
19 : #include "Utilities/Literals.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /*!
23 : * \ingroup NumericalAlgorithmsGroup
24 : * \brief Compute the weak form divergence of fluxes
25 : *
26 : * In a discontinuous Galerkin scheme we integrate the equations against the
27 : * basis functions over the element. For the flux divergence term this gives:
28 : *
29 : * \f{align*}{
30 : * \int_{\Omega}d^n x \phi_{\breve{\imath}}\partial_i F^i,
31 : * \f}
32 : *
33 : * where the basis functions are denoted by \f$\phi_{\breve{\imath}}\f$.
34 : *
35 : * Integrating by parts we get
36 : *
37 : * \f{align*}{
38 : * \int_{\Omega}d^n x\, \phi_{\breve{\imath}}\partial_i F^i
39 : * = -\int_{\Omega}d^n x\, F^i \partial_i \phi_{\breve{\imath}}
40 : * + \int_{\partial\Omega} d^{(n-1)}\Sigma\, n_i F^i \phi_{\breve{\imath}}
41 : * \f}
42 : *
43 : * Next we expand the flux \f$F^i\f$ in terms of the basis functions, yielding
44 : *
45 : * \f{align*}{
46 : * - \int_{\Omega}d^n x\,F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \partial_i
47 : * \phi_{\breve{\imath}}
48 : * + \int_{\partial\Omega} d^{(n-1)}\Sigma\, n_i F^i_{\breve{\jmath}}
49 : * \phi_{\breve{\jmath}} \phi_{\breve{\imath}}
50 : * \f}
51 : *
52 : * This function computes the volume term:
53 : *
54 : * \f{align*}{
55 : * \frac{1}{w_\breve{\imath}} \int_{\Omega}d^n x \,
56 : * F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \partial_i \phi_{\breve{\imath}}
57 : * \f}
58 : *
59 : * \note When using Gauss-Lobatto points the numerical values of the
60 : * divergence are the same for the strong and weak divergence at the interior
61 : * points. When using Gauss points they are only the same at the central grid
62 : * point and only when an odd number of grid points is used.
63 : */
64 : template <typename... ResultTags, typename... FluxTags, size_t Dim>
65 1 : void weak_divergence(
66 : const gsl::not_null<Variables<tmpl::list<ResultTags...>>*>
67 : divergence_of_fluxes,
68 : const Variables<tmpl::list<FluxTags...>>& fluxes, const Mesh<Dim>& mesh,
69 : const InverseJacobian<DataVector, Dim, Frame::ElementLogical,
70 : Frame::Inertial>& det_jac_times_inverse_jacobian) {
71 : if (UNLIKELY(divergence_of_fluxes->number_of_grid_points() !=
72 : fluxes.number_of_grid_points())) {
73 : divergence_of_fluxes->initialize(fluxes.number_of_grid_points());
74 : }
75 :
76 : const auto apply_matrix_in_first_dim =
77 : [](double* result, const double* const input, const Matrix& matrix,
78 : const size_t size, const bool add_to_result) {
79 : dgemm_<true>(
80 : 'N', 'N',
81 : matrix.rows(), // rows of matrix and result
82 : size / matrix.columns(), // columns of result and input
83 : matrix.columns(), // columns of matrix and rows of input
84 : 1.0, // overall multiplier
85 : matrix.data(), // matrix
86 : matrix.spacing(), // rows of matrix including padding
87 : input, // input
88 : matrix.columns(), // rows of input
89 : add_to_result
90 : ? 1.0
91 : : 0.0, // 1.0 means add to result, 0.0 means overwrite result
92 : result, // result
93 : matrix.rows()); // rows of result
94 : };
95 :
96 : if constexpr (Dim == 1) {
97 : (void)det_jac_times_inverse_jacobian; // is identically 1.0 in 1d
98 :
99 : apply_matrix_in_first_dim(divergence_of_fluxes->data(), fluxes.data(),
100 : Spectral::weak_flux_differentiation_matrix(mesh),
101 : fluxes.size(), false);
102 : } else {
103 : // Multiplies the flux by det_jac_time_inverse_jacobian.
104 : const auto transform_to_logical_frame =
105 : [&det_jac_times_inverse_jacobian, &fluxes](
106 : auto flux_tag_v, auto div_tag_v,
107 : const gsl::not_null<Variables<tmpl::list<ResultTags...>>*>
108 : result_buffer,
109 : auto logical_index_of_jacobian) {
110 : using flux_tag = tmpl::type_from<decltype(flux_tag_v)>;
111 : using div_tag = tmpl::type_from<decltype(div_tag_v)>;
112 : using TensorIndexType =
113 : std::array<size_t, div_tag::type::num_tensor_indices>;
114 : using FluxIndexType =
115 : std::array<size_t, flux_tag::type::num_tensor_indices>;
116 :
117 : auto& result = get<div_tag>(*result_buffer);
118 : const auto& flux = get<flux_tag>(fluxes);
119 :
120 : for (size_t result_storage_index = 0;
121 : result_storage_index < result.size(); ++result_storage_index) {
122 : const TensorIndexType result_tensor_index =
123 : result.get_tensor_index(result_storage_index);
124 : const FluxIndexType flux_x_tensor_index =
125 : prepend(result_tensor_index, 0_st);
126 : const FluxIndexType flux_y_tensor_index =
127 : prepend(result_tensor_index, 1_st);
128 : if constexpr (Dim == 2) {
129 : result[result_storage_index] =
130 : get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
131 : 0>(det_jac_times_inverse_jacobian) *
132 : flux.get(flux_x_tensor_index) +
133 : get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
134 : 1>(det_jac_times_inverse_jacobian) *
135 : flux.get(flux_y_tensor_index);
136 : } else {
137 : const FluxIndexType flux_z_tensor_index =
138 : prepend(result_tensor_index, 2_st);
139 : result[result_storage_index] =
140 : get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
141 : 0>(det_jac_times_inverse_jacobian) *
142 : flux.get(flux_x_tensor_index) +
143 : get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
144 : 1>(det_jac_times_inverse_jacobian) *
145 : flux.get(flux_y_tensor_index) +
146 : get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
147 : 2>(det_jac_times_inverse_jacobian) *
148 : flux.get(flux_z_tensor_index);
149 : }
150 : }
151 : };
152 :
153 : if constexpr (Dim == 2) {
154 : Variables<tmpl::list<ResultTags...>> data_buffer{
155 : divergence_of_fluxes->number_of_grid_points()};
156 :
157 : const Matrix& eta_weak_div_matrix =
158 : Spectral::weak_flux_differentiation_matrix(mesh.slice_through(1));
159 : const Matrix& xi_weak_div_matrix =
160 : Spectral::weak_flux_differentiation_matrix(mesh.slice_through(0));
161 :
162 : // Compute the eta divergence term. Since that also needs a transpose,
163 : // copy into result, then transpose into `data_buffer`
164 : EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
165 : tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
166 : make_not_null(&data_buffer), std::integral_constant<size_t, 1>{}));
167 : double* div_ptr = divergence_of_fluxes->data();
168 : raw_transpose(make_not_null(div_ptr), data_buffer.data(),
169 : xi_weak_div_matrix.rows(),
170 : divergence_of_fluxes->size() / xi_weak_div_matrix.rows());
171 : apply_matrix_in_first_dim(data_buffer.data(),
172 : divergence_of_fluxes->data(),
173 : eta_weak_div_matrix, data_buffer.size(), false);
174 :
175 : const size_t chunk_size = Variables<tmpl::list<Tags::div<FluxTags>...>>::
176 : number_of_independent_components *
177 : eta_weak_div_matrix.rows();
178 : raw_transpose(make_not_null(div_ptr), data_buffer.data(), chunk_size,
179 : data_buffer.size() / chunk_size);
180 :
181 : // Now compute xi divergence and *add* to eta divergence
182 : EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
183 : tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
184 : make_not_null(&data_buffer), std::integral_constant<size_t, 0>{}));
185 : apply_matrix_in_first_dim(divergence_of_fluxes->data(),
186 : data_buffer.data(), xi_weak_div_matrix,
187 : data_buffer.size(), true);
188 : } else if constexpr (Dim == 3) {
189 : Variables<tmpl::list<ResultTags...>> data_buffer0{
190 : divergence_of_fluxes->number_of_grid_points()};
191 : Variables<tmpl::list<ResultTags...>> data_buffer1{
192 : divergence_of_fluxes->number_of_grid_points()};
193 : constexpr size_t number_of_independent_components =
194 : decltype(data_buffer1)::number_of_independent_components;
195 :
196 : const Matrix& zeta_weak_div_matrix =
197 : Spectral::weak_flux_differentiation_matrix(mesh.slice_through(2));
198 : const Matrix& eta_weak_div_matrix =
199 : Spectral::weak_flux_differentiation_matrix(mesh.slice_through(1));
200 : const Matrix& xi_weak_div_matrix =
201 : Spectral::weak_flux_differentiation_matrix(mesh.slice_through(0));
202 :
203 : // Compute the zeta divergence term. Since that also needs a transpose,
204 : // copy into data_buffer0, then transpose into `data_buffer1`.
205 : EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
206 : tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
207 : make_not_null(&data_buffer0), std::integral_constant<size_t, 2>{}));
208 : size_t chunk_size =
209 : xi_weak_div_matrix.rows() * eta_weak_div_matrix.rows();
210 : double* result_ptr = data_buffer1.data();
211 : raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
212 : data_buffer0.size() / chunk_size);
213 : apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
214 : zeta_weak_div_matrix, data_buffer0.size(),
215 : false);
216 : chunk_size =
217 : number_of_independent_components * zeta_weak_div_matrix.rows();
218 : result_ptr = divergence_of_fluxes->data();
219 : raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
220 : data_buffer1.size() / chunk_size);
221 :
222 : // Compute the eta divergence term. Since that also needs a transpose,
223 : // copy into data_buffer0, then transpose into `data_buffer1`.
224 : EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
225 : tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
226 : make_not_null(&data_buffer0), std::integral_constant<size_t, 1>{}));
227 : chunk_size = xi_weak_div_matrix.rows();
228 : result_ptr = data_buffer1.data();
229 : raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
230 : data_buffer1.size() / chunk_size);
231 : apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
232 : eta_weak_div_matrix, data_buffer0.size(),
233 : false);
234 : chunk_size = number_of_independent_components *
235 : eta_weak_div_matrix.rows() * zeta_weak_div_matrix.rows();
236 : result_ptr = data_buffer1.data();
237 : raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
238 : data_buffer0.size() / chunk_size);
239 : *divergence_of_fluxes += data_buffer1;
240 :
241 : // Now compute xi divergence and *add* to eta divergence
242 : EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
243 : tmpl::type_<FluxTags>{}, tmpl::type_<ResultTags>{},
244 : make_not_null(&data_buffer0), std::integral_constant<size_t, 0>{}));
245 : apply_matrix_in_first_dim(divergence_of_fluxes->data(),
246 : data_buffer0.data(), xi_weak_div_matrix,
247 : data_buffer0.size(), true);
248 : } else {
249 : static_assert(Dim == 1 or Dim == 2 or Dim == 3,
250 : "Weak divergence only implemented in 1d, 2d, and 3d.");
251 : }
252 : }
253 : }
254 :
255 : /// @{
256 : /*!
257 : * \brief Compute the weak divergence of fluxes in logical coordinates
258 : *
259 : * Applies the transpose logical differentiation matrix to the fluxes in each
260 : * dimension and sums over dimensions.
261 : *
262 : * \see weak_divergence
263 : */
264 : template <typename ResultTensor, typename FluxTensor, size_t Dim>
265 1 : void logical_weak_divergence(gsl::not_null<ResultTensor*> div_flux,
266 : const FluxTensor& flux, const Mesh<Dim>& mesh);
267 :
268 : /*!
269 : * \param divergence_of_fluxes Result buffer. Will be resized if necessary.
270 : * \param fluxes The fluxes to take the divergence of. Their first index must be
271 : * an upper spatial index in element-logical coordinates.
272 : * \param mesh The element mesh.
273 : * \param add_to_buffer If `true`, add the divergence to the existing values in
274 : * `divergence_of_fluxes`. Otherwise, overwrite the existing values.
275 : */
276 : template <typename... ResultTags, typename... FluxTags, size_t Dim>
277 1 : void logical_weak_divergence(
278 : const gsl::not_null<Variables<tmpl::list<ResultTags...>>*>
279 : divergence_of_fluxes,
280 : const Variables<tmpl::list<FluxTags...>>& fluxes, const Mesh<Dim>& mesh,
281 : const bool add_to_buffer = false) {
282 : static_assert(
283 : (... and
284 : std::is_same_v<tmpl::front<typename FluxTags::type::index_list>,
285 : SpatialIndex<Dim, UpLo::Up, Frame::ElementLogical>>),
286 : "The first index of each flux must be an upper spatial index in "
287 : "element-logical coordinates.");
288 : static_assert(
289 : (... and
290 : std::is_same_v<
291 : typename ResultTags::type,
292 : TensorMetafunctions::remove_first_index<typename FluxTags::type>>),
293 : "The result tensors must have the same type as the flux tensors with "
294 : "their first index removed.");
295 : if (not add_to_buffer) {
296 : divergence_of_fluxes->initialize(mesh.number_of_grid_points(), 0.);
297 : }
298 : EXPAND_PACK_LEFT_TO_RIGHT(logical_weak_divergence(
299 : make_not_null(&get<ResultTags>(*divergence_of_fluxes)),
300 : get<FluxTags>(fluxes), mesh));
301 : }
302 : /// @}
|