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