Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <cstddef>
7 : #include <memory>
8 : #include <optional>
9 :
10 : #include "DataStructures/DataBox/Tag.hpp"
11 : #include "DataStructures/Tensor/TypeAliases.hpp"
12 : #include "Domain/FaceNormal.hpp"
13 : #include "Evolution/BoundaryCorrection.hpp"
14 : #include "Evolution/Systems/GeneralizedHarmonic/Tags.hpp"
15 : #include "NumericalAlgorithms/DiscontinuousGalerkin/Formulation.hpp"
16 : #include "Options/String.hpp"
17 : #include "PointwiseFunctions/GeneralRelativity/GeneralizedHarmonic/ConstraintDampingTags.hpp"
18 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
19 : #include "Utilities/Serialization/CharmPupable.hpp"
20 : #include "Utilities/TMPL.hpp"
21 :
22 : /// \cond
23 : class DataVector;
24 : namespace gsl {
25 : template <typename T>
26 : class not_null;
27 : } // namespace gsl
28 : namespace PUP {
29 : class er;
30 : } // namespace PUP
31 : /// \endcond
32 :
33 1 : namespace gh::BoundaryCorrections {
34 : /*!
35 : * \brief Computes the generalized harmonic upwind multipenalty boundary
36 : * correction using the averaged fields across the interface.
37 : *
38 : * The correction is given by
39 : *
40 : * \f{equation}{
41 : * D_\beta =
42 : * T^{-1}_{\beta\hat{\beta}}(u^{\text{avg}})
43 : * \Lambda^-_{\hat{\beta}\hat{\alpha}}(u^{\text{avg}})
44 : * T_{\hat{\alpha}\alpha}(u^{\text{avg}})
45 : * \Delta u_\alpha,
46 : * \f}
47 : *
48 : * where $u_\alpha$ are the evolved fields, $\Lambda^-$ is the
49 : * diagonal matrix of incoming characteristic speeds,
50 : * $T_{\hat{\alpha}\alpha}$ transforms from evolved to characteristic
51 : * fields, and $\Delta u_\alpha = u_\alpha^{\text{ext}} -
52 : * u_\alpha^{\text{int}}$ is the discontinuity in the evolved fields.
53 : * All factors except the boundary discontinuity are calculated using
54 : * the averaged of the internal and external evolved fields:
55 : *
56 : * \f{equation}
57 : * u^{\text{avg}} = \frac{u^{\text{int}} + u^{\text{ext}}}{2}.
58 : * \f}
59 : *
60 : * For the first-order generalized harmonic system the correction is:
61 : *
62 : * \f{align}{
63 : * D_{g_{ab}} &= \lambda_{v^g}^- \Delta g_{ab} \\
64 : * D_{\Pi_{ab}}
65 : * &= \left(\lambda_{v^g}^-
66 : * - \frac{\lambda_{v^+}^- + \lambda_{v^-}^-}{2}\right)
67 : * \gamma_2 \Delta g_{ab}
68 : * + \frac{\lambda_{v^+}^- + \lambda_{v^-}^-}{2} \Delta \Pi_{ab}
69 : * + \frac{\lambda_{v^+}^- - \lambda_{v^-}^-}{2} n^j \Delta \Phi_{jab} \\
70 : * D_{\Phi_{iab}}
71 : * &= \lambda_{v^0}^- (\delta_i^j - n_i n^j) \Delta \Phi_{jab}
72 : * + \frac{\lambda_{v^+}^- - \lambda_{v^-}^-}{2} n_i
73 : * (\Delta \Pi_{ab} - \gamma_2 \Delta g_{ab})
74 : * + \frac{\lambda_{v^+}^- + \lambda_{v^-}^-}{2} n_i n^j \Delta \Phi_{jab}.
75 : * \f}
76 : *
77 : * In the above expressions, the outgoing face normal $n_i$ is
78 : * normalized using the average metric from the two sides. Quantities
79 : * such as $\gamma_2$ that should analytically be the same on both
80 : * sides of the interface are also averaged, as they can differ
81 : * because of truncation errors during projection.
82 : *
83 : * The characteristic speeds are given by
84 : *
85 : * \f{align}{
86 : * \lambda_{v^g} &= - (1 + \gamma_1) (\beta^i + v^i_g) n_i, \\
87 : * \lambda_{v^0} &= - (\beta^i + v^i_g) n_i, \\
88 : * \lambda_{v^\pm} &= \pm \alpha - (\beta^i + v^i_g) n_i,
89 : * \f}
90 : *
91 : * where \f$v_g^i\f$ is the mesh velocity, with the incoming speeds
92 : *
93 : * \f{equation}{
94 : * \lambda^-_{\hat{\alpha}} =
95 : * \begin{cases}
96 : * \lambda_{\hat{\alpha}} & \lambda_{\hat{\alpha}} < 0 \\
97 : * 0 & \text{otherwise.}
98 : * \end{cases}
99 : * \f}
100 : */
101 : template <size_t Dim>
102 1 : class AveragedUpwindPenalty final : public evolution::BoundaryCorrection {
103 : public:
104 0 : using options = tmpl::list<>;
105 0 : static constexpr Options::String help = {
106 : "Computes the upwind penalty boundary correction term for the "
107 : "generalized harmonic system using the averaged fields from the two "
108 : "sides of the interface."};
109 :
110 0 : AveragedUpwindPenalty() = default;
111 0 : AveragedUpwindPenalty(const AveragedUpwindPenalty&) = default;
112 0 : AveragedUpwindPenalty& operator=(const AveragedUpwindPenalty&) = default;
113 0 : AveragedUpwindPenalty(AveragedUpwindPenalty&&) = default;
114 0 : AveragedUpwindPenalty& operator=(AveragedUpwindPenalty&&) = default;
115 0 : ~AveragedUpwindPenalty() override = default;
116 :
117 : /// \cond
118 : explicit AveragedUpwindPenalty(CkMigrateMessage* msg);
119 : using PUP::able::register_constructor;
120 : WRAPPED_PUPable_decl_template(AveragedUpwindPenalty); // NOLINT
121 : /// \endcond
122 0 : void pup(PUP::er& p) override; // NOLINT
123 :
124 0 : std::unique_ptr<BoundaryCorrection> get_clone() const override;
125 :
126 0 : struct MeshVelocity : db::SimpleTag {
127 0 : using type = tnsr::I<DataVector, Dim, Frame::Inertial>;
128 : };
129 :
130 0 : using dg_package_field_tags =
131 : tmpl::list<gr::Tags::SpacetimeMetric<DataVector, Dim>,
132 : gh::Tags::Pi<DataVector, Dim>, gh::Tags::Phi<DataVector, Dim>,
133 : gh::Tags::ConstraintGamma1, gh::Tags::ConstraintGamma2,
134 : domain::Tags::UnnormalizedFaceNormal<Dim>, MeshVelocity>;
135 : // We don't need lapse and shift, but this list must be exactly this
136 : // or the boundary condition code doesn't compile.
137 0 : using dg_package_data_temporary_tags =
138 : tmpl::list<gh::Tags::ConstraintGamma1, gh::Tags::ConstraintGamma2,
139 : gr::Tags::Lapse<DataVector>, gr::Tags::Shift<DataVector, Dim>>;
140 0 : using dg_package_data_primitive_tags = tmpl::list<>;
141 0 : using dg_package_data_volume_tags = tmpl::list<>;
142 0 : using dg_boundary_terms_volume_tags = tmpl::list<>;
143 :
144 0 : double dg_package_data(
145 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
146 : packaged_spacetime_metric,
147 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*> packaged_pi,
148 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*> packaged_phi,
149 : gsl::not_null<Scalar<DataVector>*> packaged_constraint_gamma1,
150 : gsl::not_null<Scalar<DataVector>*> packaged_constraint_gamma2,
151 : gsl::not_null<tnsr::i<DataVector, Dim, Frame::Inertial>*> packaged_normal,
152 : gsl::not_null<tnsr::I<DataVector, Dim, Frame::Inertial>*>
153 : packaged_mesh_velocity,
154 :
155 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric,
156 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi,
157 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi,
158 :
159 : const Scalar<DataVector>& constraint_gamma1,
160 : const Scalar<DataVector>& constraint_gamma2,
161 : const Scalar<DataVector>& lapse, const tnsr::I<DataVector, Dim>& shift,
162 :
163 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_covector,
164 : const tnsr::I<DataVector, Dim, Frame::Inertial>& normal_vector,
165 : const std::optional<tnsr::I<DataVector, Dim, Frame::Inertial>>&
166 : mesh_velocity,
167 : const std::optional<Scalar<DataVector>>& normal_dot_mesh_velocity) const;
168 :
169 0 : void dg_boundary_terms(
170 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
171 : boundary_correction_spacetime_metric,
172 : gsl::not_null<tnsr::aa<DataVector, Dim, Frame::Inertial>*>
173 : boundary_correction_pi,
174 : gsl::not_null<tnsr::iaa<DataVector, Dim, Frame::Inertial>*>
175 : boundary_correction_phi,
176 :
177 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric_int,
178 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi_int,
179 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi_int,
180 : const Scalar<DataVector>& constraint_gamma1_int,
181 : const Scalar<DataVector>& constraint_gamma2_int,
182 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_int,
183 : const tnsr::I<DataVector, Dim, Frame::Inertial>& mesh_velocity_int,
184 :
185 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& spacetime_metric_ext,
186 : const tnsr::aa<DataVector, Dim, Frame::Inertial>& pi_ext,
187 : const tnsr::iaa<DataVector, Dim, Frame::Inertial>& phi_ext,
188 : const Scalar<DataVector>& constraint_gamma1_ext,
189 : const Scalar<DataVector>& constraint_gamma2_ext,
190 : const tnsr::i<DataVector, Dim, Frame::Inertial>& normal_ext,
191 : const tnsr::I<DataVector, Dim, Frame::Inertial>& mesh_velocity_ext,
192 :
193 : dg::Formulation dg_formulation) const;
194 : };
195 :
196 : template <size_t Dim>
197 0 : bool operator==(const AveragedUpwindPenalty<Dim>& lhs,
198 : const AveragedUpwindPenalty<Dim>& rhs);
199 : template <size_t Dim>
200 0 : bool operator!=(const AveragedUpwindPenalty<Dim>& lhs,
201 : const AveragedUpwindPenalty<Dim>& rhs);
202 : } // namespace gh::BoundaryCorrections
|
