WeakDivergence.hpp
2 // See LICENSE.txt for details.
3
4 #pragma once
5
6 #include <cstddef>
7 #include <type_traits>
8
13
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  * \int_{\Omega}d^n x\,F^i_{\breve{\jmath}} \phi_{\breve{\jmath}} \partial_i
56  * \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... FluxTags, size_t Dim>
66  const gsl::not_null<Variables<tmpl::list<Tags::div<FluxTags>...>>*>
67  divergence_of_fluxes,
68  const Variables<tmpl::list<FluxTags...>>& fluxes, const Mesh<Dim>& mesh,
69  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
70  det_jac_times_inverse_jacobian) noexcept {
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) noexcept {
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
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(),
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,
108  result_buffer,
109  auto logical_index_of_jacobian) noexcept {
110  using flux_tag = tmpl::type_from<decltype(flux_tag_v)>;
111  using div_tag = Tags::div<flux_tag>;
112
113  auto& result = get<div_tag>(*result_buffer);
114  const auto& flux = get<flux_tag>(fluxes);
115
116  for (size_t result_storage_index = 0;
117  result_storage_index < result.size(); ++result_storage_index) {
118  const auto result_tensor_index =
119  result.get_tensor_index(result_storage_index);
120  const auto flux_x_tensor_index = prepend(result_tensor_index, 0_st);
121  const auto flux_y_tensor_index = prepend(result_tensor_index, 1_st);
122  if constexpr (Dim == 2) {
123  result[result_storage_index] =
124  get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
125  0>(det_jac_times_inverse_jacobian) *
126  flux.get(flux_x_tensor_index) +
127  get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
128  1>(det_jac_times_inverse_jacobian) *
129  flux.get(flux_y_tensor_index);
130  } else {
131  const auto flux_z_tensor_index =
132  prepend(result_tensor_index, 2_st);
133  result[result_storage_index] =
134  get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
135  0>(det_jac_times_inverse_jacobian) *
136  flux.get(flux_x_tensor_index) +
137  get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
138  1>(det_jac_times_inverse_jacobian) *
139  flux.get(flux_y_tensor_index) +
140  get<std::decay_t<decltype(logical_index_of_jacobian)>::value,
141  2>(det_jac_times_inverse_jacobian) *
142  flux.get(flux_z_tensor_index);
143  }
144  }
145  };
146
147  if constexpr (Dim == 2) {
148  Variables<tmpl::list<Tags::div<FluxTags>...>> data_buffer{
149  divergence_of_fluxes->number_of_grid_points()};
150
151  const Matrix& eta_weak_div_matrix =
152  Spectral::weak_flux_differentiation_matrix(mesh.slice_through(1));
153  const Matrix& xi_weak_div_matrix =
154  Spectral::weak_flux_differentiation_matrix(mesh.slice_through(0));
155
156  // Compute the eta divergence term. Since that also needs a transpose,
157  // copy into result, then transpose into data_buffer
158  EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
159  tmpl::type_<FluxTags>{}, make_not_null(&data_buffer),
161  double* div_ptr = divergence_of_fluxes->data();
162  raw_transpose(make_not_null(div_ptr), data_buffer.data(),
163  xi_weak_div_matrix.rows(),
164  divergence_of_fluxes->size() / xi_weak_div_matrix.rows());
165  apply_matrix_in_first_dim(data_buffer.data(),
166  divergence_of_fluxes->data(),
167  eta_weak_div_matrix, data_buffer.size(), false);
168
169  const size_t chunk_size = Variables<tmpl::list<Tags::div<FluxTags>...>>::
170  number_of_independent_components *
171  eta_weak_div_matrix.rows();
172  raw_transpose(make_not_null(div_ptr), data_buffer.data(), chunk_size,
173  data_buffer.size() / chunk_size);
174
175  // Now compute xi divergence and *add* to eta divergence
176  EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
177  tmpl::type_<FluxTags>{}, make_not_null(&data_buffer),
179  apply_matrix_in_first_dim(divergence_of_fluxes->data(),
180  data_buffer.data(), xi_weak_div_matrix,
181  data_buffer.size(), true);
182  } else if constexpr (Dim == 3) {
183  Variables<tmpl::list<Tags::div<FluxTags>...>> data_buffer0{
184  divergence_of_fluxes->number_of_grid_points()};
185  Variables<tmpl::list<Tags::div<FluxTags>...>> data_buffer1{
186  divergence_of_fluxes->number_of_grid_points()};
187  constexpr size_t number_of_independent_components =
188  decltype(data_buffer1)::number_of_independent_components;
189
190  const Matrix& zeta_weak_div_matrix =
191  Spectral::weak_flux_differentiation_matrix(mesh.slice_through(2));
192  const Matrix& eta_weak_div_matrix =
193  Spectral::weak_flux_differentiation_matrix(mesh.slice_through(1));
194  const Matrix& xi_weak_div_matrix =
195  Spectral::weak_flux_differentiation_matrix(mesh.slice_through(0));
196
197  // Compute the zeta divergence term. Since that also needs a transpose,
198  // copy into data_buffer0, then transpose into data_buffer1.
199  EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
200  tmpl::type_<FluxTags>{}, make_not_null(&data_buffer0),
202  size_t chunk_size =
203  xi_weak_div_matrix.rows() * eta_weak_div_matrix.rows();
204  double* result_ptr = data_buffer1.data();
205  raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
206  data_buffer0.size() / chunk_size);
207  apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
208  zeta_weak_div_matrix, data_buffer0.size(),
209  false);
210  chunk_size =
211  number_of_independent_components * zeta_weak_div_matrix.rows();
212  result_ptr = divergence_of_fluxes->data();
213  raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
214  data_buffer1.size() / chunk_size);
215
216  // Compute the eta divergence term. Since that also needs a transpose,
217  // copy into data_buffer0, then transpose into data_buffer1.
218  EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
219  tmpl::type_<FluxTags>{}, make_not_null(&data_buffer0),
221  chunk_size = xi_weak_div_matrix.rows();
222  result_ptr = data_buffer1.data();
223  raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
224  data_buffer1.size() / chunk_size);
225  apply_matrix_in_first_dim(data_buffer0.data(), data_buffer1.data(),
226  eta_weak_div_matrix, data_buffer0.size(),
227  false);
228  chunk_size = number_of_independent_components *
229  eta_weak_div_matrix.rows() * zeta_weak_div_matrix.rows();
230  result_ptr = data_buffer1.data();
231  raw_transpose(make_not_null(result_ptr), data_buffer0.data(), chunk_size,
232  data_buffer0.size() / chunk_size);
233  *divergence_of_fluxes += data_buffer1;
234
235  // Now compute xi divergence and *add* to eta divergence
236  EXPAND_PACK_LEFT_TO_RIGHT(transform_to_logical_frame(
237  tmpl::type_<FluxTags>{}, make_not_null(&data_buffer0),
239  apply_matrix_in_first_dim(divergence_of_fluxes->data(),
240  data_buffer0.data(), xi_weak_div_matrix,
241  data_buffer0.size(), true);
242  } else {
243  static_assert(Dim == 1 or Dim == 2 or Dim == 3,
244  "Weak divergence only implemented in 1d, 2d, and 3d.");
245  }
246  }
247 }
Matrix
A dynamically sized matrix of doubles with column-major storage.
Definition: Matrix.hpp:19
std::integral_constant
EXPAND_PACK_LEFT_TO_RIGHT
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:601
UNLIKELY
#define UNLIKELY(x)
Definition: Gsl.hpp:73
Divergence.hpp
Literals.hpp
Spectral::weak_flux_differentiation_matrix
const Matrix & weak_flux_differentiation_matrix(size_t num_points) noexcept
Matrix used to compute the divergence of the flux in weak form.
Transpose.hpp
db::get
const auto & get(const DataBox< TagList > &box) noexcept
Retrieve the item with tag Tag from the DataBox.
Definition: DataBox.hpp:791
prepend
constexpr std::array< T, Dim+1 > prepend(const std::array< T, Dim > &a, T value) noexcept
Construct an array from an existing array prepending a value.
Definition: StdArrayHelpers.hpp:132
Spectral.hpp
dgemm_< true >
void dgemm_< true >(const char &TRANSA, const char &TRANSB, const size_t &M, const size_t &N, const size_t &K, const double &ALPHA, const double *A, const size_t &LDA, const double *B, const size_t &LDB, const double &BETA, double *C, const size_t &LDC)
Perform a matrix-matrix multiplication.
Definition: Blas.hpp:123
Tags::div
Prefix indicating the divergence.
Definition: Divergence.hpp:44
weak_divergence
void weak_divergence(const gsl::not_null< Variables< tmpl::list< Tags::div< FluxTags >... >> * > divergence_of_fluxes, const Variables< tmpl::list< FluxTags... >> &fluxes, const Mesh< Dim > &mesh, const InverseJacobian< DataVector, Dim, Frame::Logical, Frame::Inertial > &det_jac_times_inverse_jacobian) noexcept
Compute the weak form divergence of fluxes.
Definition: WeakDivergence.hpp:65
cstddef
std::decay_t
raw_transpose
void raw_transpose(const gsl::not_null< T * > result, const T *const data, const size_t chunk_size, const size_t number_of_chunks) noexcept
Function to compute transposed data.
Definition: Transpose.hpp:25
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:49
Gsl.hpp
Tensor.hpp
Matrix.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
Blas.hpp
type_traits
TMPL.hpp
Mesh.hpp
gsl::not_null
Require a pointer to not be a nullptr