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 : * When a shape map is specified, an excised sphere uses it on the innermost
161 : * radial shell by default. A filled sphere with multiple radial shells uses it
162 : * on at least the two innermost shells: the first shell turns on the
163 : * deformation away from the central cube, and the subsequent selected shells
164 : * roll it off. At least one outer shell must remain without a shape map, so
165 : * this filled configuration requires at least three radial shells. The
166 : * `NumberOfRadialShellsWithShapeMap` option can be used to choose more inner
167 : * shells. With only one radial shell, a filled sphere does not roll off the
168 : * shape map, while an excised sphere rolls it off at the outer boundary.
169 : *
170 : * ##### None
171 : * To not have any time dependent maps, pass a `std::nullopt` to the appropriate
172 : * argument in the constructor. In the input file, simply have
173 : * `TimeDependentMaps: None`.
174 : *
175 : */
176 1 : class Sphere : public DomainCreator<3> {
177 : private:
178 0 : using Affine = CoordinateMaps::Affine;
179 0 : using Affine3D = CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
180 0 : using Equiangular = CoordinateMaps::Equiangular;
181 0 : using Equiangular3D =
182 : CoordinateMaps::ProductOf3Maps<Equiangular, Equiangular, Equiangular>;
183 0 : using BulgedCube = CoordinateMaps::BulgedCube;
184 :
185 : public:
186 0 : using maps_list = tmpl::append<
187 : tmpl::list<
188 : // Inner cube
189 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
190 : BulgedCube>,
191 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D>,
192 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
193 : Equiangular3D>,
194 : // Wedges
195 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
196 : CoordinateMaps::Wedge<3>>,
197 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
198 : CoordinateMaps::Wedge<3>,
199 : CoordinateMaps::EquatorialCompression>>,
200 : typename sphere::TimeDependentMapOptions::maps_list>;
201 :
202 0 : struct InnerRadius {
203 0 : using type = double;
204 0 : static constexpr Options::String help = {
205 : "Radius circumscribing the inner cube or the excision."};
206 : };
207 :
208 0 : struct OuterRadius {
209 0 : using type = double;
210 0 : static constexpr Options::String help = {"Radius of the sphere."};
211 : };
212 :
213 0 : using Excision = detail::Excision;
214 0 : using InnerCube = detail::InnerCube;
215 :
216 0 : struct Interior {
217 0 : using type = std::variant<Excision, InnerCube>;
218 0 : static constexpr Options::String help = {
219 : "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
220 : "excise the interior of the sphere, leaving a spherical shell "
221 : "(or just 'Excise' if boundary conditions are disabled). "
222 : "Or specify 'CubeWithSphericity' to fill the interior."};
223 : };
224 :
225 0 : struct InitialRefinement {
226 0 : using type =
227 : std::variant<size_t, std::array<size_t, 3>,
228 : std::vector<std::array<size_t, 3>>,
229 : std::unordered_map<std::string, std::array<size_t, 3>>>;
230 0 : static constexpr Options::String help = {
231 : "Initial refinement level. Specify one of: a single number, a "
232 : "list representing [phi, theta, r], or such a list for every block "
233 : "in the domain. The central cube always uses the value for 'theta' "
234 : "in both y- and z-direction."};
235 : };
236 :
237 0 : struct InitialGridPoints {
238 0 : using type =
239 : std::variant<size_t, std::array<size_t, 3>,
240 : std::vector<std::array<size_t, 3>>,
241 : std::unordered_map<std::string, std::array<size_t, 3>>>;
242 0 : static constexpr Options::String help = {
243 : "Initial number of grid points. Specify one of: a single number, a "
244 : "list representing [phi, theta, r], or such a list for every block "
245 : "in the domain. The central cube always uses the value for 'theta' "
246 : "in both y- and z-direction."};
247 : };
248 :
249 0 : struct UseEquiangularMap {
250 0 : using type = bool;
251 0 : static constexpr Options::String help = {
252 : "Use equiangular instead of equidistant coordinates. Equiangular "
253 : "coordinates give better gridpoint spacings in the angular "
254 : "directions, while equidistant coordinates give better gridpoint "
255 : "spacings in the inner cube."};
256 : };
257 :
258 : /// Options for the EquatorialCompression map
259 1 : struct EquatorialCompressionOptions {
260 0 : static constexpr Options::String help = {
261 : "Options for the EquatorialCompression map."};
262 0 : struct AspectRatio {
263 0 : using type = double;
264 0 : static constexpr Options::String help = {
265 : "An aspect ratio greater than 1 moves grid points toward the "
266 : "equator, and an aspect ratio smaller than 1 moves grid points "
267 : "toward the poles."};
268 0 : static double lower_bound() { return 0.0; }
269 : };
270 0 : struct IndexPolarAxis {
271 0 : using type = size_t;
272 0 : static constexpr Options::String help = {
273 : "The index (0, 1, or 2) of the axis along which equatorial "
274 : "compression is applied, where 0 is x, 1 is y, and 2 is z."};
275 0 : static size_t upper_bound() { return 2; }
276 : };
277 0 : using options = tmpl::list<AspectRatio, IndexPolarAxis>;
278 :
279 0 : double aspect_ratio;
280 0 : size_t index_polar_axis;
281 : };
282 :
283 0 : struct EquatorialCompression {
284 0 : using type =
285 : Options::Auto<EquatorialCompressionOptions, Options::AutoLabel::None>;
286 0 : static constexpr Options::String help = {
287 : "Apply an equatorial compression map to focus resolution on the "
288 : "equator or on the poles. The equatorial compression is an angular "
289 : "redistribution of grid points and will preserve the spherical shape "
290 : "of the inner and outer boundaries."};
291 : };
292 :
293 0 : struct RadialPartitioning {
294 0 : using type = std::vector<double>;
295 0 : static constexpr Options::String help = {
296 : "Radial coordinates of the boundaries splitting the spherical shell "
297 : "between InnerRadius and OuterRadius. They must be given in ascending "
298 : "order. This should be used if boundaries need to be set at specific "
299 : "radii. If the number but not the specific locations of the boundaries "
300 : "are important, use InitialRefinement instead."};
301 : };
302 :
303 0 : struct RadialDistribution {
304 0 : using type =
305 : std::variant<domain::CoordinateMaps::Distribution,
306 : std::vector<domain::CoordinateMaps::Distribution>>;
307 0 : static constexpr Options::String help = {
308 : "Select the radial distribution of grid points in each spherical "
309 : "shell. There must be N+1 radial distributions specified for N radial "
310 : "partitions. If the interior of the sphere is filled with a cube, the "
311 : "innermost shell must have a 'Linear' distribution because it changes "
312 : "in sphericity. You can also specify just a single radial distribution "
313 : "(not in a vector) which will use the same distribution for all "
314 : "partitions."};
315 : };
316 :
317 0 : struct WhichWedges {
318 0 : using type = ShellWedges;
319 0 : static constexpr Options::String help = {
320 : "Which wedges to include in the shell."};
321 0 : static constexpr type suggested_value() { return ShellWedges::All; }
322 : };
323 :
324 0 : using TimeDepOptionType = std::variant<
325 : sphere::TimeDependentMapOptions,
326 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>>;
327 :
328 0 : struct TimeDependentMaps {
329 0 : using type = Options::Auto<TimeDepOptionType, Options::AutoLabel::None>;
330 0 : static constexpr Options::String help = {
331 : "The options for time dependent maps. This can either be a "
332 : "TimeDependence or hard coded time dependent options. Specify `None` "
333 : "for no time dependent maps."};
334 : };
335 :
336 : template <typename BoundaryConditionsBase>
337 0 : struct OuterBoundaryCondition {
338 0 : static constexpr Options::String help =
339 : "Options for the boundary conditions at the outer radius.";
340 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
341 : };
342 :
343 0 : using basic_options =
344 : tmpl::list<InnerRadius, OuterRadius, Interior, InitialRefinement,
345 : InitialGridPoints, UseEquiangularMap, EquatorialCompression,
346 : RadialPartitioning, RadialDistribution, WhichWedges,
347 : TimeDependentMaps>;
348 :
349 : template <typename Metavariables>
350 0 : using options = tmpl::conditional_t<
351 : domain::BoundaryConditions::has_boundary_conditions_base_v<
352 : typename Metavariables::system>,
353 : tmpl::push_back<
354 : basic_options,
355 : OuterBoundaryCondition<
356 : domain::BoundaryConditions::get_boundary_conditions_base<
357 : typename Metavariables::system>>>,
358 : basic_options>;
359 :
360 0 : static constexpr Options::String help{
361 : "A 3D cubed sphere. Six wedges surround an interior region, which is "
362 : "either excised or filled in with a seventh block. The interior region "
363 : "is a (possibly deformed) sphere when excised, or a (possibly deformed) "
364 : "cube when filled in. Additional spherical shells, each composed of six "
365 : "wedges, can be added with the 'RadialPartitioning' option."};
366 :
367 0 : Sphere(
368 : double inner_radius, double outer_radius,
369 : std::variant<Excision, InnerCube> interior,
370 : const typename InitialRefinement::type& initial_refinement,
371 : const typename InitialGridPoints::type& initial_number_of_grid_points,
372 : bool use_equiangular_map,
373 : std::optional<EquatorialCompressionOptions> equatorial_compression = {},
374 : std::vector<double> radial_partitioning = {},
375 : const typename RadialDistribution::type& radial_distribution =
376 : domain::CoordinateMaps::Distribution::Linear,
377 : ShellWedges which_wedges = ShellWedges::All,
378 : std::optional<TimeDepOptionType> time_dependent_options = std::nullopt,
379 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
380 : outer_boundary_condition = nullptr,
381 : const Options::Context& context = {});
382 :
383 0 : Sphere() = default;
384 0 : Sphere(const Sphere&) = delete;
385 0 : Sphere(Sphere&&) = default;
386 0 : Sphere& operator=(const Sphere&) = delete;
387 0 : Sphere& operator=(Sphere&&) = default;
388 0 : ~Sphere() override = default;
389 :
390 0 : Domain<3> create_domain() const override;
391 :
392 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
393 1 : grid_anchors() const override {
394 : return grid_anchors_;
395 : }
396 :
397 : std::vector<DirectionMap<
398 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
399 1 : external_boundary_conditions() const override;
400 :
401 1 : std::vector<std::array<size_t, 3>> initial_extents() const override {
402 : return initial_number_of_grid_points_;
403 : }
404 :
405 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels()
406 : const override {
407 : return initial_refinement_;
408 : }
409 :
410 1 : std::vector<std::string> block_names() const override { return block_names_; }
411 :
412 : std::unordered_map<std::string, std::unordered_set<std::string>>
413 1 : block_groups() const override {
414 : return block_groups_;
415 : }
416 :
417 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
418 : initial_expiration_times = {}) const
419 : -> std::unordered_map<
420 : std::string,
421 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
422 :
423 : private:
424 0 : double inner_radius_{};
425 0 : double outer_radius_{};
426 0 : std::variant<Excision, InnerCube> interior_{};
427 0 : bool fill_interior_ = false;
428 0 : std::vector<std::array<size_t, 3>> initial_refinement_{};
429 0 : std::vector<std::array<size_t, 3>> initial_number_of_grid_points_{};
430 0 : bool use_equiangular_map_ = false;
431 0 : std::optional<EquatorialCompressionOptions> equatorial_compression_{};
432 0 : std::vector<double> radial_partitioning_{};
433 0 : std::vector<domain::CoordinateMaps::Distribution> radial_distribution_{};
434 0 : ShellWedges which_wedges_ = ShellWedges::All;
435 0 : std::optional<TimeDepOptionType> time_dependent_options_{};
436 0 : bool use_hard_coded_maps_{false};
437 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
438 0 : outer_boundary_condition_;
439 0 : size_t num_shells_{};
440 0 : size_t num_blocks_;
441 0 : size_t num_blocks_per_shell_;
442 0 : std::vector<std::string> block_names_{};
443 : std::unordered_map<std::string, std::unordered_set<std::string>>
444 0 : block_groups_{};
445 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
446 0 : grid_anchors_{};
447 : };
448 :
449 : } // namespace domain::creators
|
