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