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 
15 #include "Domain/FaceNormal.hpp"
16 #include "Domain/Tags.hpp"
17 #include "NumericalAlgorithms/DiscontinuousGalerkin/NormalDotFlux.hpp"
19 #include "Options/Options.hpp"
20 #include "Utilities/Gsl.hpp"
21 #include "Utilities/TMPL.hpp"
22 
23 /// \cond
24 class DataVector;
25 namespace Tags {
26 template <typename>
27 struct Normalized;
28 } // namespace Tags
29 /// \endcond
30 
31 namespace elliptic {
32 namespace dg {
33 namespace NumericalFluxes {
34 
35 /*!
36  * \ingroup DiscontinuousGalerkinGroup
37  * \ingroup NumericalFluxesGroup
38  * \brief The internal penalty flux for first-order elliptic equations.
39  *
40  * \details Computes the internal penalty numerical flux (see e.g.
41  * \cite HesthavenWarburton, section 7.2) dotted with the interface unit normal.
42  *
43  * We implement here a suggested generalization of the internal penalty flux
44  * for any set of elliptic PDEs up to second order in derivatives. It is
45  * designed for fluxes (i.e. principal parts of the PDEs) that may depend on the
46  * dynamic fields, but do so at most linearly. This is the case for the velocity
47  * potential equation of binary neutron stars, for example. The linearity is
48  * only necessary to impose inhomogeneous boundary conditions as contributions
49  * to the fixed source (see
50  * `elliptic::dg::Actions::ImposeInhomogeneousBoundaryConditionsOnSource`).
51  *
52  * We reduce the second-order elliptic PDEs to first-order form by introducing
53  * an _auxiliary_ variable \f$v\f$ for each _primal_ variable \f$u\f$ (i.e. for
54  * each variable whose second derivative appears in the equations). Then, the
55  * equations take the _flux-form_
56  *
57  * \f{align}
58  * -\partial_i F^i_A + S_A = f_A
59  * \f}
60  *
61  * where the fluxes \f$F^i_A(u,v)\f$ and sources \f$S_A(u,v)\f$ are indexed
62  * (with capital letters) by the variables and we have defined \f$f_A(x)\f$ as
63  * those sources that are independent of the variables. Note that the fluxes are
64  * functions of the primal and auxiliary variables; e.g. for a Poisson system
65  * \f$\{u,v_i\}\f$ on a spatial metric \f$\gamma_{ij}\f$ they are simply
66  * \f$F^i_u(v)=\sqrt{\gamma}\gamma^{ij} v_j\f$ and
67  * \f$F^i_{v_j}(u)=u\delta^i_j\f$ (see `Poisson::FirstOrderSystem`) or for an
68  * Elasticity system \f$\{\xi^i,S_{ij}\}\f$ with Young's tensor
69  * \f$Y^{ijkl}\f$ they are \f$F^i_{\xi^j}(S)=Y^{ijkl}S_{kl}\f$ and
70  * \f$F^i_{S_{jk}}(\xi)=\delta^i_{(j}\xi_{k)}\f$. Since each primal flux is
71  * commonly only a function of the corresponding auxiliary variable and
72  * vice-versa, we omit the other function arguments here to lighten the
73  * notational load.
74  *
75  * We now choose the internal penalty numerical fluxes \f$(n_i F^i_A)^*\f$ as
76  * follows for each primal variable \f$u\f$ and their corresponding auxiliary
77  * variable \f$v\f$:
78  *
79  * \f{align}
80  * (n_i F^i_u)^* &= \frac{1}{2} n_i \left( F^i_u(\partial_j
81  * F^j_v(u^\mathrm{int})) + F^i_u(\partial_j F^j_v(u^\mathrm{ext}))
82  * \right) - \sigma n_i \left(F^i_u(n_j F^j_v(u^\mathrm{int})) - F^i_u(
83  * n_j F^j_v(u^\mathrm{ext}))\right) \\
84  * (n_i F^i_v)^* &= \frac{1}{2} n_i \left(F^i_v(u^\mathrm{int}) +
85  * F^i_v(u^\mathrm{ext})\right)
86  * \f}
87  *
88  * Note that we have assumed \f$n^\mathrm{ext}_i=-n_i\f$ here, i.e. face normals
89  * don't depend on the dynamic variables (which may be discontinuous on element
90  * faces). This is the case for the problems we are expecting to solve, because
91  * those will be on fixed background metrics (e.g. a conformal metric for the
92  * XCTS system).
93  *
94  * Also note that the numerical fluxes intentionally don't depend on the
95  * auxiliary field values \f$v\f$. This property allows us to use the numerical
96  * fluxes also for the second-order (or _primal_) DG formulation, where we
97  * remove the need for an auxiliary variable. For the first-order system we
98  * could replace the divergence in \f$(n_i F^i_u)^*\f$ with \f$v\f$, which would
99  * result in a generalized "stabilized central flux" that is slightly less
100  * sparse than the internal penalty flux (see e.g. \cite HesthavenWarburton,
101  * section 7.2). We could also choose to ignore the fluxes in the penalty term,
102  * but preliminary tests suggest that this may hurt convergence.
103  *
104  * For a Poisson system (see above) this numeric flux reduces to the standard
105  * internal penalty flux (see e.g. \cite HesthavenWarburton, section 7.2, or
106  * \cite Arnold2002)
107  *
108  * \f{align}
109  * (n_i F^i_u)^* &= n_i v_i^* = \frac{1}{2} n_i \left(\partial_i u^\mathrm{int}
110  * + \partial_i u^\mathrm{ext}\right) - \sigma \left(u^\mathrm{int} -
111  * u^\mathrm{ext}\right) \\
112  * (n_i F^i_{v_j})^* &= n_j u^* = \frac{1}{2} n_j \left(u^\mathrm{int} +
113  * u^\mathrm{ext}\right)
114  * \f}
115  *
116  * where a sum over repeated indices is assumed, since the equation is
117  * formulated on a Euclidean geometry.
118  *
119  * This generalization of the internal penalty flux is based on unpublished work
120  * by Nils L. Fischer (nils.fischer@aei.mpg.de).
121  *
122  * The penalty factor \f$\sigma\f$ is responsible for removing zero eigenmodes
123  * and impacts the conditioning of the linear operator to be solved. It can be
124  * chosen as \f$\sigma=C\frac{N_\mathrm{points}^2}{h}\f$ where
125  * \f$N_\mathrm{points}\f$ is the number of collocation points (i.e. the
126  * polynomial degree plus 1), \f$h\f$ is a measure of the element size in
127  * inertial coordinates (both measured perpendicular to the interface under
128  * consideration) and \f$C\geq 1\f$ is a free parameter (see e.g.
129  * \cite HesthavenWarburton, section 7.2).
130  */
131 template <size_t Dim, typename FluxesComputerTag, typename FieldTagsList,
132  typename AuxiliaryFieldTagsList,
133  typename FluxesComputer = db::item_type<FluxesComputerTag>,
134  typename FluxesArgs = typename FluxesComputer::argument_tags>
136 
137 template <size_t Dim, typename FluxesComputerTag, typename... FieldTags,
138  typename... AuxiliaryFieldTags, typename FluxesComputer,
139  typename... FluxesArgs>
140 struct FirstOrderInternalPenalty<Dim, FluxesComputerTag,
141  tmpl::list<FieldTags...>,
142  tmpl::list<AuxiliaryFieldTags...>,
143  FluxesComputer, tmpl::list<FluxesArgs...>> {
144  private:
145  using fluxes_computer_tag = FluxesComputerTag;
146 
147  template <typename Tag>
148  struct NormalDotDivFlux : db::PrefixTag, db::SimpleTag {
149  using tag = Tag;
151  };
152 
153  template <typename Tag>
154  struct NormalDotNormalDotFlux : db::PrefixTag, db::SimpleTag {
155  using tag = Tag;
157  };
158 
159  public:
160  struct PenaltyParameter {
161  using type = double;
162  // Currently this is used as the full prefactor `sigma` to the penalty term.
163  // This means it needs to be chosen large enough so that the scheme is
164  // stable everywhere. A good estimate is the the largest
165  // `sigma > N_points^2 / h` (see class documentation) over all elements in
166  // the domain. Choosing `sigma` the same everywhere means it is an
167  // overestimate on non-uniform meshes where elements are large or polynomial
168  // degrees are small, but this only affects the conditioning of the scheme,
169  // i.e. it will converge slower but to the same solution. When it becomes
170  // possible to communicate non-tensors (the element size and the number of
171  // collocation points on either side of the mortar), this option will be
172  // changed to be just the free parameter `C` multiplying `N_points^2 / h`.
173  // This will improve conditioning on non-uniform meshes. Currently, the
174  // `packaged_data` is a `Variables` which can only hold tensors.
175  static constexpr OptionString help = {
176  "The prefactor to the penalty term of the flux."};
177  };
178  using options = tmpl::list<PenaltyParameter>;
179  static constexpr OptionString help = {
180  "The internal penalty flux for elliptic systems."};
181  static std::string name() noexcept { return "InternalPenalty"; }
182 
183  FirstOrderInternalPenalty() = default;
184  explicit FirstOrderInternalPenalty(const double penalty_parameter) noexcept
185  : penalty_parameter_(penalty_parameter) {}
186 
187  // clang-tidy: non-const reference
188  void pup(PUP::er& p) noexcept { p | penalty_parameter_; } // NOLINT
189 
190  // These tags are sliced to the interface of the element and passed to
191  // `package_data` to provide the data needed to compute the numerical fluxes.
192  using argument_tags =
193  tmpl::list<::Tags::NormalDotFlux<AuxiliaryFieldTags>...,
194  ::Tags::div<::Tags::Flux<AuxiliaryFieldTags, tmpl::size_t<Dim>,
195  Frame::Inertial>>...,
196  fluxes_computer_tag, FluxesArgs...,
198  using volume_tags = tmpl::list<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  void package_data(
212  const gsl::not_null<Variables<package_tags>*> packaged_data,
214  AuxiliaryFieldTags>>&... normal_dot_auxiliary_field_fluxes,
215  const db::item_type<::Tags::div<
216  ::Tags::Flux<AuxiliaryFieldTags, tmpl::size_t<Dim>,
217  Frame::Inertial>>>&... div_auxiliary_field_fluxes,
218  const FluxesComputer& fluxes_computer,
219  const db::item_type<FluxesArgs>&... fluxes_args,
220  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal)
221  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  void compute_dirichlet_boundary(
340  FieldTags>>*>... numerical_flux_for_fields,
342  AuxiliaryFieldTags>>*>... numerical_flux_for_auxiliary_fields,
343  const db::item_type<FieldTags>&... dirichlet_fields,
344  const FluxesComputer& fluxes_computer,
345  const db::item_type<FluxesArgs>&... fluxes_args,
346  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal)
347  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:122
The normalized (co)vector represented by Tag.
Definition: Magnitude.hpp:119
#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:27
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: DataBoxTag.hpp:64
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:29
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:54
Defines class Variables.
Definition: DataBoxTag.hpp:29
The internal penalty flux for first-order elliptic equations.
Definition: FirstOrderInternalPenalty.hpp:135
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:475
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: DataBoxTag.hpp:111
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:136
Defines functions and tags for taking a divergence.
Require a pointer to not be a nullptr
Definition: Gsl.hpp:182