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 <unordered_set>
13 : #include <variant>
14 : #include <vector>
15 :
16 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
17 : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
18 : #include "Domain/CoordinateMaps/BulgedCube.hpp"
19 : #include "Domain/CoordinateMaps/Equiangular.hpp"
20 : #include "Domain/CoordinateMaps/ProductMaps.hpp"
21 : #include "Domain/Creators/DomainCreator.hpp"
22 : #include "Domain/Creators/Sphere.hpp"
23 : #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
24 : #include "Domain/Creators/TimeDependentOptions/Sphere.hpp"
25 : #include "Options/Auto.hpp"
26 : #include "Options/Context.hpp"
27 : #include "Options/String.hpp"
28 : #include "Utilities/TMPL.hpp"
29 :
30 : /// \cond
31 : template <size_t Dim, typename T>
32 : class DirectionMap;
33 : template <size_t Dim>
34 : class Domain;
35 : namespace domain {
36 : namespace CoordinateMaps {
37 : class Affine;
38 : template <size_t Dim>
39 : class Identity;
40 : class Interval;
41 : template <typename Map1, typename Map2>
42 : class ProductOf2Maps;
43 : class SphericalToCartesianPfaffian;
44 : template <size_t Dim>
45 : class Wedge;
46 : } // namespace CoordinateMaps
47 :
48 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
49 : class CoordinateMap;
50 : } // namespace domain
51 : /// \endcond
52 :
53 : namespace domain::creators {
54 : /*!
55 : * \brief A set of non-conforming concentric spherical shells
56 : *
57 : * \details The inner spherical shells are decomposed into six wedges
58 : * surrounding an optionally excised interior region. The outer spherical
59 : * shells will use a spherical harmonic basis which cannot be used with subcell.
60 : *
61 : * This domain creator offers one grid anchor "Center" at the origin.
62 : *
63 : */
64 1 : class NonconformingSphericalShells : public DomainCreator<3> {
65 : private:
66 0 : using Affine = CoordinateMaps::Affine;
67 0 : using Affine3D = CoordinateMaps::ProductOf3Maps<Affine, Affine, Affine>;
68 0 : using Equiangular = CoordinateMaps::Equiangular;
69 0 : using Equiangular3D =
70 : CoordinateMaps::ProductOf3Maps<Equiangular, Equiangular, Equiangular>;
71 0 : using BulgedCube = CoordinateMaps::BulgedCube;
72 :
73 : public:
74 0 : using maps_list = tmpl::list<
75 : // Inner cube
76 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, BulgedCube>,
77 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial, Affine3D>,
78 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
79 : Equiangular3D>,
80 : // Wedges
81 : domain::CoordinateMap<Frame::BlockLogical, Frame::Inertial,
82 : CoordinateMaps::Wedge<3>>,
83 : // Spherical shells
84 : domain::CoordinateMap<
85 : Frame::BlockLogical, Frame::Inertial,
86 : domain::CoordinateMaps::ProductOf2Maps<
87 : domain::CoordinateMaps::Interval,
88 : domain::CoordinateMaps::Identity<2>>,
89 : domain::CoordinateMaps::SphericalToCartesianPfaffian>,
90 : typename sphere::TimeDependentMapOptions::maps_list>;
91 :
92 0 : struct InnerRadius {
93 0 : using type = double;
94 0 : static constexpr Options::String help = {
95 : "Inner radius of the inner wedges."};
96 : };
97 :
98 0 : struct InterfaceRadius {
99 0 : using type = double;
100 0 : static constexpr Options::String help = {
101 : "Radius of interface between the inner wedges and the outer spherical "
102 : "shells."};
103 : };
104 :
105 0 : struct OuterRadius {
106 0 : using type = double;
107 0 : static constexpr Options::String help = {
108 : "Outer radius of the outer spherical shell."};
109 : };
110 :
111 0 : using Excision = detail::Excision;
112 0 : using InnerCube = detail::InnerCube;
113 :
114 0 : struct Interior {
115 0 : using type = std::variant<Excision, InnerCube>;
116 0 : static constexpr Options::String help = {
117 : "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
118 : "excise the interior of the sphere, leaving a spherical shell "
119 : "(or just 'Excise' if boundary conditions are disabled). "
120 : "Or specify 'FillWithSphericity' to fill the interior."};
121 : };
122 :
123 0 : struct InitialCubeRefinement {
124 0 : using type =
125 : std::variant<std::array<size_t, 2>, std::vector<std::array<size_t, 2>>,
126 : std::unordered_map<std::string, std::array<size_t, 2>>>;
127 0 : static constexpr Options::String help = {
128 : "Initial cube refinement level. Specify one of: a "
129 : "list representing [angular, r], or such a list for every block "
130 : "in the domain. The central cube always uses the angular value for all "
131 : "directions."};
132 : };
133 :
134 0 : struct InitialSHRefinement {
135 0 : using type = std::variant<size_t, std::vector<size_t>,
136 : std::unordered_map<std::string, size_t>>;
137 0 : static constexpr Options::String help = {
138 : "Initial spherical harmonic shell radial refinement level. Specify one "
139 : "of: a single number, or such a number for every block in the domain."};
140 : };
141 :
142 0 : struct InitialCubeGridPoints {
143 0 : using type =
144 : std::variant<std::array<size_t, 2>, std::vector<std::array<size_t, 2>>,
145 : std::unordered_map<std::string, std::array<size_t, 2>>>;
146 0 : static constexpr Options::String help = {
147 : "Initial number of grid points for the cube region. Specify one of: a "
148 : "list representing [angular, r], or such a list for every block "
149 : "in the domain. The central cube always uses the angular value for all "
150 : "directions."};
151 : };
152 :
153 0 : struct InitialSHGridPoints {
154 0 : using type =
155 : std::variant<std::array<size_t, 2>, std::vector<std::array<size_t, 2>>,
156 : std::unordered_map<std::string, std::array<size_t, 2>>>;
157 0 : static constexpr Options::String help = {
158 : "Initial number of grid points for the spherical harmonic shells. "
159 : "Specify one of: a list representing [l_max, r], or such a list for "
160 : "every block in the domain."};
161 : };
162 :
163 0 : struct RadialPartitioning {
164 0 : using type = std::array<std::vector<double>, 2>;
165 0 : static constexpr Options::String help = {
166 : "Radial coordinates of the boundaries splitting the spherical shell "
167 : "between InnerRadius and InterfaceRadius and then the InterfaceRadius "
168 : "and OuterRadius. They must be given in ascending order."};
169 : };
170 :
171 0 : struct RadialDistribution {
172 0 : using type =
173 : std::array<std::vector<domain::CoordinateMaps::Distribution>, 2>;
174 0 : static constexpr Options::String help = {
175 : "Select the radial distribution of grid points in each spherical "
176 : "shell. There must be N+1 radial distributions specified for N radial "
177 : "partitions for both the wedges and spherical shells. If the interior "
178 : "of the sphere is filled with a cube, the innermost shell must have a "
179 : "'Linear' distribution because it changes in sphericity."};
180 : };
181 :
182 0 : struct UseEquiangularMap {
183 0 : using type = bool;
184 0 : static constexpr Options::String help = {
185 : "Use equiangular instead of equidistant coordinates. Equiangular "
186 : "coordinates give better gridpoint spacings in the angular "
187 : "directions, while equidistant coordinates give better gridpoint "
188 : "spacings in the inner cube."};
189 : };
190 :
191 0 : using TimeDepOptionType = std::variant<
192 : sphere::TimeDependentMapOptions,
193 : std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>>;
194 :
195 0 : struct TimeDependentMaps {
196 0 : using type = Options::Auto<TimeDepOptionType, Options::AutoLabel::None>;
197 0 : static constexpr Options::String help = {
198 : "The options for time dependent maps. This can either be a "
199 : "TimeDependence or hard coded time dependent options. Specify `None` "
200 : "for no time dependent maps."};
201 : };
202 :
203 : template <typename BoundaryConditionsBase>
204 0 : struct OuterBoundaryCondition {
205 0 : static constexpr Options::String help =
206 : "Options for the boundary conditions at the outer radius.";
207 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
208 : };
209 :
210 0 : using basic_options =
211 : tmpl::list<InnerRadius, InterfaceRadius, OuterRadius, Interior,
212 : InitialCubeRefinement, InitialSHRefinement,
213 : InitialCubeGridPoints, InitialSHGridPoints, RadialPartitioning,
214 : RadialDistribution, UseEquiangularMap, TimeDependentMaps>;
215 :
216 : template <typename Metavariables>
217 0 : using options = tmpl::conditional_t<
218 : domain::BoundaryConditions::has_boundary_conditions_base_v<
219 : typename Metavariables::system>,
220 : tmpl::push_back<
221 : basic_options,
222 : OuterBoundaryCondition<
223 : domain::BoundaryConditions::get_boundary_conditions_base<
224 : typename Metavariables::system>>>,
225 : basic_options>;
226 :
227 0 : static constexpr Options::String help{
228 : "A set of concentric spherical shells centered at the origin."};
229 :
230 0 : NonconformingSphericalShells(
231 : double inner_radius, double interface_radius, double outer_radius,
232 : std::variant<Excision, InnerCube> interior,
233 : const typename InitialCubeRefinement::type& initial_cube_refinement,
234 : const typename InitialSHRefinement::type& initial_sh_refinement,
235 : const typename InitialCubeGridPoints::type& initial_cube_grid_points,
236 : const typename InitialSHGridPoints::type& initial_sh_grid_points,
237 : std::array<std::vector<double>, 2> radial_partitioning = {},
238 : std::array<std::vector<domain::CoordinateMaps::Distribution>, 2>
239 : radial_distribution =
240 : {std::vector<domain::CoordinateMaps::Distribution>{
241 : domain::CoordinateMaps::Distribution::Linear},
242 : std::vector<domain::CoordinateMaps::Distribution>{
243 : domain::CoordinateMaps::Distribution::Linear}},
244 : bool use_equiangular_map = true,
245 : std::optional<TimeDepOptionType> time_dependent_options = std::nullopt,
246 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
247 : outer_boundary_condition = nullptr,
248 : const Options::Context& context = {});
249 :
250 0 : NonconformingSphericalShells() = default;
251 0 : NonconformingSphericalShells(const NonconformingSphericalShells&) = delete;
252 0 : NonconformingSphericalShells(NonconformingSphericalShells&&) = default;
253 0 : NonconformingSphericalShells& operator=(const NonconformingSphericalShells&) =
254 : delete;
255 0 : NonconformingSphericalShells& operator=(NonconformingSphericalShells&&) =
256 : default;
257 0 : ~NonconformingSphericalShells() override = default;
258 :
259 0 : Domain<3> create_domain() const override;
260 :
261 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
262 1 : grid_anchors() const override;
263 :
264 : std::vector<DirectionMap<
265 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
266 1 : external_boundary_conditions() const override;
267 :
268 1 : std::vector<std::string> block_names() const override { return block_names_; }
269 :
270 : std::unordered_map<std::string, std::unordered_set<std::string>>
271 1 : block_groups() const override {
272 : return block_groups_;
273 : }
274 :
275 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
276 :
277 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
278 :
279 1 : auto functions_of_time(const std::unordered_map<std::string, double>&
280 : initial_expiration_times = {}) const
281 : -> std::unordered_map<
282 : std::string,
283 : std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
284 :
285 : private:
286 0 : double inner_radius_{};
287 0 : double interface_radius_{};
288 0 : double outer_radius_{};
289 0 : std::variant<Excision, InnerCube> interior_{};
290 0 : bool fill_interior_ = false;
291 0 : std::vector<std::array<size_t, 2>> initial_cube_refinement_{};
292 0 : std::vector<size_t> initial_sh_refinement_{};
293 0 : std::vector<std::array<size_t, 2>> initial_cube_grid_points_{};
294 0 : std::vector<std::array<size_t, 2>> initial_sh_grid_points_{};
295 0 : std::array<std::vector<double>, 2> radial_partitioning_;
296 : std::array<std::vector<domain::CoordinateMaps::Distribution>, 2>
297 0 : radial_distribution_;
298 0 : bool use_equiangular_map_ = false;
299 0 : std::optional<TimeDepOptionType> time_dependent_options_{};
300 0 : bool use_hard_coded_maps_{false};
301 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
302 0 : outer_boundary_condition_{};
303 0 : std::vector<std::string> block_names_;
304 : std::unordered_map<std::string, std::unordered_set<std::string>>
305 0 : block_groups_;
306 0 : size_t num_blocks_{};
307 0 : size_t num_cube_shells_{};
308 0 : size_t num_sh_shells_{};
309 : std::unordered_map<std::string, tnsr::I<double, 3, Frame::Grid>>
310 0 : grid_anchors_{};
311 : };
312 : } // namespace domain::creators
|
