Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <array>
7 : #include <cstddef>
8 : #include <memory>
9 : #include <string>
10 : #include <unordered_map>
11 : #include <unordered_set>
12 : #include <vector>
13 :
14 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
15 : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
16 : #include "Domain/Creators/DomainCreator.hpp"
17 : #include "Domain/Creators/Sphere.hpp"
18 : #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
19 : #include "Domain/Domain.hpp"
20 : #include "Domain/Structure/DirectionMap.hpp"
21 : #include "Options/Context.hpp"
22 : #include "Options/String.hpp"
23 : #include "Utilities/TMPL.hpp"
24 :
25 : /// \cond
26 : namespace domain {
27 : namespace CoordinateMaps {
28 : class Affine;
29 : class Equiangular;
30 : template <size_t Dim>
31 : class Identity;
32 : template <typename Map1, typename Map2>
33 : class ProductOf2Maps;
34 : template <typename Map1, typename Map2, typename Map3>
35 : class ProductOf3Maps;
36 : template <size_t Dim>
37 : class Wedge;
38 : template <size_t Dim>
39 : class DiscreteRotation;
40 : } // namespace CoordinateMaps
41 :
42 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
43 : class CoordinateMap;
44 : } // namespace domain
45 : /// \endcond
46 :
47 : namespace domain::creators::detail{
48 :
49 : /// Options for filling the interior of the sphere with a cube
50 : struct InnerSquare {
51 : static constexpr Options::String help = {
52 : "Fill the interior of the 2D sphere with a half-square."};
53 : struct Sphericity {
54 : static std::string name() { return "FillWithSphericity"; }
55 : using type = double;
56 : static constexpr Options::String help = {
57 : "Sphericity of the inner half-square. Only a sphericity of "
58 : "0.0 is currently implemented."};
59 : static double lower_bound() { return 0.0; }
60 : static double upper_bound() { return 0.0; }
61 : };
62 : using options = tmpl::list<Sphericity>;
63 : InnerSquare() = default;
64 : explicit InnerSquare(double sphericity_in) : sphericity(sphericity_in) {}
65 : double sphericity = std::numeric_limits<double>::signaling_NaN();
66 : };
67 : } // namespace domain::creators::detail
68 :
69 : namespace domain::creators {
70 : /// Create a 3D Domain with a half-disk computational domain employing axial
71 : /// symmetry. The third dimension uses a Cartoon basis with Killing vector
72 : /// along the $\phi$ direction.
73 1 : class CartoonSphere2D : public DomainCreator<3> {
74 : public:
75 0 : using maps_list = tmpl::list<
76 : domain::CoordinateMap<
77 : Frame::BlockLogical, Frame::Inertial,
78 : CoordinateMaps::ProductOf2Maps<
79 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Affine,
80 : CoordinateMaps::Affine>,
81 : CoordinateMaps::Identity<1>>>,
82 : domain::CoordinateMap<
83 : Frame::BlockLogical, Frame::Inertial,
84 : CoordinateMaps::ProductOf2Maps<
85 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Equiangular,
86 : CoordinateMaps::Equiangular>,
87 : CoordinateMaps::Identity<1>>>,
88 : domain::CoordinateMap<
89 : Frame::BlockLogical, Frame::Inertial,
90 : CoordinateMaps::DiscreteRotation<3>,
91 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Wedge<2>,
92 : CoordinateMaps::Identity<1>>>>;
93 :
94 0 : struct InnerRadius {
95 0 : using type = double;
96 0 : static constexpr Options::String help = {
97 : "Radius of the circle circumscribing the inner half-square, or the "
98 : "radius of the inner boundary if ExciseCenter = true."};
99 : };
100 :
101 0 : struct OuterRadius {
102 0 : using type = double;
103 0 : static constexpr Options::String help = {"Outer radius of the half-disk."};
104 : };
105 :
106 0 : struct InitialRefinement {
107 0 : using type =
108 : std::variant<std::array<size_t, 2>, std::vector<std::array<size_t, 2>>>;
109 0 : static constexpr Options::String help = {
110 : "Initial refinement level in [r, theta]. If one pair is given, it will "
111 : "be applied to all blocks, otherwise each shell can be specified.\n"
112 : "Note: the half-wedges will have their theta values decremented by "
113 : "one.\n"
114 : "Note: the inner square, if included, will have refinement set to "
115 : "the theta value for the center half-circle for both dimensions "
116 : "(decremented by one in the halved dimension)."};
117 : };
118 :
119 0 : struct InitialGridPoints {
120 0 : using type =
121 : std::variant<std::array<size_t, 2>, std::vector<std::array<size_t, 2>>>;
122 0 : static constexpr Options::String help = {
123 : "Initial number of grid points in [r,theta]. If one pair is given, it "
124 : "will be applied to all blocks, otherwise each shell can be "
125 : "specified, from innermost to outermost.\n"
126 : "Note: if included, the inner square will have both dimensions'"
127 : "number of grid points set to the theta value of the surrounding "
128 : "wedges."};
129 : };
130 :
131 0 : struct RadialPartitioning {
132 0 : using type = std::vector<double>;
133 0 : static constexpr Options::String help = {
134 : "Radial coordinates of the boundaries splitting the radial shells. "
135 : "Leave emtpy for only one layer."};
136 : };
137 :
138 0 : struct UseEquiangularMap {
139 0 : using type = bool;
140 0 : static constexpr Options::String help = {
141 : "Use equiangular instead of equidistant coordinates in wedges."};
142 : };
143 :
144 0 : using Excision = detail::Excision;
145 0 : using InnerSquare = detail::InnerSquare;
146 :
147 0 : struct Interior {
148 0 : using type = std::variant<Excision, InnerSquare>;
149 0 : static constexpr Options::String help = {
150 : "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
151 : "excise the interior of the sphere, leaving a spherical shell "
152 : "(or just 'Excise' if boundary conditions are disabled). "
153 : "Or specify 'FillWithSphericity' to fill the interior."};
154 : };
155 :
156 : template <typename BoundaryConditionsBase>
157 0 : struct YAxisBoundaryCondition {
158 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
159 0 : static constexpr Options::String help = {
160 : "The boundary condition to impose at the x=z=0 boundary."};
161 : };
162 :
163 : template <typename BoundaryConditionsBase>
164 0 : struct OuterBoundaryCondition {
165 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
166 0 : static constexpr Options::String help = {
167 : "The boundary condition to impose at the outer boundary of the "
168 : "domain."};
169 : };
170 :
171 0 : struct TimeDependence {
172 0 : using type =
173 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>;
174 0 : static constexpr Options::String help = {
175 : "The time dependence of the moving mesh domain. Specify `None` for no "
176 : "time dependant maps."};
177 : };
178 :
179 0 : using basic_options =
180 : tmpl::list<InnerRadius, OuterRadius, InitialRefinement, InitialGridPoints,
181 : RadialPartitioning, UseEquiangularMap, Interior,
182 : TimeDependence>;
183 :
184 : template <typename Metavariables>
185 0 : using options = tmpl::conditional_t<
186 : domain::BoundaryConditions::has_boundary_conditions_base_v<
187 : typename Metavariables::system>,
188 : tmpl::push_back<
189 : basic_options,
190 : YAxisBoundaryCondition<
191 : domain::BoundaryConditions::get_boundary_conditions_base<
192 : typename Metavariables::system>>,
193 : OuterBoundaryCondition<
194 : domain::BoundaryConditions::get_boundary_conditions_base<
195 : typename Metavariables::system>>>,
196 : basic_options>;
197 :
198 0 : static constexpr Options::String help{
199 : "Creates a sphere that requires/enforces axial-symmetry.\n"
200 : "The computational domain is a 2D half-disk, consisting of an inner\n"
201 : "half-square with two half-wedges above and below and a full wedge to\n"
202 : "the right. This can be extended to have mulitple shells of wedges by\n"
203 : "specifying a radial partition. The inner half-square can be excised,\n"
204 : "in which case the inner sphericity is set to 1.\n"
205 : "Equiangular coordinates give better gridpoint spacings in the angular\n"
206 : "direction, while equidistant coordinates give better gridpoint\n"
207 : "spacings in the center half-square."};
208 :
209 0 : CartoonSphere2D(
210 : double inner_radius, double outer_radius,
211 : typename InitialRefinement::type&& initial_refinement,
212 : typename InitialGridPoints::type&& initial_number_of_grid_points,
213 : std::vector<double> radial_partitioning, bool use_equiangular_map,
214 : std::variant<Excision, InnerSquare> interior,
215 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
216 : time_dependence = nullptr,
217 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
218 : y_axis_boundary_condition = nullptr,
219 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
220 : outer_boundary_condition = nullptr,
221 : const Options::Context& context = {});
222 :
223 0 : CartoonSphere2D() = default;
224 0 : CartoonSphere2D(const CartoonSphere2D&) = delete;
225 0 : CartoonSphere2D(CartoonSphere2D&&) = default;
226 0 : CartoonSphere2D& operator=(const CartoonSphere2D&) = delete;
227 0 : CartoonSphere2D& operator=(CartoonSphere2D&&) = default;
228 0 : ~CartoonSphere2D() override = default;
229 :
230 0 : Domain<3> create_domain() const override;
231 :
232 : std::vector<DirectionMap<
233 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
234 1 : external_boundary_conditions() const override;
235 :
236 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
237 :
238 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
239 :
240 1 : std::vector<std::string> block_names() const override { return block_names_; }
241 :
242 : std::unordered_map<std::string, std::unordered_set<std::string>>
243 1 : block_groups() const override {
244 : return block_groups_;
245 : }
246 :
247 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
248 : initial_expiration_times = {}) const
249 : -> std::unordered_map<
250 : std::string,
251 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
252 :
253 : private:
254 0 : double inner_radius_{};
255 0 : double outer_radius_{};
256 0 : std::vector<std::array<size_t, 2>> initial_refinement_{};
257 0 : std::vector<std::array<size_t, 2>> initial_number_of_grid_points_{};
258 0 : std::vector<double> radial_partitioning_{};
259 0 : bool use_equiangular_map_{false};
260 0 : std::variant<Excision, InnerSquare> interior_{};
261 0 : bool fill_interior_{false};
262 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
263 0 : time_dependence_ = nullptr;
264 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
265 0 : y_axis_boundary_condition_ = nullptr;
266 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
267 0 : outer_boundary_condition_ = nullptr;
268 0 : std::vector<std::string> block_names_;
269 : std::unordered_map<std::string, std::unordered_set<std::string>>
270 0 : block_groups_;
271 0 : size_t num_shells_{};
272 0 : size_t num_blocks_{};
273 : };
274 : } // namespace domain::creators
|
