Equations.hpp
1 // Distributed under the MIT License.
2 // See LICENSE.txt for details.
3 
4 #pragma once
5 
6 #include <cstddef>
7 #include <string>
8 
11 #include "DataStructures/Variables.hpp" // IWYU pragma: keep
13 #include "Domain/FaceNormal.hpp"
14 #include "Domain/Tags.hpp"
16 #include "Options/Options.hpp"
17 #include "Utilities/Gsl.hpp"
18 #include "Utilities/TMPL.hpp"
19 // IWYU pragma: no_forward_declare Tags::deriv
20 // IWYU pragma: no_forward_declare Tensor
21 // IWYU pragma: no_forward_declare Variables
22 
23 /// \cond
24 class DataVector;
25 template <size_t>
26 class Mesh;
27 namespace Tags {
28 template <typename>
29 struct Normalized;
30 } // namespace Tags
31 namespace LinearSolver {
32 namespace Tags {
33 template <typename>
34 struct Operand;
35 } // namespace Tags
36 } // namespace LinearSolver
37 namespace Poisson {
38 struct Field;
39 template <size_t>
40 struct AuxiliaryField;
41 } // namespace Poisson
42 namespace PUP {
43 class er;
44 } // namespace PUP
45 /// \endcond
46 
47 namespace Poisson {
48 
49 /*!
50  * \brief The bulk contribution to the linear operator action for the
51  * first order formulation of the Poisson equation.
52  *
53  * \details The bulk contribution for the equation sourced by \f$f(x)\f$ is
54  * \f$-\nabla \cdot \boldsymbol{v}(x)\f$ and the one for the auxiliary equation
55  * is \f$\nabla u(x) - \boldsymbol{v}(x)\f$ (see `Poisson::FirstOrderSystem`).
56  */
57 template <size_t Dim>
59  using argument_tags =
60  tmpl::list<Tags::deriv<LinearSolver::Tags::Operand<Field>,
61  tmpl::size_t<Dim>, Frame::Inertial>,
66  static void apply(
67  gsl::not_null<Scalar<DataVector>*> operator_for_field_source,
68  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
69  operator_for_auxiliary_field_source,
70  const tnsr::i<DataVector, Dim, Frame::Inertial>& grad_field,
71  const tnsr::I<DataVector, Dim, Frame::Inertial>& auxiliary_field,
72  const Mesh<Dim>& mesh,
73  const InverseJacobian<DataVector, Dim, Frame::Logical, Frame::Inertial>&
74  inverse_jacobian) noexcept;
75 };
76 
77 /*!
78  * \brief The interface normal dotted into the fluxes for the first order
79  * formulation of the Poisson equation.
80  *
81  * \details For the sourced equation this is \f$-\boldsymbol{n} \cdot
82  * \boldsymbol{v}(x)\f$ and for the auxiliary equation it is
83  * \f$\boldsymbol{n} u\f$ (see `Poisson::FirstOrderSystem` and `dg::lift_flux`).
84  */
85 template <size_t Dim>
87  using argument_tags =
88  tmpl::list<LinearSolver::Tags::Operand<Field>,
91  static void apply(
92  gsl::not_null<Scalar<DataVector>*> normal_dot_flux_for_field,
93  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
94  normal_dot_flux_for_auxiliary_field,
95  const Scalar<DataVector>& field,
96  const tnsr::I<DataVector, Dim, Frame::Inertial>& auxiliary_field,
97  const tnsr::i<DataVector, Dim, Frame::Inertial>&
98  interface_unit_normal) noexcept;
99 };
100 
101 /*!
102  * \brief The internal penalty flux for the first oder formulation of the
103  * Poisson equation.
104  *
105  * \details For the sourced equation this is \f$-\frac{\nabla u_\mathrm{int} +
106  * \nabla u_\mathrm{ext}}{2} + \sigma \left(u_\mathrm{int} -
107  * u_\mathrm{ext}\right)\f$ and for the auxiliary equation it is \f$\frac{
108  * u_\mathrm{int} + u_\mathrm{ext}}{2}\f$. The penalty factor \f$\sigma\f$ is
109  * responsible for removing zero eigenmodes and impacts the conditioning of the
110  * linear operator to be solved. It can be chosen as
111  * \f$\sigma=C\frac{N_\mathrm{points}^2}{h}\f$ where \f$N_\mathrm{points}\f$ is
112  * the number of collocation points (i.e. the polynomial degree plus 1),
113  * \f$h\f$ is a measure of the element size in inertial coordinates and \f$C\leq
114  * 1\f$ is a free parameter (see e.g. \cite HesthavenWarburton, section 7.2).
115  */
116 template <size_t Dim>
118  public:
120  using type = double;
121  // Currently this is used as the full prefactor to the penalty term. When it
122  // becomes possible to compute a measure of the size $h$ of an element and
123  // the number of collocation points $p$ on both sides of the mortar, this
124  // should be changed to be just the parameter multiplying $\frac{p^2}{h}$.
125  static constexpr OptionString help = {
126  "The prefactor to the penalty term of the flux."};
127  };
128  using options = tmpl::list<PenaltyParameter>;
129  static constexpr OptionString help = {
130  "Computes the internal penalty flux for a Poisson system."};
131 
132  FirstOrderInternalPenaltyFlux() = default;
133  explicit FirstOrderInternalPenaltyFlux(double penalty_parameter)
134  : penalty_parameter_(penalty_parameter) {}
135 
136  // clang-tidy: non-const reference
137  void pup(PUP::er& p) noexcept { p | penalty_parameter_; } // NOLINT
138 
140  using type = tnsr::i<DataVector, Dim, Frame::Inertial>;
141  static std::string name() noexcept { return "NormalTimesFieldFlux"; }
142  };
143 
145  using type = Scalar<DataVector>;
146  static std::string name() noexcept { return "NormalDotGradFieldFlux"; }
147  };
148 
149  // These tags are sliced to the interface of the element and passed to
150  // `package_data` to provide the data needed to compute the numerical fluxes.
151  using argument_tags =
152  tmpl::list<LinearSolver::Tags::Operand<Field>,
153  Tags::deriv<LinearSolver::Tags::Operand<Field>,
154  tmpl::size_t<Dim>, Frame::Inertial>,
156 
157  // This is the data needed to compute the numerical flux.
158  // `SendBoundaryFluxes` calls `package_data` to store these tags in a
159  // Variables. Local and remote values of this data are then combined in the
160  // `()` operator.
161  using package_tags = tmpl::list<LinearSolver::Tags::Operand<Field>,
163 
164  // Following the packaged_data pointer, this function expects as arguments the
165  // types in `argument_tags`.
166  void package_data(gsl::not_null<Variables<package_tags>*> packaged_data,
167  const Scalar<DataVector>& field,
168  const tnsr::i<DataVector, Dim, Frame::Inertial>& grad_field,
169  const tnsr::i<DataVector, Dim, Frame::Inertial>&
170  interface_unit_normal) const noexcept;
171 
172  // This function combines local and remote data to the numerical fluxes.
173  // The numerical fluxes as not-null pointers are the first arguments. The
174  // other arguments are the packaged types for the interior side followed by
175  // the packaged types for the exterior side.
176  void operator()(
177  gsl::not_null<Scalar<DataVector>*> numerical_flux_for_field,
178  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
179  numerical_flux_for_auxiliary_field,
180  const Scalar<DataVector>& field_interior,
181  const tnsr::i<DataVector, Dim, Frame::Inertial>&
182  normal_times_field_interior,
183  const Scalar<DataVector>& normal_dot_grad_field_interior,
184  const Scalar<DataVector>& field_exterior,
185  const tnsr::i<DataVector, Dim, Frame::Inertial>&
186  minus_normal_times_field_exterior,
187  const Scalar<DataVector>& minus_normal_dot_grad_field_exterior) const
188  noexcept;
189 
190  // This function computes the boundary contributions from Dirichlet boundary
191  // conditions. This data is what remains to be added to the boundaries when
192  // homogeneous (i.e. zero) boundary conditions are assumed in the calculation
193  // of the numerical fluxes, but we wish to impose inhomogeneous (i.e. nonzero)
194  // boundary conditions. Since this contribution does not depend on the
195  // numerical field values, but only on the Dirichlet boundary data, it may be
196  // added as contribution to the source of the elliptic systems. Then, it
197  // remains to solve the homogeneous problem with the modified source.
198  // The first arguments to this function are the boundary contributions to
199  // compute as not-null pointers, in the order they appear in the
200  // `system::fields_tag`. They are followed by the field values of the tags in
201  // `system::impose_boundary_conditions_on_fields`. The last argument is the
202  // normalized unit covector to the element face.
203  void compute_dirichlet_boundary(
204  gsl::not_null<Scalar<DataVector>*> numerical_flux_for_field,
205  gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
206  numerical_flux_for_auxiliary_field,
207  const Scalar<DataVector>& dirichlet_field,
208  const tnsr::i<DataVector, Dim, Frame::Inertial>& interface_unit_normal)
209  const noexcept;
210 
211  private:
212  double penalty_parameter_{};
213 };
214 
215 } // namespace Poisson
The bulk contribution to the linear operator action for the first order formulation of the Poisson eq...
Definition: Equations.hpp:58
Definition: Strahlkorper.hpp:14
The internal penalty flux for the first oder formulation of the Poisson equation. ...
Definition: Equations.hpp:117
The normalized (co)vector represented by Tag.
Definition: Magnitude.hpp:93
Functionality for solving linear systems of equations.
Definition: TerminateIfConverged.hpp:22
Holds the number of grid points, basis, and quadrature in each direction of the computational grid...
Definition: Mesh.hpp:49
Defines classes and functions for making classes creatable from input files.
Items related to solving a Poisson equation .
Definition: Actions.hpp:6
Tags for the DataBox inherit from this type.
Definition: DataBoxTag.hpp:65
The interface normal dotted into the fluxes for the first order formulation of the Poisson equation...
Definition: Equations.hpp:86
constexpr auto apply(F &&f, const DataBox< BoxTags > &box, Args &&... args)
Apply the function f with argument Tags TagsList from DataBox box
Definition: DataBox.hpp:1595
The operand that the local linear operator is applied to.
Definition: Tags.hpp:40
const char *const OptionString
The string used in option structs.
Definition: Options.hpp:27
Computes the inverse Jacobian of the map held by MapTag at the coordinates held by SourceCoordsTag...
Definition: Tags.hpp:125
Defines functions computing partial derivatives.
Defines class Variables.
Definition: DataBoxTag.hpp:29
Defines a list of useful type aliases for tensors.
Declares function unnormalized_face_normal.
The coordinates in a given frame.
Definition: Tags.hpp:95
Stores a collection of function values.
Definition: DataVector.hpp:46
Wraps the template metaprogramming library used (brigand)
The computational grid of the Element in the DataBox.
Definition: Tags.hpp:75
Defines functions and classes from the GSL.
Defines tags related to domain quantities.
Definition: IndexType.hpp:44
Defines helper functions for use with Variables class.
Defines classes SimpleTag, PrefixTag, ComputeTag and several functions for retrieving tag info...
Tensor< T, Symmetry<>, index_list<> > Scalar
Scalar type.
Definition: TypeAliases.hpp:21
Require a pointer to not be a nullptr
Definition: ConservativeFromPrimitive.hpp:12