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 <variant>
13 : #include <vector>
14 :
15 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
16 : #include "Domain/BoundaryConditions/Cartoon.hpp"
17 : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
18 : #include "Domain/CoordinateMaps/Distribution.hpp"
19 : #include "Domain/Creators/DomainCreator.hpp"
20 : #include "Domain/Creators/Sphere.hpp"
21 : #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
22 : #include "Domain/Domain.hpp"
23 : #include "Domain/Structure/DirectionMap.hpp"
24 : #include "Options/Context.hpp"
25 : #include "Options/String.hpp"
26 : #include "Utilities/TMPL.hpp"
27 :
28 : /// \cond
29 : namespace domain {
30 : namespace CoordinateMaps {
31 : class Affine;
32 : class Equiangular;
33 : template <size_t Dim>
34 : class Identity;
35 : template <typename Map1, typename Map2>
36 : class ProductOf2Maps;
37 : template <typename Map1, typename Map2, typename Map3>
38 : class ProductOf3Maps;
39 : template <size_t Dim>
40 : class Wedge;
41 : template <size_t Dim>
42 : class DiscreteRotation;
43 : } // namespace CoordinateMaps
44 :
45 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
46 : class CoordinateMap;
47 : } // namespace domain
48 : /// \endcond
49 :
50 : namespace domain::creators::detail {
51 : /// Options for filling the interior of the disk with a square
52 : struct InnerSquare {
53 : static constexpr Options::String help = {
54 : "Fill the interior of the 2D sphere with a half-square."};
55 : struct Sphericity {
56 : static std::string name() { return "FillWithSphericity"; }
57 : using type = double;
58 : static constexpr Options::String help = {
59 : "Sphericity of the inner half-square. Only a sphericity of "
60 : "0.0 is currently implemented."};
61 : static double lower_bound() { return 0.0; }
62 : static double upper_bound() { return 0.0; }
63 : };
64 : using options = tmpl::list<Sphericity>;
65 : InnerSquare() = default;
66 : explicit InnerSquare(double sphericity_in) : sphericity(sphericity_in) {}
67 : double sphericity = std::numeric_limits<double>::signaling_NaN();
68 : };
69 : } // namespace domain::creators::detail
70 :
71 : namespace domain::creators {
72 : /// Create a 3D Domain with a half-disk computational domain employing axial
73 : /// symmetry. The third dimension uses a Cartoon basis with Killing vector
74 : /// along the $\phi$ direction.
75 1 : class CartoonSphere2D : public DomainCreator<3> {
76 : public:
77 0 : using maps_list = tmpl::list<
78 : domain::CoordinateMap<
79 : Frame::BlockLogical, Frame::Inertial,
80 : CoordinateMaps::ProductOf2Maps<
81 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Affine,
82 : CoordinateMaps::Affine>,
83 : CoordinateMaps::Identity<1>>>,
84 : domain::CoordinateMap<
85 : Frame::BlockLogical, Frame::Inertial,
86 : CoordinateMaps::ProductOf2Maps<
87 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Equiangular,
88 : CoordinateMaps::Equiangular>,
89 : CoordinateMaps::Identity<1>>>,
90 : domain::CoordinateMap<
91 : Frame::BlockLogical, Frame::Inertial,
92 : CoordinateMaps::DiscreteRotation<3>,
93 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Wedge<2>,
94 : CoordinateMaps::Identity<1>>>>;
95 :
96 0 : struct InnerRadius {
97 0 : using type = double;
98 0 : static constexpr Options::String help = {
99 : "Radius of the circle circumscribing the inner half-square, or the "
100 : "radius of the inner boundary if ExciseCenter = true."};
101 : };
102 :
103 0 : struct OuterRadius {
104 0 : using type = double;
105 0 : static constexpr Options::String help = {"Outer radius of the half-disk."};
106 : };
107 :
108 0 : struct InitialAngularRefinement {
109 0 : using type = size_t;
110 0 : static constexpr Options::String help = {
111 : "Initial refinement level in theta. It will be applied to all "
112 : "blocks because h-refinement is not supported for the ZernikeB1 "
113 : "basis.\n"
114 : "Note: the half-wedges will have their theta values decremented by "
115 : "one.\n"
116 : "Note: the inner square, if included, will have refinement set to "
117 : "the theta value for the center half-circle for both dimensions "
118 : "(decremented by one in the halved dimension)."};
119 : };
120 :
121 0 : struct InitialRadialRefinement {
122 0 : using type = std::variant<size_t, std::vector<size_t>>;
123 0 : static constexpr Options::String help = {
124 : "Initial refinement level in r. If one value is given, it will be "
125 : "applied to all blocks, otherwise each shell can be specificed, "
126 : "from innermost to outermost.\n"};
127 : };
128 :
129 0 : struct InitialGridPoints {
130 0 : using type = std::array<size_t, 2>;
131 0 : static constexpr Options::String help = {
132 : "Initial number of grid points in [r,theta]. It will be applied to all "
133 : "blocks because p-refinement is not supported for the ZernikeB1 "
134 : "basis.\n"
135 : "Note: if included, the inner square will have both dimensions'"
136 : "number of grid points set to the theta value of the surrounding "
137 : "wedges."};
138 : };
139 :
140 0 : struct RadialPartitioning {
141 0 : using type = std::vector<double>;
142 0 : static constexpr Options::String help = {
143 : "Radial coordinates of the boundaries splitting the radial shells. "
144 : "Leave empty for only one layer."};
145 : };
146 :
147 0 : struct UseEquiangularMap {
148 0 : using type = bool;
149 0 : static constexpr Options::String help = {
150 : "Use equiangular instead of equidistant coordinates in wedges."};
151 : };
152 :
153 0 : using Excision = detail::Excision;
154 0 : using InnerSquare = detail::InnerSquare;
155 :
156 0 : struct Interior {
157 0 : using type = std::variant<Excision, InnerSquare>;
158 0 : static constexpr Options::String help = {
159 : "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
160 : "excise the interior of the sphere, leaving a spherical shell "
161 : "(or just 'Excise' if boundary conditions are disabled). "
162 : "Or specify 'FillWithSphericity' to fill the interior."};
163 : };
164 :
165 : template <typename BoundaryConditionsBase>
166 0 : struct OuterBoundaryCondition {
167 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
168 0 : static constexpr Options::String help = {
169 : "The boundary condition to impose at the outer boundary of the "
170 : "domain."};
171 : };
172 :
173 0 : struct RadialDistribution {
174 0 : using type = std::variant<CoordinateMaps::Distribution,
175 : std::vector<CoordinateMaps::Distribution>>;
176 0 : static constexpr Options::String help = {
177 : "Distribution of grid points for the radial direction of the wedges. "
178 : "A single value will be applied to all shells, or every shell can be "
179 : "specified individually, in which case for N partitions there must be "
180 : "N+1 distributions. For anything but linear in the innermost shell, "
181 : "the center must be excised."};
182 : };
183 :
184 0 : struct TimeDependence {
185 0 : using type =
186 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>;
187 0 : static constexpr Options::String help = {
188 : "The time dependence of the moving mesh domain. Specify `None` for no "
189 : "time dependant maps."};
190 : };
191 :
192 0 : using basic_options =
193 : tmpl::list<InnerRadius, OuterRadius, InitialAngularRefinement,
194 : InitialRadialRefinement, InitialGridPoints, RadialPartitioning,
195 : UseEquiangularMap, Interior, RadialDistribution,
196 : TimeDependence>;
197 :
198 : template <typename Metavariables>
199 0 : using options = tmpl::conditional_t<
200 : domain::BoundaryConditions::has_boundary_conditions_base_v<
201 : typename Metavariables::system>,
202 : tmpl::push_back<
203 : basic_options,
204 : OuterBoundaryCondition<
205 : domain::BoundaryConditions::get_boundary_conditions_base<
206 : typename Metavariables::system>>>,
207 : basic_options>;
208 :
209 0 : static constexpr Options::String help{
210 : "Creates a sphere that requires/enforces axial-symmetry.\n"
211 : "The computational domain is a 2D half-disk, consisting of an inner\n"
212 : "half-square with two half-wedges above and below and a full wedge to\n"
213 : "the right. This can be extended to have mulitple shells of wedges by\n"
214 : "specifying a radial partition. The inner half-square can be excised,\n"
215 : "in which case the inner sphericity is set to 1.\n"
216 : "Equiangular coordinates give better gridpoint spacings in the angular\n"
217 : "direction, while equidistant coordinates give better gridpoint\n"
218 : "spacings in the center half-square.\n"
219 : "Elements touching the x=0 axis use ZernikeB1 bases in the x direction\n"
220 : "and automatically apply a system's Cartoon-type boundary condition.\n"
221 : "Angular refinement must be globally set as the required mortar\n"
222 : "projection has not been implemented."};
223 :
224 0 : CartoonSphere2D(
225 : double inner_radius, double outer_radius,
226 : size_t initial_angular_refinement,
227 : const typename InitialRadialRefinement::type& initial_radial_refinement,
228 : std::array<size_t, 2> initial_number_of_grid_points,
229 : std::vector<double> radial_partitioning, bool use_equiangular_map,
230 : std::variant<Excision, InnerSquare> interior,
231 : const typename RadialDistribution::type& radial_distribution,
232 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
233 : time_dependence,
234 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
235 : outer_boundary_condition,
236 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
237 : cartoon_boundary_condition,
238 : const Options::Context& context);
239 :
240 0 : CartoonSphere2D() = default;
241 0 : CartoonSphere2D(const CartoonSphere2D&) = delete;
242 0 : CartoonSphere2D(CartoonSphere2D&&) = default;
243 0 : CartoonSphere2D& operator=(const CartoonSphere2D&) = delete;
244 0 : CartoonSphere2D& operator=(CartoonSphere2D&&) = default;
245 0 : ~CartoonSphere2D() override = default;
246 :
247 0 : Domain<3> create_domain() const override;
248 :
249 : std::vector<DirectionMap<
250 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
251 1 : external_boundary_conditions() const override;
252 :
253 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
254 :
255 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
256 :
257 : /// The block names are Shell0_LowerY, Shell0_UpperX, Shell0_UpperY,
258 : /// Shell1_LowerY, and so on. The naming and numbering goes from outer-most
259 : /// shell, starting from bottom and going counterclockwise, going to inside
260 : /// neighboring shell. If the center is included, the "center" half-circle
261 : /// follows same numbering with the _HalfSquare being the last block.
262 1 : std::vector<std::string> block_names() const override { return block_names_; }
263 :
264 : /// The block groups are Shell0, Shell1, ... starting from the outermost
265 : /// partition and working in.
266 : std::unordered_map<std::string, std::unordered_set<std::string>>
267 1 : block_groups() const override {
268 : return block_groups_;
269 : }
270 :
271 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
272 : initial_expiration_times = {}) const
273 : -> std::unordered_map<
274 : std::string,
275 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
276 :
277 : private:
278 0 : double inner_radius_{};
279 0 : double outer_radius_{};
280 0 : size_t initial_angular_refinement_{};
281 0 : std::vector<size_t> initial_radial_refinement_{};
282 0 : std::array<size_t, 2> initial_number_of_grid_points_{};
283 0 : std::vector<double> radial_partitioning_{};
284 0 : bool use_equiangular_map_{false};
285 0 : std::variant<Excision, InnerSquare> interior_{};
286 0 : bool fill_interior_{false};
287 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
288 0 : time_dependence_ = nullptr;
289 0 : std::vector<CoordinateMaps::Distribution> radial_distributions_{};
290 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
291 0 : outer_boundary_condition_ = nullptr;
292 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
293 0 : cartoon_boundary_condition_ = nullptr;
294 0 : std::vector<std::string> block_names_;
295 : std::unordered_map<std::string, std::unordered_set<std::string>>
296 0 : block_groups_;
297 0 : size_t num_shells_{};
298 0 : size_t num_blocks_{};
299 : };
300 : } // namespace domain::creators
301 :
302 : namespace domain::creators::detail {
303 : /// \brief Helper struct for CartoonSphere2D options parsing so the internal
304 : /// cartoon boundary condition at $x = 0$ is not a required argument.
305 : ///
306 : /// \details To get the cartoon-type boundary condition from a system we need
307 : /// access to the system's metavariables, which is only present when parsing
308 : /// options. The point of this design is so that the input file does not require
309 : /// a dummy value for an InnerBoundaryCondition that is always the same for a
310 : /// given system.
311 : ///
312 : /// This helper's constructor does not take an InnerBoundaryCondition, which
313 : /// means if we parse a CartoonSphere2D as a CartonSphere2DOptionsHelper
314 : /// by having an options parsing specialization (so the input file values are
315 : /// the helper's construtor, not CartoonSphere2D's), we can detect the
316 : /// cartoon-type boundary condition for the given system and only then call the
317 : /// CartoonSphere2D constructor with the extra information.
318 : struct CartoonSphere2DOptionsHelper {
319 : // Inherit the options template from CartoonSphere2D
320 : template <typename Metavariables>
321 : using options = typename domain::creators::CartoonSphere2D::template options<
322 : Metavariables>;
323 :
324 : using InitialRadialRefinement =
325 : domain::creators::CartoonSphere2D::InitialRadialRefinement;
326 : using Interior = domain::creators::CartoonSphere2D::Interior;
327 :
328 : template <typename BoundaryConditionsBase>
329 : using OuterBoundaryCondition =
330 : domain::creators::CartoonSphere2D::OuterBoundaryCondition<
331 : BoundaryConditionsBase>;
332 :
333 : static constexpr Options::String help = {"CartoonSphere2DOptionsHelper."};
334 :
335 : // Default constructor required by Options system
336 : CartoonSphere2DOptionsHelper() = default;
337 :
338 : // Does not take cartoon BC
339 : CartoonSphere2DOptionsHelper(
340 : double inner_radius, double outer_radius,
341 : size_t initial_angular_refinement,
342 : typename InitialRadialRefinement::type&& initial_radial_refinement,
343 : std::array<size_t, 2> initial_number_of_grid_points,
344 : std::vector<double> radial_partitioning, bool use_equiangular_map,
345 : std::variant<Excision, InnerSquare> interior,
346 : typename domain::creators::CartoonSphere2D::RadialDistribution::type
347 : radial_distribution = CoordinateMaps::Distribution::Linear,
348 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
349 : time_dependence = nullptr,
350 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
351 : outer_boundary_condition = nullptr,
352 : Options::Context&& context = {});
353 :
354 : // Same members as in CartoonSphere2D; public, to be extracted
355 : double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
356 : double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
357 : size_t initial_angular_refinement_{0};
358 : typename InitialRadialRefinement::type initial_radial_refinement_{};
359 : std::array<size_t, 2> initial_number_of_grid_points_{};
360 : std::vector<double> radial_partitioning_{};
361 : bool use_equiangular_map_{false};
362 : typename Interior::type interior_{};
363 : typename domain::creators::CartoonSphere2D::RadialDistribution::type
364 : radial_distribution_{CoordinateMaps::Distribution::Linear};
365 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
366 : time_dependence_{nullptr};
367 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
368 : outer_boundary_condition_{nullptr};
369 : Options::Context context_{};
370 : };
371 : } // namespace domain::creators::detail
372 :
373 : // Options parsing specialization to automate CartoonSphere2D's cartoon
374 : // boundary condition
375 : template <>
376 0 : struct Options::create_from_yaml<domain::creators::CartoonSphere2D> {
377 : template <typename Metavariables>
378 0 : static domain::creators::CartoonSphere2D create(
379 : const Options::Option& options) {
380 : auto helper =
381 : options.parse_as<domain::creators::detail::CartoonSphere2DOptionsHelper,
382 : Metavariables>();
383 :
384 : // Create cartoon BC if system supports it, if not will throw parse error in
385 : // real constructor
386 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
387 : cartoon_boundary_condition = nullptr;
388 : if constexpr (domain::BoundaryConditions::has_boundary_conditions_base_v<
389 : typename Metavariables::system>) {
390 : if constexpr (domain::BoundaryConditions::system_has_cartoon_bc_v<
391 : Metavariables>) {
392 : cartoon_boundary_condition =
393 : domain::BoundaryConditions::make_cartoon_boundary_condition<
394 : Metavariables>();
395 : }
396 : }
397 :
398 : // Construct CartoonSphere2D using the helper's parsed data + cartoon BC
399 : return domain::creators::CartoonSphere2D(
400 : helper.inner_radius_, helper.outer_radius_,
401 : helper.initial_angular_refinement_,
402 : std::move(helper.initial_radial_refinement_),
403 : std::move(helper.initial_number_of_grid_points_),
404 : std::move(helper.radial_partitioning_), helper.use_equiangular_map_,
405 : std::move(helper.interior_), helper.radial_distribution_,
406 : std::move(helper.time_dependence_),
407 : std::move(helper.outer_boundary_condition_),
408 : std::move(cartoon_boundary_condition), helper.context_);
409 : }
410 : };
|
