Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <limits>
7 : #include <memory>
8 : #include <optional>
9 :
10 : #include "DataStructures/DataBox/Prefixes.hpp"
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Domain/Tags.hpp"
13 : #include "Evolution/BoundaryCorrection.hpp"
14 : #include "Evolution/Systems/SecondOrderScalarWave/Tags.hpp"
15 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
16 : #include "Options/String.hpp"
17 : #include "Utilities/Serialization/CharmPupable.hpp"
18 : #include "Utilities/TMPL.hpp"
19 :
20 : /// \cond
21 : class DataVector;
22 : namespace gsl {
23 : template <typename T>
24 : class not_null;
25 : } // namespace gsl
26 : namespace PUP {
27 : class er;
28 : } // namespace PUP
29 : /// \endcond
30 :
31 : namespace SecondOrderScalarWave::BoundaryCorrections {
32 : /*!
33 : * \brief A Lax-Friedrichs-type boundary correction for the second-order
34 : * scalar wave system, providing the numerical fluxes for the local
35 : * discontinuous Galerkin (LDG) scheme.
36 : *
37 : * Each evaluation of the system's time derivative proceeds in two passes: the
38 : * auxiliary pass computes \f$\Phi_i = \partial_i \Psi\f$, then the physical
39 : * pass evolves \f$\partial_t \Pi = -\partial_i \Phi^i\f$ and
40 : * \f$\partial_t \Psi = -\Pi\f$. This class supplies the numerical flux for
41 : * each pass: the `dg_auxiliary_*` interface couples neighboring elements in
42 : * the auxiliary pass, the standard `dg_*` interface in the physical pass.
43 : *
44 : * Below, \f$n^i\f$ denotes the interior element's outward-directed unit face
45 : * normal. \f$\{\{\Psi\}\}\f$ denotes the central average of \f$\Psi\f$ across
46 : * the interface, and
47 : * \f$[[\Psi]]^i=n^i\left(\Psi^\text{int}-\Psi^\text{ext}\right)\f$
48 : * denotes the jump of \f$\Psi\f$ across the interface.
49 : *
50 : * ### Auxiliary pass
51 : *
52 : * Integrating \f$\Phi_i = \partial_i \Psi\f$ against a test
53 : * function \f$l\f$ by parts twice over an element \f$K\f$ gives
54 : *
55 : * \f{align*}{
56 : * \int_K \Phi_i l = \int_K (\partial_i \Psi)\, l
57 : * + \oint_{\partial K} l n_i \left(\Psi^* - \Psi\right),
58 : * \f}
59 : *
60 : * where
61 : *
62 : * \f{align*}{
63 : * \Psi^* = \{\{ \Psi \}\}
64 : * \f}
65 : *
66 : * is the central flux.
67 : *
68 : * ### Physical pass
69 : *
70 : * Integrating \f$\partial_t \Pi = -\partial_i \Phi^i\f$ by parts twice gives
71 : *
72 : * \f{align*}{
73 : * \int_K (\partial_t \Pi)\, l = -\int_K (\partial_i f_\Phi^i)\, l
74 : * - \oint_{\partial K}
75 : * \left[(f_\Phi^i)^* - f_\Phi^i\right] l n_i ,
76 : * \f}
77 : *
78 : * where
79 : *
80 : * \f{align*}{
81 : * f_\Phi^i = \Phi^i,
82 : * \f}
83 : *
84 : * \f{align*}{
85 : * (f_\Phi^i)^* = \{\{ f_\Phi^i \}\} + \frac{\tau}{2} [[ \Pi ]]^i
86 : * \f}
87 : *
88 : * \f$\partial_t \Psi = -\Pi\f$ contains no spatial derivative so receives
89 : * no boundary correction.
90 : */
91 : template <size_t Dim>
92 1 : class LaxFriedrichs final : public evolution::BoundaryCorrection {
93 : public:
94 0 : struct Tau {
95 0 : using type = double;
96 0 : static constexpr Options::String help = {
97 : "The penalty parameter tau for the physical numerical flux."};
98 : };
99 :
100 0 : using options = tmpl::list<Tau>;
101 0 : static constexpr Options::String help = {
102 : "Boundary correction to the LDG method using the LaxFriedrichs numerical "
103 : "flux."};
104 :
105 0 : LaxFriedrichs() = default;
106 0 : explicit LaxFriedrichs(double tau);
107 0 : LaxFriedrichs(const LaxFriedrichs&) = default;
108 0 : LaxFriedrichs& operator=(const LaxFriedrichs&) = default;
109 0 : LaxFriedrichs(LaxFriedrichs&&) = default;
110 0 : LaxFriedrichs& operator=(LaxFriedrichs&&) = default;
111 0 : ~LaxFriedrichs() override = default;
112 :
113 : /// \cond
114 : explicit LaxFriedrichs(CkMigrateMessage* msg);
115 : using PUP::able::register_constructor;
116 : WRAPPED_PUPable_decl_template(LaxFriedrichs); // NOLINT
117 : /// \endcond
118 0 : void pup(PUP::er& p) override; // NOLINT
119 :
120 0 : std::unique_ptr<BoundaryCorrection> get_clone() const override;
121 :
122 0 : using dg_package_field_tags = tmpl::list<Tags::Pi, Tags::NormalDotPhi>;
123 0 : using dg_package_data_temporary_tags = tmpl::list<>;
124 0 : using dg_package_data_volume_tags = tmpl::list<>;
125 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
126 :
127 0 : double dg_package_data(
128 : gsl::not_null<Scalar<DataVector>*> packaged_pi,
129 : gsl::not_null<Scalar<DataVector>*> packaged_normal_dot_phi,
130 :
131 : const Scalar<DataVector>& /*psi*/, const Scalar<DataVector>& pi,
132 : const tnsr::i<DataVector, Dim, Frame::Inertial>& phi,
133 :
134 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
135 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
136 : /*mesh_velocity*/,
137 : const std::optional<Scalar<DataVector>>& /*normal_dot_mesh_velocity*/)
138 : const;
139 :
140 0 : void dg_boundary_terms(
141 : gsl::not_null<Scalar<DataVector>*> psi_boundary_correction,
142 : gsl::not_null<Scalar<DataVector>*> pi_boundary_correction,
143 :
144 : const Scalar<DataVector>& pi_int,
145 : const Scalar<DataVector>& normal_dot_phi_int,
146 :
147 : const Scalar<DataVector>& pi_ext,
148 : const Scalar<DataVector>& normal_dot_phi_ext,
149 :
150 : dg::Formulation /*dg_formulation*/) const;
151 :
152 0 : using dg_auxiliary_package_field_tags = tmpl::list<Tags::PsiTimesNormal<Dim>>;
153 0 : using dg_auxiliary_package_data_temporary_tags = tmpl::list<>;
154 0 : using dg_auxiliary_package_data_volume_tags = tmpl::list<>;
155 0 : using dg_auxiliary_boundary_terms_volume_tags = tmpl::list<>;
156 :
157 0 : double dg_auxiliary_package_data(
158 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
159 : psi_times_normal,
160 :
161 : const Scalar<DataVector>& psi, const Scalar<DataVector>& /*pi*/,
162 :
163 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
164 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
165 : /*mesh_velocity*/,
166 : const std::optional<Scalar<DataVector>>& /*normal_dot_mesh_velocity*/)
167 : const;
168 :
169 0 : void dg_auxiliary_boundary_terms(
170 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*>
171 : phi_boundary_correction,
172 :
173 : const tnsr::i<DataVector, Dim, Frame::Inertial>& psi_times_normal_int,
174 :
175 : const tnsr::i<DataVector, Dim, Frame::Inertial>& psi_times_normal_ext,
176 :
177 : dg::Formulation /*dg_formulation*/) const;
178 :
179 : private:
180 : template <size_t LocalDim>
181 : // NOLINTNEXTLINE(readability-redundant-declaration)
182 0 : friend bool operator==(const LaxFriedrichs<LocalDim>& lhs,
183 : const LaxFriedrichs<LocalDim>& rhs);
184 :
185 0 : double tau_ = std::numeric_limits<double>::signaling_NaN();
186 : };
187 :
188 : template <size_t Dim>
189 0 : bool operator!=(const LaxFriedrichs<Dim>& lhs, const LaxFriedrichs<Dim>& rhs);
190 : } // namespace SecondOrderScalarWave::BoundaryCorrections
|
