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 <pup.h>
9 : #include <string>
10 : #include <vector>
11 :
12 : #include "DataStructures/Tensor/TypeAliases.hpp"
13 : #include "Domain/Tags.hpp"
14 : #include "Domain/Tags/FaceNormal.hpp"
15 : #include "Elliptic/BoundaryConditions/BoundaryCondition.hpp"
16 : #include "Elliptic/BoundaryConditions/BoundaryConditionType.hpp"
17 : #include "Elliptic/Systems/Xcts/FluxesAndSources.hpp"
18 : #include "Options/String.hpp"
19 : #include "Utilities/Gsl.hpp"
20 : #include "Utilities/Serialization/CharmPupable.hpp"
21 : #include "Utilities/TMPL.hpp"
22 :
23 : /// \cond
24 : class DataVector;
25 : /// \endcond
26 :
27 : namespace Xcts::BoundaryConditions {
28 :
29 : /*!
30 : * \brief Impose Robin boundary conditions at the outer boundary
31 : *
32 : * Impose $\partial_r(r\psi)=0$, $\partial_r(\alpha\psi)=0$, and
33 : * $\partial_r(r\beta_\mathrm{excess}^j)=0$ on the boundary. These Robin
34 : * boundary conditions incur an error of order $1/R^2$, where $R$ is the outer
35 : * radius of the domain. This allows to place the outer boundary much closer
36 : * to the center than for Dirichlet boundary conditions
37 : * (Xcts::BoundaryConditions::Flatness), which incur an error of order $1/R$.
38 : *
39 : * The Robin boundary conditions are imposed as Neumann-type as follows:
40 : *
41 : * \begin{align*}
42 : * (n^i \partial_i \psi)_N &= -\frac{\psi}{r}, \\
43 : * (n^i \partial_i (\alpha\psi))_N &= -\frac{\alpha\psi}{r}, \\
44 : * (n^i (L\beta)^{ij})_N = n^i (L\beta)^{ij} - \left[
45 : * \frac{\beta_\mathrm{excess}^j}{r}
46 : * + \frac{1}{3} n^j \frac{n^i \beta_\mathrm{excess}^i}{r}
47 : * + n^i \partial_i \beta_\mathrm{excess}^j
48 : * + \frac{1}{3} n^j n^i n_k \partial_i \beta_\mathrm{excess}^k
49 : * \right]
50 : * \end{align*}
51 : *
52 : * Here, the condition on the longitudinal shift is derived by imposing the
53 : * Robin boundary condition
54 : * $n^i\partial_i\beta_\mathrm{excess}^j=-\beta_\mathrm{excess}^j/r$ only on the
55 : * normal component of the shift gradient. To do this we can use the projection
56 : * operator $P_{ij}=\delta_{ij}-n_i n_j$ to set the shift gradient to
57 : * $\partial_i\beta_\mathrm{excess}^j=P_{ik}\partial_k\beta_\mathrm{excess}^j
58 : * -n_i \beta_\mathrm{excess}^j/r$ and then apply the longitudinal operator.
59 : *
60 : * \tparam EnabledEquations The subset of XCTS equations that are being solved
61 : */
62 : template <Xcts::Equations EnabledEquations>
63 1 : class Robin : public elliptic::BoundaryConditions::BoundaryCondition<3> {
64 : private:
65 0 : using Base = elliptic::BoundaryConditions::BoundaryCondition<3>;
66 :
67 : public:
68 0 : using options = tmpl::list<>;
69 0 : static constexpr Options::String help =
70 : "Impose Robin boundary conditions at the outer boundary. "
71 : "They incur an error of order 1/R^2, where R is the outer radius.";
72 :
73 0 : Robin() = default;
74 0 : Robin(const Robin&) = default;
75 0 : Robin& operator=(const Robin&) = default;
76 0 : Robin(Robin&&) = default;
77 0 : Robin& operator=(Robin&&) = default;
78 0 : ~Robin() override = default;
79 :
80 : /// \cond
81 : explicit Robin(CkMigrateMessage* m) : Base(m) {}
82 : using PUP::able::register_constructor;
83 : WRAPPED_PUPable_decl_template(Robin);
84 : /// \endcond
85 :
86 0 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition> get_clone()
87 : const override;
88 :
89 0 : std::vector<elliptic::BoundaryConditionType> boundary_condition_types()
90 : const override;
91 :
92 0 : using argument_tags =
93 : tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>,
94 : domain::Tags::FaceNormal<3, Frame::Inertial>>;
95 0 : using volume_tags = tmpl::list<>;
96 :
97 0 : void apply(gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
98 : gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
99 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
100 : const tnsr::I<DataVector, 3>& x,
101 : const tnsr::i<DataVector, 3>& face_normal) const;
102 :
103 0 : void apply(
104 : gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
105 : gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
106 : gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
107 : gsl::not_null<Scalar<DataVector>*>
108 : n_dot_lapse_times_conformal_factor_gradient,
109 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
110 : const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
111 : const tnsr::I<DataVector, 3>& x,
112 : const tnsr::i<DataVector, 3>& face_normal) const;
113 :
114 0 : void apply(
115 : gsl::not_null<Scalar<DataVector>*> conformal_factor_minus_one,
116 : gsl::not_null<Scalar<DataVector>*> lapse_times_conformal_factor_minus_one,
117 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess,
118 : gsl::not_null<Scalar<DataVector>*> n_dot_conformal_factor_gradient,
119 : gsl::not_null<Scalar<DataVector>*>
120 : n_dot_lapse_times_conformal_factor_gradient,
121 : gsl::not_null<tnsr::I<DataVector, 3>*> n_dot_longitudinal_shift_excess,
122 : const tnsr::i<DataVector, 3>& deriv_conformal_factor,
123 : const tnsr::i<DataVector, 3>& deriv_lapse_times_conformal_factor,
124 : const tnsr::iJ<DataVector, 3>& deriv_shift_excess,
125 : const tnsr::I<DataVector, 3>& x,
126 : const tnsr::i<DataVector, 3>& face_normal) const;
127 :
128 0 : using argument_tags_linearized =
129 : tmpl::list<domain::Tags::Coordinates<3, Frame::Inertial>,
130 : domain::Tags::FaceNormal<3, Frame::Inertial>>;
131 0 : using volume_tags_linearized = tmpl::list<>;
132 :
133 0 : void apply_linearized(
134 : gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
135 : gsl::not_null<Scalar<DataVector>*>
136 : n_dot_conformal_factor_gradient_correction,
137 : const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
138 : const tnsr::I<DataVector, 3>& x,
139 : const tnsr::i<DataVector, 3>& face_normal) const;
140 :
141 0 : void apply_linearized(
142 : gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
143 : gsl::not_null<Scalar<DataVector>*>
144 : lapse_times_conformal_factor_correction,
145 : gsl::not_null<Scalar<DataVector>*>
146 : n_dot_conformal_factor_gradient_correction,
147 : gsl::not_null<Scalar<DataVector>*>
148 : n_dot_lapse_times_conformal_factor_gradient_correction,
149 : const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
150 : const tnsr::i<DataVector, 3>&
151 : deriv_lapse_times_conformal_factor_correction,
152 : const tnsr::I<DataVector, 3>& x,
153 : const tnsr::i<DataVector, 3>& face_normal) const;
154 :
155 0 : void apply_linearized(
156 : gsl::not_null<Scalar<DataVector>*> conformal_factor_correction,
157 : gsl::not_null<Scalar<DataVector>*>
158 : lapse_times_conformal_factor_correction,
159 : gsl::not_null<tnsr::I<DataVector, 3>*> shift_excess_correction,
160 : gsl::not_null<Scalar<DataVector>*>
161 : n_dot_conformal_factor_gradient_correction,
162 : gsl::not_null<Scalar<DataVector>*>
163 : n_dot_lapse_times_conformal_factor_gradient_correction,
164 : gsl::not_null<tnsr::I<DataVector, 3>*>
165 : n_dot_longitudinal_shift_excess_correction,
166 : const tnsr::i<DataVector, 3>& deriv_conformal_factor_correction,
167 : const tnsr::i<DataVector, 3>&
168 : deriv_lapse_times_conformal_factor_correction,
169 : const tnsr::iJ<DataVector, 3>& deriv_shift_excess_correction,
170 : const tnsr::I<DataVector, 3>& x,
171 : const tnsr::i<DataVector, 3>& face_normal) const;
172 : };
173 :
174 : template <Xcts::Equations LocalEnabledEquations>
175 0 : bool operator==(const Robin<LocalEnabledEquations>& lhs,
176 : const Robin<LocalEnabledEquations>& rhs);
177 :
178 : template <Xcts::Equations EnabledEquations>
179 0 : bool operator!=(const Robin<EnabledEquations>& lhs,
180 : const Robin<EnabledEquations>& rhs);
181 :
182 : } // namespace Xcts::BoundaryConditions
|
