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 <optional>
10 : #include <string>
11 : #include <unordered_map>
12 : #include <variant>
13 : #include <vector>
14 :
15 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
16 : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
17 : #include "Domain/CoordinateMaps/Distribution.hpp"
18 : #include "Domain/Creators/DomainCreator.hpp"
19 : #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
20 : #include "Domain/Creators/TimeDependentOptions/Sphere.hpp"
21 : #include "Domain/Domain.hpp"
22 : #include "Domain/Structure/DirectionMap.hpp"
23 : #include "Options/Auto.hpp"
24 : #include "Options/Context.hpp"
25 : #include "Options/Options.hpp"
26 : #include "Options/ParseError.hpp"
27 : #include "Options/String.hpp"
28 : #include "Utilities/TMPL.hpp"
29 :
30 : /// \cond
31 : namespace domain {
32 : namespace CoordinateMaps {
33 : class Affine;
34 : class BulgedCube;
35 : class EquatorialCompression;
36 : class Equiangular;
37 : template <typename Map1, typename Map2, typename Map3>
38 : class ProductOf3Maps;
39 : template <size_t Dim>
40 : class Wedge;
41 : } // namespace CoordinateMaps
42 :
43 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
44 : class CoordinateMap;
45 : } // namespace domain
46 : /// \endcond
47 :
48 : namespace domain::creators::detail {
49 :
50 : /// Options for excising the interior of the sphere. This class parses as the
51 : /// `ExcisionFromOptions` subclass if boundary conditions are enabled, and as a
52 : /// plain string if boundary conditions are disabled.
53 : struct Excision {
54 : Excision() = default;
55 : Excision(std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
56 : boundary_condition);
57 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
58 : boundary_condition = nullptr;
59 : };
60 :
61 : struct ExcisionFromOptions : Excision {
62 : static constexpr Options::String help = {
63 : "Excise the interior of the sphere, leaving a spherical shell."};
64 : template <typename BoundaryConditionsBase>
65 : struct BoundaryCondition {
66 : static std::string name() { return "ExciseWithBoundaryCondition"; }
67 : using type = std::unique_ptr<BoundaryConditionsBase>;
68 : static constexpr Options::String help = {
69 : "The boundary condition to impose on the excision surface."};
70 : };
71 : template <typename Metavariables>
72 : using options = tmpl::list<BoundaryCondition<
73 : domain::BoundaryConditions::get_boundary_conditions_base<
74 : typename Metavariables::system>>>;
75 : using Excision::Excision;
76 : };
77 :
78 : /// Options for filling the interior of the sphere with a cube
79 : struct InnerCube {
80 : static constexpr Options::String help = {
81 : "Fill the interior of the sphere with a cube."};
82 : struct Sphericity {
83 : static std::string name() { return "FillWithSphericity"; }
84 : using type = double;
85 : static constexpr Options::String help = {
86 : "Sphericity of the inner cube. A sphericity of 0 uses a product "
87 : "of 1D maps as the map in the center. A sphericity > 0 uses a "
88 : "BulgedCube. A sphericity of exactly 1 is not allowed. See "
89 : "BulgedCube docs for why."};
90 : static double lower_bound() { return 0.0; }
91 : static double upper_bound() { return 1.0; }
92 : };
93 : using options = tmpl::list<Sphericity>;
94 : InnerCube() = default;
95 : explicit InnerCube(double sphericity_in) : sphericity(sphericity_in) {}
96 : double sphericity = std::numeric_limits<double>::signaling_NaN();
97 : };
98 :
99 : } // namespace domain::creators::detail
100 :
101 : template <>
102 : struct Options::create_from_yaml<domain::creators::detail::Excision> {
103 : template <typename Metavariables>
104 : static domain::creators::detail::Excision create(
105 : const Options::Option& options) {
106 : if constexpr (domain::BoundaryConditions::has_boundary_conditions_base_v<
107 : typename Metavariables::system>) {
108 : // Boundary conditions are enabled. Parse with a nested option.
109 : return options.parse_as<domain::creators::detail::ExcisionFromOptions,
110 : Metavariables>();
111 : } else {
112 : // Boundary conditions are disabled. Parse as a plain string.
113 : if (options.parse_as<std::string>() == "Excise") {
114 : return domain::creators::detail::Excision{};
115 : } else {
116 : PARSE_ERROR(options.context(), "Parse error");
117 : }
118 : }
119 : }
120 : };
121 :
122 : namespace domain::creators {
123 :
124 : /*!
125 : * \brief A 3D cubed sphere.
126 : *
127 : * Six wedges surround an interior region, which is either excised or filled in
128 : * with a seventh block. The interior region is a (possibly deformed) sphere
129 : * when excised, or a (possibly deformed) cube when filled in. Additional
130 : * spherical shells, each composed of six wedges, can be added with the
131 : * 'RadialPartitioning' option.
132 : *
133 : * \image html WedgeOrientations.png "The orientation of each wedge in a cubed
134 : * sphere."
135 : *
136 : * This domain creator offers one grid anchor "Center" at the origin.
137 : *
138 : * #### Inner cube sphericity
139 : * The inner cube is a BulgedCube except if the inner cube sphericity is
140 : * exactly 0. Then an Equiangular or Affine map is used (depending on if it's
141 : * equiangular or not) to avoid a root find in the BulgedCube map.
142 : *
143 : * #### Time dependent maps
144 : * There are two ways to add time dependent maps to the Sphere domain
145 : * creator. In the input file, these are specified under the
146 : * `TimeDependentMaps:` block.
147 : *
148 : * ##### TimeDependence
149 : * You can use a simple TimeDependence (e.g.
150 : * `domain::creators::time_dependence::UniformTranslation` or
151 : * `domain::creators::time_dependence::RotationAboutZAxis`) to add time
152 : * dependent maps. This method will add the same maps to all blocks in the
153 : * domain. This method can be used with an inner cube or with an excision
154 : * surface.
155 : *
156 : * ##### Hard-coded time dependent maps
157 : * The Sphere domain creator also has the option to use some hard coded time
158 : * dependent maps that may be useful in certain scenarios. This method adds the
159 : * maps in `domain::creators::sphere::TimeDependentMapOptions` to the domain.
160 : * Currently, the first (inner-most) shell has maps between `Frame::Grid`,
161 : * `Frame::Distorted`, and `Frame::Inertial` while all subsequent shells only
162 : * have maps between `Frame::Grid` and `Frame::Inertial`.
163 : *
164 : * \note You can only use hard-coded time dependent maps if you have an excision
165 : * surface. You cannot have a inner cube.
166 : *
167 : * ##### None
168 : * To not have any time dependent maps, pass a `std::nullopt` to appropriate
169 : * argument in the constructor. In the input file, simple have
170 : * `TimeDependentMaps: None`.
171 : *
172 : */
173 1 : class Sphere : public DomainCreator<3> {
174 : private:
175 0 : using Affine = CoordinateMaps::Affine;
176 0 : using Affine3D = CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
177 0 : using Equiangular = CoordinateMaps::Equiangular;
178 0 : using Equiangular3D =
179 : CoordinateMaps::ProductOf3Maps<Equiangular, Equiangular, Equiangular>;
180 0 : using BulgedCube = CoordinateMaps::BulgedCube;
181 :
182 : public:
183 0 : using maps_list = tmpl::append<
184 : tmpl::list<
185 : // Inner cube
186 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
187 : BulgedCube>,
188 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D>,
189 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
190 : Equiangular3D>,
191 : // Wedges
192 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
193 : CoordinateMaps::Wedge<3>>,
194 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
195 : CoordinateMaps::Wedge<3>,
196 : CoordinateMaps::EquatorialCompression>>,
197 : typename sphere::TimeDependentMapOptions::maps_list>;
198 :
199 0 : struct InnerRadius {
200 0 : using type = double;
201 0 : static constexpr Options::String help = {
202 : "Radius circumscribing the inner cube or the excision."};
203 : };
204 :
205 0 : struct OuterRadius {
206 0 : using type = double;
207 0 : static constexpr Options::String help = {"Radius of the sphere."};
208 : };
209 :
210 0 : using Excision = detail::Excision;
211 0 : using InnerCube = detail::InnerCube;
212 :
213 0 : struct Interior {
214 0 : using type = std::variant<Excision, InnerCube>;
215 0 : static constexpr Options::String help = {
216 : "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
217 : "excise the interior of the sphere, leaving a spherical shell "
218 : "(or just 'Excise' if boundary conditions are disabled). "
219 : "Or specify 'CubeWithSphericity' to fill the interior."};
220 : };
221 :
222 0 : struct InitialRefinement {
223 0 : using type =
224 : std::variant<size_t, std::array<size_t, 3>,
225 : std::vector<std::array<size_t, 3>>,
226 : std::unordered_map<std::string, std::array<size_t, 3>>>;
227 0 : static constexpr Options::String help = {
228 : "Initial refinement level. Specify one of: a single number, a "
229 : "list representing [phi, theta, r], or such a list for every block "
230 : "in the domain. The central cube always uses the value for 'theta' "
231 : "in both y- and z-direction."};
232 : };
233 :
234 0 : struct InitialGridPoints {
235 0 : using type =
236 : std::variant<size_t, std::array<size_t, 3>,
237 : std::vector<std::array<size_t, 3>>,
238 : std::unordered_map<std::string, std::array<size_t, 3>>>;
239 0 : static constexpr Options::String help = {
240 : "Initial number of grid points. Specify one of: a single number, a "
241 : "list representing [phi, theta, r], or such a list for every block "
242 : "in the domain. The central cube always uses the value for 'theta' "
243 : "in both y- and z-direction."};
244 : };
245 :
246 0 : struct UseEquiangularMap {
247 0 : using type = bool;
248 0 : static constexpr Options::String help = {
249 : "Use equiangular instead of equidistant coordinates. Equiangular "
250 : "coordinates give better gridpoint spacings in the angular "
251 : "directions, while equidistant coordinates give better gridpoint "
252 : "spacings in the inner cube."};
253 : };
254 :
255 : /// Options for the EquatorialCompression map
256 1 : struct EquatorialCompressionOptions {
257 0 : static constexpr Options::String help = {
258 : "Options for the EquatorialCompression map."};
259 0 : struct AspectRatio {
260 0 : using type = double;
261 0 : static constexpr Options::String help = {
262 : "An aspect ratio greater than 1 moves grid points toward the "
263 : "equator, and an aspect ratio smaller than 1 moves grid points "
264 : "toward the poles."};
265 0 : static double lower_bound() { return 0.0; }
266 : };
267 0 : struct IndexPolarAxis {
268 0 : using type = size_t;
269 0 : static constexpr Options::String help = {
270 : "The index (0, 1, or 2) of the axis along which equatorial "
271 : "compression is applied, where 0 is x, 1 is y, and 2 is z."};
272 0 : static size_t upper_bound() { return 2; }
273 : };
274 0 : using options = tmpl::list<AspectRatio, IndexPolarAxis>;
275 :
276 0 : double aspect_ratio;
277 0 : size_t index_polar_axis;
278 : };
279 :
280 0 : struct EquatorialCompression {
281 0 : using type =
282 : Options::Auto<EquatorialCompressionOptions, Options::AutoLabel::None>;
283 0 : static constexpr Options::String help = {
284 : "Apply an equatorial compression map to focus resolution on the "
285 : "equator or on the poles. The equatorial compression is an angular "
286 : "redistribution of grid points and will preserve the spherical shape "
287 : "of the inner and outer boundaries."};
288 : };
289 :
290 0 : struct RadialPartitioning {
291 0 : using type = std::vector<double>;
292 0 : static constexpr Options::String help = {
293 : "Radial coordinates of the boundaries splitting the spherical shell "
294 : "between InnerRadius and OuterRadius. They must be given in ascending "
295 : "order. This should be used if boundaries need to be set at specific "
296 : "radii. If the number but not the specific locations of the boundaries "
297 : "are important, use InitialRefinement instead."};
298 : };
299 :
300 0 : struct RadialDistribution {
301 0 : using type =
302 : std::variant<domain::CoordinateMaps::Distribution,
303 : std::vector<domain::CoordinateMaps::Distribution>>;
304 0 : static constexpr Options::String help = {
305 : "Select the radial distribution of grid points in each spherical "
306 : "shell. There must be N+1 radial distributions specified for N radial "
307 : "partitions. If the interior of the sphere is filled with a cube, the "
308 : "innermost shell must have a 'Linear' distribution because it changes "
309 : "in sphericity. You can also specify just a single radial distribution "
310 : "(not in a vector) which will use the same distribution for all "
311 : "partitions."};
312 : };
313 :
314 0 : struct WhichWedges {
315 0 : using type = ShellWedges;
316 0 : static constexpr Options::String help = {
317 : "Which wedges to include in the shell."};
318 0 : static constexpr type suggested_value() { return ShellWedges::All; }
319 : };
320 :
321 0 : using TimeDepOptionType = std::variant<
322 : sphere::TimeDependentMapOptions,
323 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>>;
324 :
325 0 : struct TimeDependentMaps {
326 0 : using type = Options::Auto<TimeDepOptionType, Options::AutoLabel::None>;
327 0 : static constexpr Options::String help = {
328 : "The options for time dependent maps. This can either be a "
329 : "TimeDependence or hard coded time dependent options. Specify `None` "
330 : "for no time dependent maps."};
331 : };
332 :
333 : template <typename BoundaryConditionsBase>
334 0 : struct OuterBoundaryCondition {
335 0 : static constexpr Options::String help =
336 : "Options for the boundary conditions at the outer radius.";
337 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
338 : };
339 :
340 0 : using basic_options =
341 : tmpl::list<InnerRadius, OuterRadius, Interior, InitialRefinement,
342 : InitialGridPoints, UseEquiangularMap, EquatorialCompression,
343 : RadialPartitioning, RadialDistribution, WhichWedges,
344 : TimeDependentMaps>;
345 :
346 : template <typename Metavariables>
347 0 : using options = tmpl::conditional_t<
348 : domain::BoundaryConditions::has_boundary_conditions_base_v<
349 : typename Metavariables::system>,
350 : tmpl::push_back<
351 : basic_options,
352 : OuterBoundaryCondition<
353 : domain::BoundaryConditions::get_boundary_conditions_base<
354 : typename Metavariables::system>>>,
355 : basic_options>;
356 :
357 0 : static constexpr Options::String help{
358 : "A 3D cubed sphere. Six wedges surround an interior region, which is "
359 : "either excised or filled in with a seventh block. The interior region "
360 : "is a (possibly deformed) sphere when excised, or a (possibly deformed) "
361 : "cube when filled in. Additional spherical shells, each composed of six "
362 : "wedges, can be added with the 'RadialPartitioning' option."};
363 :
364 0 : Sphere(
365 : double inner_radius, double outer_radius,
366 : std::variant<Excision, InnerCube> interior,
367 : const typename InitialRefinement::type& initial_refinement,
368 : const typename InitialGridPoints::type& initial_number_of_grid_points,
369 : bool use_equiangular_map,
370 : std::optional<EquatorialCompressionOptions> equatorial_compression = {},
371 : std::vector<double> radial_partitioning = {},
372 : const typename RadialDistribution::type& radial_distribution =
373 : domain::CoordinateMaps::Distribution::Linear,
374 : ShellWedges which_wedges = ShellWedges::All,
375 : std::optional<TimeDepOptionType> time_dependent_options = std::nullopt,
376 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
377 : outer_boundary_condition = nullptr,
378 : const Options::Context& context = {});
379 :
380 0 : Sphere() = default;
381 0 : Sphere(const Sphere&) = delete;
382 0 : Sphere(Sphere&&) = default;
383 0 : Sphere& operator=(const Sphere&) = delete;
384 0 : Sphere& operator=(Sphere&&) = default;
385 0 : ~Sphere() override = default;
386 :
387 0 : Domain<3> create_domain() const override;
388 :
389 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
390 1 : grid_anchors() const override {
391 : return grid_anchors_;
392 : }
393 :
394 : std::vector<DirectionMap<
395 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
396 1 : external_boundary_conditions() const override;
397 :
398 1 : std::vector<std::array<size_t, 3>> initial_extents() const override {
399 : return initial_number_of_grid_points_;
400 : }
401 :
402 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels()
403 : const override {
404 : return initial_refinement_;
405 : }
406 :
407 1 : std::vector<std::string> block_names() const override { return block_names_; }
408 :
409 : std::unordered_map<std::string, std::unordered_set<std::string>>
410 1 : block_groups() const override {
411 : return block_groups_;
412 : }
413 :
414 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
415 : initial_expiration_times = {}) const
416 : -> std::unordered_map<
417 : std::string,
418 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
419 :
420 : private:
421 0 : double inner_radius_{};
422 0 : double outer_radius_{};
423 0 : std::variant<Excision, InnerCube> interior_{};
424 0 : bool fill_interior_ = false;
425 0 : std::vector<std::array<size_t, 3>> initial_refinement_{};
426 0 : std::vector<std::array<size_t, 3>> initial_number_of_grid_points_{};
427 0 : bool use_equiangular_map_ = false;
428 0 : std::optional<EquatorialCompressionOptions> equatorial_compression_{};
429 0 : std::vector<double> radial_partitioning_{};
430 0 : std::vector<domain::CoordinateMaps::Distribution> radial_distribution_{};
431 0 : ShellWedges which_wedges_ = ShellWedges::All;
432 0 : std::optional<TimeDepOptionType> time_dependent_options_{};
433 0 : bool use_hard_coded_maps_{false};
434 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
435 0 : outer_boundary_condition_;
436 0 : size_t num_shells_{};
437 0 : size_t num_blocks_;
438 0 : size_t num_blocks_per_shell_;
439 0 : std::vector<std::string> block_names_{};
440 : std::unordered_map<std::string, std::unordered_set<std::string>>
441 0 : block_groups_{};
442 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
443 0 : grid_anchors_{};
444 : };
445 :
446 : } // namespace domain::creators
|
