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