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 
12 #include "DataStructures/DataBox/Tag.hpp"
16 #include "Domain/FaceNormal.hpp"
17 #include "Domain/Tags.hpp"
18 #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
20 #include "Options/Options.hpp"
21 #include "Utilities/Gsl.hpp"
22 #include "Utilities/TMPL.hpp"
23 
24 /// \cond
25 class DataVector;
26 namespace Tags {
27 template <typename>
28 struct Normalized;
29 } // namespace Tags
30 /// \endcond
31 
32 namespace elliptic {
33 namespace dg {
34 namespace NumericalFluxes {
35 
36 /*!
37  * \ingroup DiscontinuousGalerkinGroup
38  * \ingroup NumericalFluxesGroup
39  * \brief The internal penalty flux for first-order elliptic equations.
40  *
41  * \details Computes the internal penalty numerical flux (see e.g.
42  * \cite HesthavenWarburton, section 7.2) dotted with the interface unit normal.
43  *
44  * We implement here a suggested generalization of the internal penalty flux
45  * for any set of elliptic PDEs up to second order in derivatives. It is
46  * designed for fluxes (i.e. principal parts of the PDEs) that may depend on the
47  * dynamic fields, but do so at most linearly. This is the case for the velocity
48  * potential equation of binary neutron stars, for example. The linearity is
49  * only necessary to impose inhomogeneous boundary conditions as contributions
50  * to the fixed source (see
51  * `elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource`).
52  *
53  * We reduce the second-order elliptic PDEs to first-order form by introducing
54  * an _auxiliary_ variable \f$v\f$ for each _primal_ variable \f$u\f$ (i.e. for
55  * each variable whose second derivative appears in the equations). Then, the
56  * equations take the _flux-form_
57  *
58  * \f{align}
59  * -\partial_i F^i_A + S_A = f_A
60  * \f}
61  *
62  * where the fluxes \f$F^i_A(u,v)\f$ and sources \f$S_A(u,v)\f$ are indexed
63  * (with capital letters) by the variables and we have defined \f$f_A(x)\f$ as
64  * those sources that are independent of the variables. Note that the fluxes are
65  * functions of the primal and auxiliary variables; e.g. for a Poisson system
66  * \f$\{u,v_i\}\f$ on a spatial metric \f$\gamma_{ij}\f$ they are simply
67  * \f$F^i_u(v)=\sqrt{\gamma}\gamma^{ij} v_j\f$ and
68  * \f$F^i_{v_j}(u)=u\delta^i_j\f$ (see `Poisson::FirstOrderSystem`) or for an
69  * Elasticity system \f$\{\xi^i,S_{ij}\}\f$ with Young's tensor
70  * \f$Y^{ijkl}\f$ they are \f$F^i_{\xi^j}(S)=Y^{ijkl}S_{kl}\f$ and
71  * \f$F^i_{S_{jk}}(\xi)=\delta^i_{(j}\xi_{k)}\f$. Since each primal flux is
72  * commonly only a function of the corresponding auxiliary variable and
73  * vice-versa, we omit the other function arguments here to lighten the
74  * notational load.
75  *
76  * We now choose the internal penalty numerical fluxes \f$(n_i F^i_A)^*\f$ as
77  * follows for each primal variable \f$u\f$ and their corresponding auxiliary
78  * variable \f$v\f$:
79  *
80  * \f{align}
81  * (n_i F^i_u)^* &= \frac{1}{2} n_i \left( F^i_u(\partial_j
82  * F^j_v(u^\mathrm{int})) + F^i_u(\partial_j F^j_v(u^\mathrm{ext}))
83  * \right) - \sigma n_i \left(F^i_u(n_j F^j_v(u^\mathrm{int})) - F^i_u(
84  * n_j F^j_v(u^\mathrm{ext}))\right) \\
85  * (n_i F^i_v)^* &= \frac{1}{2} n_i \left(F^i_v(u^\mathrm{int}) +
86  * F^i_v(u^\mathrm{ext})\right)
87  * \f}
88  *
89  * Note that we have assumed \f$n^\mathrm{ext}_i=-n_i\f$ here, i.e. face normals
90  * don't depend on the dynamic variables (which may be discontinuous on element
91  * faces). This is the case for the problems we are expecting to solve, because
92  * those will be on fixed background metrics (e.g. a conformal metric for the
93  * XCTS system).
94  *
95  * Also note that the numerical fluxes intentionally don't depend on the
96  * auxiliary field values \f$v\f$. This property allows us to use the numerical
97  * fluxes also for the second-order (or _primal_) DG formulation, where we
98  * remove the need for an auxiliary variable. For the first-order system we
99  * could replace the divergence in \f$(n_i F^i_u)^*\f$ with \f$v\f$, which would
100  * result in a generalized "stabilized central flux" that is slightly less
101  * sparse than the internal penalty flux (see e.g. \cite HesthavenWarburton,
102  * section 7.2). We could also choose to ignore the fluxes in the penalty term,
103  * but preliminary tests suggest that this may hurt convergence.
104  *
105  * For a Poisson system (see above) this numeric flux reduces to the standard
106  * internal penalty flux (see e.g. \cite HesthavenWarburton, section 7.2, or
107  * \cite Arnold2002)
108  *
109  * \f{align}
110  * (n_i F^i_u)^* &= n_i v_i^* = \frac{1}{2} n_i \left(\partial_i u^\mathrm{int}
111  * + \partial_i u^\mathrm{ext}\right) - \sigma \left(u^\mathrm{int} -
112  * u^\mathrm{ext}\right) \\
113  * (n_i F^i_{v_j})^* &= n_j u^* = \frac{1}{2} n_j \left(u^\mathrm{int} +
114  * u^\mathrm{ext}\right)
115  * \f}
116  *
117  * where a sum over repeated indices is assumed, since the equation is
118  * formulated on a Euclidean geometry.
119  *
120  * This generalization of the internal penalty flux is based on unpublished work
121  * by Nils L. Fischer (nils.fischer@aei.mpg.de).
122  *
123  * The penalty factor \f$\sigma\f$ is responsible for removing zero eigenmodes
124  * and impacts the conditioning of the linear operator to be solved. It can be
125  * chosen as \f$\sigma=C\frac{N_\mathrm{points}^2}{h}\f$ where
126  * \f$N_\mathrm{points}\f$ is the number of collocation points (i.e. the
127  * polynomial degree plus 1), \f$h\f$ is a measure of the element size in
128  * inertial coordinates (both measured perpendicular to the interface under
129  * consideration) and \f$C\geq 1\f$ is a free parameter (see e.g.
130  * \cite HesthavenWarburton, section 7.2).
131  */
132 template <size_t Dim, typename FluxesComputerTag, typename FieldTagsList,
133  typename AuxiliaryFieldTagsList>
135 
136 template <size_t Dim, typename FluxesComputerTag, typename... FieldTags,
137  typename... AuxiliaryFieldTags>
138 struct FirstOrderInternalPenalty<Dim, FluxesComputerTag,
139  tmpl::list<FieldTags...>,
140  tmpl::list<AuxiliaryFieldTags...>> {
141  private:
142  using fluxes_computer_tag = FluxesComputerTag;
144 
145  template <typename Tag>
146  struct NormalDotDivFlux : db::PrefixTag, db::SimpleTag {
147  using tag = Tag;
149  };
150 
151  template <typename Tag>
152  struct NormalDotNormalDotFlux : db::PrefixTag, db::SimpleTag {
153  using tag = Tag;
155  };
156 
157  public:
158  struct PenaltyParameter {
159  using type = double;
160  // Currently this is used as the full prefactor `sigma` to the penalty term.
161  // This means it needs to be chosen large enough so that the scheme is
162  // stable everywhere. A good estimate is the the largest
163  // `sigma > N_points^2 / h` (see class documentation) over all elements in
164  // the domain. Choosing `sigma` the same everywhere means it is an
165  // overestimate on non-uniform meshes where elements are large or polynomial
166  // degrees are small, but this only affects the conditioning of the scheme,
167  // i.e. it will converge slower but to the same solution. When it becomes
168  // possible to communicate non-tensors (the element size and the number of
169  // collocation points on either side of the mortar), this option will be
170  // changed to be just the free parameter `C` multiplying `N_points^2 / h`.
171  // This will improve conditioning on non-uniform meshes. Currently, the
172  // `packaged_data` is a `Variables` which can only hold tensors.
173  static constexpr OptionString help = {
174  "The prefactor to the penalty term of the flux."};
175  };
176  using options = tmpl::list<PenaltyParameter>;
177  static constexpr OptionString help = {
178  "The internal penalty flux for elliptic systems."};
179  static std::string name() noexcept { return "InternalPenalty"; }
180 
181  FirstOrderInternalPenalty() = default;
182  explicit FirstOrderInternalPenalty(const double penalty_parameter) noexcept
183  : penalty_parameter_(penalty_parameter) {}
184 
185  // clang-tidy: non-const reference
186  void pup(PUP::er& p) noexcept { p | penalty_parameter_; } // NOLINT
187 
188  // These tags are sliced to the interface of the element and passed to
189  // `package_data` to provide the data needed to compute the numerical fluxes.
190  using argument_tags = tmpl::push_front<
191  typename FluxesComputer::argument_tags,
194  Frame::Inertial>>...,
196  fluxes_computer_tag>;
197  using volume_tags =
198  tmpl::push_front<get_volume_tags<FluxesComputer>, fluxes_computer_tag>;
199 
200  // This is the data needed to compute the numerical flux.
201  // `SendBoundaryFluxes` calls `package_data` to store these tags in a
202  // Variables. Local and remote values of this data are then combined in the
203  // `()` operator.
204  using package_tags =
205  tmpl::list<::Tags::NormalDotFlux<AuxiliaryFieldTags>...,
206  NormalDotDivFlux<AuxiliaryFieldTags>...,
207  NormalDotNormalDotFlux<AuxiliaryFieldTags>...>;
208 
209  // Following the packaged_data pointer, this function expects as arguments the
210  // types in `argument_tags`.
211  template <typename... FluxesArgs>
212  void package_data(
213  const gsl::not_null<Variables<package_tags>*> packaged_data,
215  AuxiliaryFieldTags>>&... normal_dot_auxiliary_field_fluxes,
217  ::Tags::Flux<AuxiliaryFieldTags, tmpl::size_t<Dim>,
218  Frame::Inertial>>>&... div_auxiliary_field_fluxes,
219  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
220  const FluxesComputer& fluxes_computer,
221  const FluxesArgs&... fluxes_args) const noexcept {
222  auto principal_div_aux_field_fluxes = make_with_value<Variables<tmpl::list<
224  interface_unit_normal, std::numeric_limits<double>::signaling_NaN());
225  auto principal_ndot_aux_field_fluxes = make_with_value<Variables<tmpl::list<
226  ::Tags::Flux<FieldTags, tmpl::size_t<Dim>, Frame::Inertial>...>>>(
227  interface_unit_normal, std::numeric_limits<double>::signaling_NaN());
228  fluxes_computer.apply(
230  &get<::Tags::Flux<FieldTags, tmpl::size_t<Dim>, Frame::Inertial>>(
231  principal_div_aux_field_fluxes))...,
232  fluxes_args..., div_auxiliary_field_fluxes...);
233  fluxes_computer.apply(
235  &get<::Tags::Flux<FieldTags, tmpl::size_t<Dim>, Frame::Inertial>>(
236  principal_ndot_aux_field_fluxes))...,
237  fluxes_args..., normal_dot_auxiliary_field_fluxes...);
238  const auto apply_normal_dot =
239  [
240  &packaged_data, &principal_div_aux_field_fluxes,
241  &principal_ndot_aux_field_fluxes, &interface_unit_normal
242  ](auto field_tag_v, auto auxiliary_field_tag_v,
243  const auto& normal_dot_auxiliary_field_flux) noexcept {
244  using field_tag = decltype(field_tag_v);
245  using auxiliary_field_tag = decltype(auxiliary_field_tag_v);
246  // Compute n.F_v(u)
247  get<::Tags::NormalDotFlux<auxiliary_field_tag>>(*packaged_data) =
248  normal_dot_auxiliary_field_flux;
249  // Compute n.F_u(div(F_v(u))) and n.F_u(n.F_v(u))
250  normal_dot_flux(
252  &get<NormalDotDivFlux<auxiliary_field_tag>>(*packaged_data)),
253  interface_unit_normal,
254  get<::Tags::Flux<field_tag, tmpl::size_t<Dim>, Frame::Inertial>>(
255  principal_div_aux_field_fluxes));
256  normal_dot_flux(
257  make_not_null(&get<NormalDotNormalDotFlux<auxiliary_field_tag>>(
258  *packaged_data)),
259  interface_unit_normal,
260  get<::Tags::Flux<field_tag, tmpl::size_t<Dim>, Frame::Inertial>>(
261  principal_ndot_aux_field_fluxes));
262  };
263  EXPAND_PACK_LEFT_TO_RIGHT(apply_normal_dot(
264  FieldTags{}, AuxiliaryFieldTags{}, normal_dot_auxiliary_field_fluxes));
265  }
266 
267  // This function combines local and remote data to the numerical fluxes.
268  // The numerical fluxes as not-null pointers are the first arguments. The
269  // other arguments are the packaged types for the interior side followed by
270  // the packaged types for the exterior side.
271  void operator()(
273  FieldTags>>*>... numerical_flux_for_fields,
275  AuxiliaryFieldTags>>*>... numerical_flux_for_auxiliary_fields,
277  AuxiliaryFieldTags>>&... normal_dot_auxiliary_flux_interiors,
278  const db::item_type<NormalDotDivFlux<
279  AuxiliaryFieldTags>>&... normal_dot_div_auxiliary_flux_interiors,
280  const db::item_type<NormalDotNormalDotFlux<
281  AuxiliaryFieldTags>>&... ndot_ndot_aux_flux_interiors,
283  AuxiliaryFieldTags>>&... minus_normal_dot_auxiliary_flux_exteriors,
284  const db::item_type<NormalDotDivFlux<
285  AuxiliaryFieldTags>>&... minus_normal_dot_div_aux_flux_exteriors,
286  const db::item_type<NormalDotDivFlux<
287  AuxiliaryFieldTags>>&... ndot_ndot_aux_flux_exteriors) const
288  noexcept {
289  // Need polynomial degrees and element size to compute this dynamically
290  const double penalty = penalty_parameter_;
291 
292  const auto assemble_numerical_fluxes = [&penalty](
293  const auto numerical_flux_for_field,
294  const auto numerical_flux_for_auxiliary_field,
295  const auto& normal_dot_auxiliary_flux_interior,
296  const auto& normal_dot_div_auxiliary_flux_interior,
297  const auto& ndot_ndot_aux_flux_interior,
298  const auto& minus_normal_dot_auxiliary_flux_exterior,
299  const auto& minus_normal_dot_div_aux_flux_exterior,
300  const auto& ndot_ndot_aux_flux_exterior) noexcept {
301  for (auto it = numerical_flux_for_auxiliary_field->begin();
302  it != numerical_flux_for_auxiliary_field->end(); it++) {
303  const auto index =
304  numerical_flux_for_auxiliary_field->get_tensor_index(it);
305  // We are working with the n.F_v(u) computed on either side of the
306  // interface, so (assuming the normal is independent of the dynamic
307  // variables) the data we get from the other element contains _minus_
308  // the normal from this element. So we cancel the minus sign when
309  // computing the average here.
310  *it = 0.5 * (normal_dot_auxiliary_flux_interior.get(index) -
311  minus_normal_dot_auxiliary_flux_exterior.get(index));
312  }
313  for (auto it = numerical_flux_for_field->begin();
314  it != numerical_flux_for_field->end(); it++) {
315  const auto index = numerical_flux_for_field->get_tensor_index(it);
316  *it = 0.5 * (normal_dot_div_auxiliary_flux_interior.get(index) -
317  minus_normal_dot_div_aux_flux_exterior.get(index)) -
318  penalty * (ndot_ndot_aux_flux_interior.get(index) -
319  ndot_ndot_aux_flux_exterior.get(index));
320  }
321  };
322  EXPAND_PACK_LEFT_TO_RIGHT(assemble_numerical_fluxes(
323  numerical_flux_for_fields, numerical_flux_for_auxiliary_fields,
324  normal_dot_auxiliary_flux_interiors,
325  normal_dot_div_auxiliary_flux_interiors, ndot_ndot_aux_flux_interiors,
326  minus_normal_dot_auxiliary_flux_exteriors,
327  minus_normal_dot_div_aux_flux_exteriors, ndot_ndot_aux_flux_exteriors));
328  }
329 
330  // This function computes the boundary contributions from Dirichlet boundary
331  // conditions. This data is what remains to be added to the boundaries when
332  // homogeneous (i.e. zero) boundary conditions are assumed in the calculation
333  // of the numerical fluxes, but we wish to impose inhomogeneous (i.e. nonzero)
334  // boundary conditions. Since this contribution does not depend on the
335  // numerical field values, but only on the Dirichlet boundary data, it may be
336  // added as contribution to the source of the elliptic systems. Then, it
337  // remains to solve the homogeneous problem with the modified source.
338  template <typename... FluxesArgs>
339  void compute_dirichlet_boundary(
341  FieldTags>>*>... numerical_flux_for_fields,
343  AuxiliaryFieldTags>>*>... numerical_flux_for_auxiliary_fields,
344  const db::const_item_type<FieldTags>&... dirichlet_fields,
345  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal,
346  const FluxesComputer& fluxes_computer,
347  const FluxesArgs&... fluxes_args) const noexcept {
348  // Need polynomial degrees and element size to compute this dynamically
349  const double penalty = penalty_parameter_;
350 
351  // Compute n.F_v(u)
352  auto dirichlet_auxiliary_field_fluxes =
353  make_with_value<Variables<tmpl::list<::Tags::Flux<
354  AuxiliaryFieldTags, tmpl::size_t<Dim>, Frame::Inertial>...>>>(
355  interface_unit_normal,
357  fluxes_computer.apply(
358  make_not_null(&get<::Tags::Flux<AuxiliaryFieldTags, tmpl::size_t<Dim>,
359  Frame::Inertial>>(
360  dirichlet_auxiliary_field_fluxes))...,
361  fluxes_args..., dirichlet_fields...);
362  const auto apply_normal_dot_aux =
363  [&interface_unit_normal, &dirichlet_auxiliary_field_fluxes ](
364  auto auxiliary_field_tag_v,
365  const auto numerical_flux_for_auxiliary_field) noexcept {
366  using auxiliary_field_tag = decltype(auxiliary_field_tag_v);
367  normal_dot_flux(
368  numerical_flux_for_auxiliary_field, interface_unit_normal,
369  get<::Tags::Flux<auxiliary_field_tag, tmpl::size_t<Dim>,
370  Frame::Inertial>>(dirichlet_auxiliary_field_fluxes));
371  };
372  EXPAND_PACK_LEFT_TO_RIGHT(apply_normal_dot_aux(
373  AuxiliaryFieldTags{}, numerical_flux_for_auxiliary_fields));
374 
375  // Compute 2 * sigma * n.F_u(n.F_v(u))
376  auto principal_dirichlet_auxiliary_field_fluxes =
377  make_with_value<Variables<tmpl::list<
379  interface_unit_normal,
381  fluxes_computer.apply(
383  &get<::Tags::Flux<FieldTags, tmpl::size_t<Dim>, Frame::Inertial>>(
384  principal_dirichlet_auxiliary_field_fluxes))...,
385  fluxes_args..., *numerical_flux_for_auxiliary_fields...);
386  const auto assemble_dirichlet_penalty = [
387  &interface_unit_normal, &penalty,
388  &principal_dirichlet_auxiliary_field_fluxes
389  ](auto field_tag_v, const auto numerical_flux_for_field) noexcept {
390  using field_tag = decltype(field_tag_v);
391  normal_dot_flux(
392  numerical_flux_for_field, interface_unit_normal,
393  get<::Tags::Flux<field_tag, tmpl::size_t<Dim>, Frame::Inertial>>(
394  principal_dirichlet_auxiliary_field_fluxes));
395  for (auto it = numerical_flux_for_field->begin();
396  it != numerical_flux_for_field->end(); it++) {
397  *it *= 2 * penalty;
398  }
399  };
401  assemble_dirichlet_penalty(FieldTags{}, numerical_flux_for_fields));
402  }
403 
404  private:
405  double penalty_parameter_ = std::numeric_limits<double>::signaling_NaN();
406 };
407 
408 } // namespace NumericalFluxes
409 } // namespace dg
410 } // namespace elliptic
Prefix indicating a boundary unit normal vector dotted into the flux.
Definition: Prefixes.hpp:124
The normalized (co)vector represented by Tag.
Definition: Magnitude.hpp:121
#define EXPAND_PACK_LEFT_TO_RIGHT(...)
Expand a parameter pack evaluating the terms from left to right.
Definition: TMPL.hpp:562
Definition: Digraph.hpp:11
T signaling_NaN(T... args)
Functionality related to solving elliptic partial differential equations.
Definition: InitializeAnalyticSolution.hpp:29
Functionality related to discontinuous Galerkin schemes.
Definition: ApplyBoundaryFluxesLocalTimeStepping.hpp:33
Defines classes and functions for making classes creatable from input files.
Define prefixes for DataBox tags.
Tags for the DataBox inherit from this type.
Definition: Tag.hpp:23
Prefix indicating the divergence.
Definition: Divergence.hpp:46
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.
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:30
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:42
Prefix indicating a flux.
Definition: Prefixes.hpp:56
Defines class Variables.
Definition: DataBoxTag.hpp:27
The internal penalty flux for first-order elliptic equations.
Definition: FirstOrderInternalPenalty.hpp:134
Defines classes for Tensor.
Declares function unnormalized_face_normal.
::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
Stores a collection of function values.
Definition: DataVector.hpp:42
Wraps the template metaprogramming library used (brigand)
typename DataBox_detail::item_type_impl< TagList, Tag >::type item_type
Get the type that can be written to the Tag. If it is a base tag then a TagList must be passed as a s...
Definition: DataBoxTag.hpp:246
Defines functions and classes from the GSL.
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, but it may be necessary to perform the conversion explicitly when type deduction is desired.
Definition: Gsl.hpp:879
Defines tags related to domain quantities.
Marks an item as being a prefix to another tag.
Definition: Tag.hpp:66
Definition: IndexType.hpp:44
Defines metafunctions used by Tensor.
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Prefix indicating a boundary unit normal vector dotted into the numerical flux.
Definition: Prefixes.hpp:138
Defines functions and tags for taking a divergence.
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182