FirstOrderInternalPenalty.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <pup.h>
8 #include <string>
9 
11 #include "DataStructures/DataBox/Tag.hpp"
12 #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
16 #include "Domain/FaceNormal.hpp"
17 #include "Domain/Tags.hpp"
18 #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
19 #include "NumericalAlgorithms/DiscontinuousGalerkin/Protocols.hpp"
21 #include "Options/Options.hpp"
22 #include "Utilities/Gsl.hpp"
23 #include "Utilities/ProtocolHelpers.hpp"
24 #include "Utilities/TMPL.hpp"
25 
26 /// \cond
27 class DataVector;
28 namespace Tags {
29 template <typename>
30 struct Normalized;
31 } // namespace Tags
32 /// \endcond
33 
34 /// Numerical fluxes for elliptic systems
36 
37 /*!
38  * \brief The penalty factor in internal penalty fluxes
39  *
40  * The penalty factor is computed as
41  *
42  * \f{equation}
43  * \sigma = C \frac{N_\text{points}^2}{h}
44  * \f}
45  *
46  * where \f$N_\text{points} = 1 + N_p\f$ is the number of points (or one plus
47  * the polynomial degree) and \f$h\f$ is a measure of the element size. Both
48  * quantities are taken perpendicular to the face of the DG element that the
49  * penalty is being computed on. \f$C\f$ is the *penalty parameter*.
50  *
51  * \see `elliptic::dg::NumericalFluxes::FirstOrderInternalPenalty` for details
52  */
53 DataVector penalty(const DataVector& element_size, size_t num_points,
54  double penalty_parameter) noexcept;
55 
56 /*!
57  * \ingroup DiscontinuousGalerkinGroup
58  * \ingroup NumericalFluxesGroup
59  * \brief The internal penalty flux for first-order elliptic equations.
60  *
61  * \details Computes the internal penalty numerical flux (see e.g.
62  * \cite HesthavenWarburton, section 7.2) dotted with the interface unit normal.
63  *
64  * We implement here a suggested generalization of the internal penalty flux
65  * for any set of elliptic PDEs up to second order in derivatives. It is
66  * designed for fluxes (i.e. principal parts of the PDEs) that may depend on the
67  * dynamic fields, but do so at most linearly. This is the case for the velocity
68  * potential equation of binary neutron stars, for example. The linearity is
69  * only necessary to impose inhomogeneous boundary conditions as contributions
70  * to the fixed source (see
71  * `elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource`).
72  *
73  * We reduce the second-order elliptic PDEs to first-order form by introducing
74  * an _auxiliary_ variable \f$v\f$ for each _primal_ variable \f$u\f$ (i.e. for
75  * each variable whose second derivative appears in the equations). Then, the
76  * equations take the _flux-form_
77  *
78  * \f{align}
79  * -\partial_i F^i_A + S_A = f_A
80  * \f}
81  *
82  * where the fluxes \f$F^i_A(u,v)\f$ and sources \f$S_A(u,v)\f$ are indexed
83  * (with capital letters) by the variables and we have defined \f$f_A(x)\f$ as
84  * those sources that are independent of the variables. Note that the fluxes are
85  * functions of the primal and auxiliary variables; e.g. for a Poisson system
86  * \f$\{u,v_i\}\f$ on a spatial metric \f$\gamma_{ij}\f$ they are simply
87  * \f$F^i_u(v)=\sqrt{\gamma}\gamma^{ij} v_j\f$ and
88  * \f$F^i_{v_j}(u)=u\delta^i_j\f$ (see `Poisson::FirstOrderSystem`) or for an
89  * Elasticity system \f$\{\xi^i,S_{ij}\}\f$ with Young's tensor
90  * \f$Y^{ijkl}\f$ they are \f$F^i_{\xi^j}(S)=Y^{ijkl}S_{kl}\f$ and
91  * \f$F^i_{S_{jk}}(\xi)=\delta^i_{(j}\xi_{k)}\f$. Since each primal flux is
92  * commonly only a function of the corresponding auxiliary variable and
93  * vice-versa, we omit the other function arguments here to lighten the
94  * notational load.
95  *
96  * We now choose the internal penalty numerical fluxes \f$(n_i F^i_A)^*\f$ as
97  * follows for each primal variable \f$u\f$ and their corresponding auxiliary
98  * variable \f$v\f$:
99  *
100  * \f{align}
101  * (n_i F^i_u)^* &= \frac{1}{2} n_i \left( F^i_u(\partial_j
102  * F^j_v(u^\mathrm{int})) + F^i_u(\partial_j F^j_v(u^\mathrm{ext}))
103  * \right) - \sigma n_i \left(F^i_u(n_j F^j_v(u^\mathrm{int})) - F^i_u(
104  * n_j F^j_v(u^\mathrm{ext}))\right) \\
105  * (n_i F^i_v)^* &= \frac{1}{2} n_i \left(F^i_v(u^\mathrm{int}) +
106  * F^i_v(u^\mathrm{ext})\right)
107  * \f}
108  *
109  * Note that we have assumed \f$n^\mathrm{ext}_i=-n_i\f$ here, i.e. face normals
110  * don't depend on the dynamic variables (which may be discontinuous on element
111  * faces). This is the case for the problems we are expecting to solve, because
112  * those will be on fixed background metrics (e.g. a conformal metric for the
113  * XCTS system).
114  *
115  * Also note that the numerical fluxes intentionally don't depend on the
116  * auxiliary field values \f$v\f$. This property allows us to use the numerical
117  * fluxes also for the second-order (or _primal_) DG formulation, where we
118  * remove the need for an auxiliary variable. For the first-order system we
119  * could replace the divergence in \f$(n_i F^i_u)^*\f$ with \f$v\f$, which would
120  * result in a generalized "stabilized central flux" that is slightly less
121  * sparse than the internal penalty flux (see e.g. \cite HesthavenWarburton,
122  * section 7.2). We could also choose to ignore the fluxes in the penalty term,
123  * but preliminary tests suggest that this may hurt convergence.
124  *
125  * For a Poisson system (see above) this numeric flux reduces to the standard
126  * internal penalty flux (see e.g. \cite HesthavenWarburton, section 7.2, or
127  * \cite Arnold2002)
128  *
129  * \f{align}
130  * (n_i F^i_u)^* &= n_i v_i^* = \frac{1}{2} n_i \left(\partial_i u^\mathrm{int}
131  * + \partial_i u^\mathrm{ext}\right) - \sigma \left(u^\mathrm{int} -
132  * u^\mathrm{ext}\right) \\
133  * (n_i F^i_{v_j})^* &= n_j u^* = \frac{1}{2} n_j \left(u^\mathrm{int} +
134  * u^\mathrm{ext}\right)
135  * \f}
136  *
137  * where a sum over repeated indices is assumed, since the equation is
138  * formulated on a Euclidean geometry.
139  *
140  * This generalization of the internal penalty flux is based on unpublished work
141  * by Nils L. Fischer (nils.fischer@aei.mpg.de).
142  *
143  * The penalty factor \f$\sigma\f$ is responsible for removing zero eigenmodes
144  * and impacts the conditioning of the linear operator to be solved. It can be
145  * chosen as \f$\sigma=C\frac{N_\mathrm{points}^2}{h}\f$ where
146  * \f$N_\mathrm{points}\f$ is the number of collocation points (i.e. the
147  * polynomial degree plus 1), \f$h\f$ is a measure of the element size in
148  * inertial coordinates (both measured perpendicular to the interface under
149  * consideration) and \f$C\geq 1\f$ is a free parameter (see e.g.
150  * \cite HesthavenWarburton, section 7.2). For the element size we choose
151  * \f$h=\frac{J_\mathrm{volume}}{J_\mathrm{face}}\f$, i.e. the ratio of Jacobi
152  * determinants from logical to inertial coordinates in the element volume and
153  * on the element face, both evaluated on the face (see \cite Vincent2019qpd).
154  * Since both \f$N_\mathrm{points}\f$ and \f$h\f$ can be different on either
155  * side of the element boundary we take the maximum of \f$N_\mathrm{points}\f$
156  * and the pointwise minimum of \f$h\f$ across the element boundary as is done
157  * in \cite Vincent2019qpd. Note that we use the number of points
158  * \f$N_\mathrm{points}\f$ where \cite Vincent2019qpd uses the polynomial degree
159  * \f$N_\mathrm{points} - 1\f$ because we found unstable configurations on
160  * curved meshes when using the polynomial degree. Optimizing the penalty on
161  * curved meshes is subject to further investigation.
162  */
163 template <size_t Dim, typename FluxesComputerTag, typename FieldTagsList,
164  typename AuxiliaryFieldTagsList>
166 
167 template <size_t Dim, typename FluxesComputerTag, typename... FieldTags,
168  typename... AuxiliaryFieldTags>
169 struct FirstOrderInternalPenalty<Dim, FluxesComputerTag,
170  tmpl::list<FieldTags...>,
171  tmpl::list<AuxiliaryFieldTags...>>
172  : tt::ConformsTo<::dg::protocols::NumericalFlux> {
173  private:
174  using fluxes_computer_tag = FluxesComputerTag;
175  using FluxesComputer = typename fluxes_computer_tag::type;
176 
177  struct PerpendicularNumPoints : db::SimpleTag {
178  using type = size_t;
179  };
180 
181  struct ElementSize : db::SimpleTag {
182  using type = Scalar<DataVector>;
183  };
184 
185  template <typename Tag>
186  struct NormalDotDivFlux : db::PrefixTag, db::SimpleTag {
187  using tag = Tag;
189  };
190 
191  template <typename Tag>
192  struct NormalDotNormalDotFlux : db::PrefixTag, db::SimpleTag {
193  using tag = Tag;
195  };
196 
197  public:
198  struct PenaltyParameter {
199  using type = double;
200  static constexpr Options::String help = {
201  "The prefactor to the penalty term of the flux. Values closer to one "
202  "lead to better-conditioned problems, but on curved meshes the penalty "
203  "parameter may need to be increased to keep the problem well-defined."};
204  static type lower_bound() noexcept { return 1.; }
205  };
206  using options = tmpl::list<PenaltyParameter>;
207  static constexpr Options::String help = {
208  "The internal penalty flux for elliptic systems."};
209  static std::string name() noexcept { return "InternalPenalty"; }
210 
211  FirstOrderInternalPenalty() = default;
212  explicit FirstOrderInternalPenalty(const double penalty_parameter) noexcept
213  : penalty_parameter_(penalty_parameter) {}
214 
215  // clang-tidy: non-const reference
216  void pup(PUP::er& p) noexcept { p | penalty_parameter_; } // NOLINT
217 
218  using variables_tags = tmpl::list<FieldTags..., AuxiliaryFieldTags...>;
219 
220  using argument_tags = tmpl::push_front<
221  typename FluxesComputer::argument_tags, domain::Tags::Mesh<Dim>,
226  Frame::Inertial>>...,
228  fluxes_computer_tag>;
229  using volume_tags =
230  tmpl::push_front<get_volume_tags<FluxesComputer>, domain::Tags::Mesh<Dim>,
231  fluxes_computer_tag>;
232 
233  using package_field_tags =
234  tmpl::list<::Tags::NormalDotFlux<AuxiliaryFieldTags>...,
235  NormalDotDivFlux<AuxiliaryFieldTags>...,
236  NormalDotNormalDotFlux<AuxiliaryFieldTags>..., ElementSize>;
237  using package_extra_tags = tmpl::list<PerpendicularNumPoints>;
238 
239  template <typename... FluxesArgs>
240  void package_data(
241  const gsl::not_null<typename ::Tags::NormalDotFlux<
242  AuxiliaryFieldTags>::type*>... packaged_n_dot_aux_fluxes,
243  const gsl::not_null<typename NormalDotDivFlux<
244  AuxiliaryFieldTags>::type*>... packaged_n_dot_div_aux_fluxes,
245  const gsl::not_null<
246  typename NormalDotNormalDotFlux<AuxiliaryFieldTags>::
247  type*>... packaged_n_dot_principal_n_dot_aux_fluxes,
248  const gsl::not_null<Scalar<DataVector>*> element_size,
249  const gsl::not_null<size_t*> perpendicular_num_points,
250  const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
251  const Scalar<DataVector>& face_normal_magnitude,
252  const typename ::Tags::NormalDotFlux<
253  AuxiliaryFieldTags>::type&... normal_dot_auxiliary_field_fluxes,
254  const typename ::Tags::div<
255  ::Tags::Flux<AuxiliaryFieldTags, tmpl::size_t<Dim>,
256  Frame::Inertial>>::type&... div_auxiliary_field_fluxes,
257  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
258  const FluxesComputer& fluxes_computer,
259  const FluxesArgs&... fluxes_args) const noexcept {
260  *perpendicular_num_points = volume_mesh.extents(direction.dimension());
261  get(*element_size) = 2. / get(face_normal_magnitude);
262  auto principal_div_aux_field_fluxes = make_with_value<Variables<tmpl::list<
264  interface_unit_normal, std::numeric_limits<double>::signaling_NaN());
265  auto principal_ndot_aux_field_fluxes = make_with_value<Variables<tmpl::list<
267  interface_unit_normal, std::numeric_limits<double>::signaling_NaN());
268  fluxes_computer.apply(
270  &get<::Tags::Flux<FieldTags, tmpl::size_t<Dim>, Frame::Inertial>>(
271  principal_div_aux_field_fluxes))...,
272  fluxes_args..., div_auxiliary_field_fluxes...);
273  fluxes_computer.apply(
275  &get<::Tags::Flux<FieldTags, tmpl::size_t<Dim>, Frame::Inertial>>(
276  principal_ndot_aux_field_fluxes))...,
277  fluxes_args..., normal_dot_auxiliary_field_fluxes...);
278  const auto apply_normal_dot =
279  [&principal_div_aux_field_fluxes, &principal_ndot_aux_field_fluxes,
280  &interface_unit_normal](
281  auto field_tag_v, auto packaged_normal_dot_auxiliary_field_flux,
282  auto packaged_normal_dot_div_auxiliary_field_flux,
283  auto packaged_normal_dot_principal_ndot_aux_field_flux,
284  const auto& normal_dot_auxiliary_field_flux) noexcept {
285  using field_tag = decltype(field_tag_v);
286  // Compute n.F_v(u)
287  *packaged_normal_dot_auxiliary_field_flux =
288  normal_dot_auxiliary_field_flux;
289  // Compute n.F_u(div(F_v(u))) and n.F_u(n.F_v(u))
290  normal_dot_flux(
291  packaged_normal_dot_div_auxiliary_field_flux,
292  interface_unit_normal,
293  get<::Tags::Flux<field_tag, tmpl::size_t<Dim>, Frame::Inertial>>(
294  principal_div_aux_field_fluxes));
295  normal_dot_flux(
296  packaged_normal_dot_principal_ndot_aux_field_flux,
297  interface_unit_normal,
298  get<::Tags::Flux<field_tag, tmpl::size_t<Dim>, Frame::Inertial>>(
299  principal_ndot_aux_field_fluxes));
300  };
301  EXPAND_PACK_LEFT_TO_RIGHT(apply_normal_dot(
302  FieldTags{}, packaged_n_dot_aux_fluxes, packaged_n_dot_div_aux_fluxes,
303  packaged_n_dot_principal_n_dot_aux_fluxes,
304  normal_dot_auxiliary_field_fluxes));
305  }
306 
307  void operator()(
308  const gsl::not_null<typename ::Tags::NormalDotNumericalFlux<
309  FieldTags>::type*>... numerical_flux_for_fields,
310  const gsl::not_null<typename ::Tags::NormalDotNumericalFlux<
311  AuxiliaryFieldTags>::type*>... numerical_flux_for_auxiliary_fields,
312  const typename ::Tags::NormalDotFlux<
313  AuxiliaryFieldTags>::type&... normal_dot_auxiliary_flux_interiors,
314  const typename NormalDotDivFlux<
315  AuxiliaryFieldTags>::type&... normal_dot_div_auxiliary_flux_interiors,
316  const typename NormalDotNormalDotFlux<
317  AuxiliaryFieldTags>::type&... ndot_ndot_aux_flux_interiors,
318  const Scalar<DataVector>& element_size_interior,
319  const size_t perpendicular_num_points_interior,
320  const typename ::Tags::NormalDotFlux<AuxiliaryFieldTags>::
321  type&... minus_normal_dot_auxiliary_flux_exteriors,
322  const typename NormalDotDivFlux<
323  AuxiliaryFieldTags>::type&... minus_normal_dot_div_aux_flux_exteriors,
324  const typename NormalDotDivFlux<
325  AuxiliaryFieldTags>::type&... ndot_ndot_aux_flux_exteriors,
326  const Scalar<DataVector>& element_size_exterior,
327  const size_t perpendicular_num_points_exterior) const noexcept {
328  const auto penalty = ::elliptic::dg::NumericalFluxes::penalty(
329  min(get(element_size_interior), get(element_size_exterior)),
330  std::max(perpendicular_num_points_interior,
331  perpendicular_num_points_exterior),
332  penalty_parameter_);
333 
334  const auto assemble_numerical_fluxes = [&penalty](
335  const auto numerical_flux_for_field,
336  const auto numerical_flux_for_auxiliary_field,
337  const auto& normal_dot_auxiliary_flux_interior,
338  const auto& normal_dot_div_auxiliary_flux_interior,
339  const auto& ndot_ndot_aux_flux_interior,
340  const auto& minus_normal_dot_auxiliary_flux_exterior,
341  const auto& minus_normal_dot_div_aux_flux_exterior,
342  const auto& ndot_ndot_aux_flux_exterior) noexcept {
343  for (auto it = numerical_flux_for_auxiliary_field->begin();
344  it != numerical_flux_for_auxiliary_field->end(); it++) {
345  const auto index =
346  numerical_flux_for_auxiliary_field->get_tensor_index(it);
347  // We are working with the n.F_v(u) computed on either side of the
348  // interface, so (assuming the normal is independent of the dynamic
349  // variables) the data we get from the other element contains _minus_
350  // the normal from this element. So we cancel the minus sign when
351  // computing the average here.
352  *it = 0.5 * (normal_dot_auxiliary_flux_interior.get(index) -
353  minus_normal_dot_auxiliary_flux_exterior.get(index));
354  }
355  for (auto it = numerical_flux_for_field->begin();
356  it != numerical_flux_for_field->end(); it++) {
357  const auto index = numerical_flux_for_field->get_tensor_index(it);
358  *it = 0.5 * (normal_dot_div_auxiliary_flux_interior.get(index) -
359  minus_normal_dot_div_aux_flux_exterior.get(index)) -
360  penalty * (ndot_ndot_aux_flux_interior.get(index) -
361  ndot_ndot_aux_flux_exterior.get(index));
362  }
363  };
364  EXPAND_PACK_LEFT_TO_RIGHT(assemble_numerical_fluxes(
365  numerical_flux_for_fields, numerical_flux_for_auxiliary_fields,
366  normal_dot_auxiliary_flux_interiors,
367  normal_dot_div_auxiliary_flux_interiors, ndot_ndot_aux_flux_interiors,
368  minus_normal_dot_auxiliary_flux_exteriors,
369  minus_normal_dot_div_aux_flux_exteriors, ndot_ndot_aux_flux_exteriors));
370  }
371 
372  // This function computes the boundary contributions from Dirichlet boundary
373  // conditions. This data is what remains to be added to the boundaries when
374  // homogeneous (i.e. zero) boundary conditions are assumed in the calculation
375  // of the numerical fluxes, but we wish to impose inhomogeneous (i.e. nonzero)
376  // boundary conditions. Since this contribution does not depend on the
377  // numerical field values, but only on the Dirichlet boundary data, it may be
378  // added as contribution to the source of the elliptic systems. Then, it
379  // remains to solve the homogeneous problem with the modified source.
380  template <typename... FluxesArgs>
381  void compute_dirichlet_boundary(
382  const gsl::not_null<typename ::Tags::NormalDotNumericalFlux<
383  FieldTags>::type*>... numerical_flux_for_fields,
384  const gsl::not_null<typename ::Tags::NormalDotNumericalFlux<
385  AuxiliaryFieldTags>::type*>... numerical_flux_for_auxiliary_fields,
386  const typename FieldTags::type&... dirichlet_fields,
387  const Mesh<Dim>& volume_mesh, const Direction<Dim>& direction,
388  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
389  const Scalar<DataVector>& face_normal_magnitude,
390  const FluxesComputer& fluxes_computer,
391  const FluxesArgs&... fluxes_args) const noexcept {
392  const auto penalty = ::elliptic::dg::NumericalFluxes::penalty(
393  2. / get(face_normal_magnitude),
394  volume_mesh.extents(direction.dimension()), penalty_parameter_);
395 
396  // Compute n.F_v(u)
397  auto dirichlet_auxiliary_field_fluxes =
398  make_with_value<Variables<tmpl::list<::Tags::Flux<
399  AuxiliaryFieldTags, tmpl::size_t<Dim>, Frame::Inertial>...>>>(
400  interface_unit_normal,
402  fluxes_computer.apply(
403  make_not_null(&get<::Tags::Flux<AuxiliaryFieldTags, tmpl::size_t<Dim>,
404  Frame::Inertial>>(
405  dirichlet_auxiliary_field_fluxes))...,
406  fluxes_args..., dirichlet_fields...);
407  const auto apply_normal_dot_aux =
408  [&interface_unit_normal, &dirichlet_auxiliary_field_fluxes ](
409  auto auxiliary_field_tag_v,
410  const auto numerical_flux_for_auxiliary_field) noexcept {
411  using auxiliary_field_tag = decltype(auxiliary_field_tag_v);
412  normal_dot_flux(
413  numerical_flux_for_auxiliary_field, interface_unit_normal,
414  get<::Tags::Flux<auxiliary_field_tag, tmpl::size_t<Dim>,
415  Frame::Inertial>>(dirichlet_auxiliary_field_fluxes));
416  };
417  EXPAND_PACK_LEFT_TO_RIGHT(apply_normal_dot_aux(
418  AuxiliaryFieldTags{}, numerical_flux_for_auxiliary_fields));
419 
420  // Compute 2 * sigma * n.F_u(n.F_v(u))
421  auto principal_dirichlet_auxiliary_field_fluxes =
422  make_with_value<Variables<tmpl::list<
424  interface_unit_normal,
426  fluxes_computer.apply(
428  &get<::Tags::Flux<FieldTags, tmpl::size_t<Dim>, Frame::Inertial>>(
429  principal_dirichlet_auxiliary_field_fluxes))...,
430  fluxes_args..., *numerical_flux_for_auxiliary_fields...);
431  const auto assemble_dirichlet_penalty = [
432  &interface_unit_normal, &penalty,
433  &principal_dirichlet_auxiliary_field_fluxes
434  ](auto field_tag_v, const auto numerical_flux_for_field) noexcept {
435  using field_tag = decltype(field_tag_v);
436  normal_dot_flux(
437  numerical_flux_for_field, interface_unit_normal,
438  get<::Tags::Flux<field_tag, tmpl::size_t<Dim>, Frame::Inertial>>(
439  principal_dirichlet_auxiliary_field_fluxes));
440  for (auto it = numerical_flux_for_field->begin();
441  it != numerical_flux_for_field->end(); it++) {
442  *it *= 2 * penalty;
443  }
444  };
446  assemble_dirichlet_penalty(FieldTags{}, numerical_flux_for_fields));
447  }
448 
449  private:
450  double penalty_parameter_ = std::numeric_limits<double>::signaling_NaN();
451 };
452 
453 } // namespace elliptic::dg::NumericalFluxes
domain::push_front
CoordinateMap< SourceFrame, TargetFrame, NewMap, Maps... > push_front(CoordinateMap< SourceFrame, TargetFrame, Maps... > old_map, NewMap new_map) noexcept
Creates a CoordinateMap by prepending the new map to the beginning of the old maps.
std::string
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:563
Frame::Inertial
Definition: IndexType.hpp:44
Divergence.hpp
get
constexpr Tag::type & get(Variables< TagList > &v) noexcept
Return Tag::type pointing into the contiguous array.
Definition: Variables.hpp:638
db::PrefixTag
Mark a struct as a prefix tag by inheriting from this.
Definition: Tag.hpp:103
elliptic::dg::NumericalFluxes
Numerical fluxes for elliptic systems.
Definition: FirstOrderInternalPenalty.cpp:11
Options.hpp
FaceNormal.hpp
Tags.hpp
domain::Tags::Mesh
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:107
db::SimpleTag
Mark a struct as a simple tag by inheriting from this.
Definition: Tag.hpp:36
Tags::NormalDotFlux
Prefix indicating a boundary unit normal vector dotted into the flux.
Definition: Prefixes.hpp:96
Tags::Flux
Prefix indicating a flux.
Definition: Prefixes.hpp:40
elliptic::dg::NumericalFluxes::penalty
DataVector penalty(const DataVector &element_size, const size_t num_points, const double penalty_parameter) noexcept
The penalty factor in internal penalty fluxes.
Definition: FirstOrderInternalPenalty.cpp:13
Direction
Definition: Direction.hpp:23
Tags::div
Prefix indicating the divergence.
Definition: Divergence.hpp:44
Tags::Normalized
Definition: Magnitude.hpp:137
make_with_value
std::remove_const_t< R > make_with_value(const T &input, const ValueType &value) noexcept
Given an object of type T, create an object of type R whose elements are initialized to value.
Definition: MakeWithValue.hpp:88
cstddef
TensorMetafunctions::remove_first_index
::Tensor< typename Tensor::type, tmpl::pop_front< typename Tensor::symmetry >, tmpl::pop_front< typename Tensor::index_list > > remove_first_index
remove the first index of a tensor
Definition: Metafunctions.hpp:113
DataVector
Stores a collection of function values.
Definition: DataVector.hpp:42
std::numeric_limits::signaling_NaN
T signaling_NaN(T... args)
Metafunctions.hpp
elliptic::dg::NumericalFluxes::FirstOrderInternalPenalty
The internal penalty flux for first-order elliptic equations.
Definition: FirstOrderInternalPenalty.hpp:165
std::min
T min(T... args)
Variables.hpp
Mesh
Holds the number of grid points, basis, and quadrature in each direction of the computational grid.
Definition: Mesh.hpp:47
Scalar
Tensor< T, Symmetry<>, index_list<> > Scalar
Definition: TypeAliases.hpp:21
Tags::Magnitude
Definition: Magnitude.hpp:98
Gsl.hpp
Options::String
const char *const String
The string used in option structs.
Definition: Options.hpp:32
domain::Tags::Direction
Definition: Tags.hpp:381
Tensor.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
Prefixes.hpp
tt::ConformsTo
Indicate a class conforms to the Protocol.
Definition: ProtocolHelpers.hpp:22
TMPL.hpp
gsl::not_null
Require a pointer to not be a nullptr
Definition: Gsl.hpp:183
string