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 <unordered_map>
10 : #include <unordered_set>
11 : #include <variant>
12 : #include <vector>
13 :
14 : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
15 : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
16 : #include "Domain/CoordinateMaps/Distribution.hpp"
17 : #include "Domain/Creators/DomainCreator.hpp" // IWYU pragma: keep
18 : #include "Domain/Domain.hpp"
19 : #include "Domain/Structure/DirectionMap.hpp"
20 : #include "Options/Context.hpp"
21 : #include "Options/String.hpp"
22 : #include "Utilities/TMPL.hpp"
23 :
24 : /// \cond
25 : namespace domain {
26 : namespace CoordinateMaps {
27 : class Interval;
28 : template <typename Map1, typename Map2>
29 : class ProductOf2Maps;
30 : template <typename Map1, typename Map2, typename Map3>
31 : class ProductOf3Maps;
32 : template <size_t Dim>
33 : class Wedge;
34 : } // namespace CoordinateMaps
35 :
36 : template <typename SourceFrame, typename TargetFrame, typename... Maps>
37 : class CoordinateMap;
38 : } // namespace domain
39 : /// \endcond
40 :
41 : namespace domain::creators {
42 : /// Create a 3D Domain in the shape of a cylinder where the cross-section
43 : /// is a square surrounded by four two-dimensional wedges (see `Wedge`).
44 : ///
45 : /// The outer shell can be split into sub-shells and the cylinder can be split
46 : /// into disks along its height.
47 : /// The block numbering starts at the inner square and goes counter-clockwise,
48 : /// starting with the eastern wedge (+x-direction), through consecutive shells,
49 : /// then repeats this pattern for all layers bottom to top.
50 : ///
51 : /// \image html Cylinder.png "The Cylinder Domain."
52 1 : class Cylinder : public DomainCreator<3> {
53 : public:
54 0 : using maps_list =
55 : tmpl::list<domain::CoordinateMap<
56 : Frame::BlockLogical, Frame::Inertial,
57 : CoordinateMaps::ProductOf3Maps<CoordinateMaps::Interval,
58 : CoordinateMaps::Interval,
59 : CoordinateMaps::Interval>>,
60 : domain::CoordinateMap<
61 : Frame::BlockLogical, Frame::Inertial,
62 : CoordinateMaps::ProductOf2Maps<CoordinateMaps::Wedge<2>,
63 : CoordinateMaps::Interval>>>;
64 :
65 0 : struct InnerRadius {
66 0 : using type = double;
67 0 : static constexpr Options::String help = {
68 : "Radius of the circle circumscribing the inner square."};
69 0 : static double lower_bound() { return 0.; }
70 : };
71 :
72 0 : struct OuterRadius {
73 0 : using type = double;
74 0 : static constexpr Options::String help = {"Radius of the cylinder."};
75 0 : static double lower_bound() { return 0.; }
76 : };
77 :
78 0 : struct LowerZBound {
79 0 : using type = double;
80 0 : static constexpr Options::String help = {
81 : "z-coordinate of the base of the cylinder."};
82 : };
83 :
84 0 : struct UpperZBound {
85 0 : using type = double;
86 0 : static constexpr Options::String help = {
87 : "z-coordinate of the top of the cylinder."};
88 : };
89 :
90 0 : struct IsPeriodicInZ {
91 0 : using type = bool;
92 0 : static constexpr Options::String help = {
93 : "True if periodic in the cylindrical z direction."};
94 : };
95 :
96 0 : struct InitialRefinement {
97 0 : using type =
98 : std::variant<size_t, std::array<size_t, 3>,
99 : std::vector<std::array<size_t, 3>>,
100 : std::unordered_map<std::string, std::array<size_t, 3>>>;
101 0 : static constexpr Options::String help = {
102 : "Initial refinement level. Specify one of: a single number, a list "
103 : "representing [r, theta, z], or such a list for every block in the "
104 : "domain. The central cube always uses the value for 'theta' in both "
105 : "x- and y-direction."};
106 : };
107 :
108 0 : struct InitialGridPoints {
109 0 : using type =
110 : std::variant<size_t, std::array<size_t, 3>,
111 : std::vector<std::array<size_t, 3>>,
112 : std::unordered_map<std::string, std::array<size_t, 3>>>;
113 0 : static constexpr Options::String help = {
114 : "Initial number of grid points. Specify one of: a single number, a "
115 : "list representing [r, theta, z], or such a list for every block in "
116 : "the domain. The central cube always uses the value for 'theta' in "
117 : "both x- and y-direction."};
118 : };
119 :
120 0 : struct UseEquiangularMap {
121 0 : using type = bool;
122 0 : static constexpr Options::String help = {
123 : "Use equiangular instead of equidistant coordinates."};
124 : };
125 :
126 0 : struct RadialPartitioning {
127 0 : using type = std::vector<double>;
128 0 : static constexpr Options::String help = {
129 : "Radial coordinates of the boundaries splitting the outer shell "
130 : "between InnerRadius and OuterRadius."};
131 : };
132 :
133 0 : struct PartitioningInZ {
134 0 : using type = std::vector<double>;
135 0 : static constexpr Options::String help = {
136 : "z-coordinates of the boundaries splitting the domain into layers "
137 : "between LowerZBound and UpperZBound."};
138 : };
139 :
140 0 : struct RadialDistribution {
141 0 : using type = std::vector<domain::CoordinateMaps::Distribution>;
142 0 : static constexpr Options::String help = {
143 : "Select the radial distribution of grid points in each cylindrical "
144 : "shell. The innermost shell must have a 'Linear' distribution because "
145 : "it changes in circularity. The 'RadialPartitioning' determines the "
146 : "number of shells."};
147 0 : static size_t lower_bound_on_size() { return 1; }
148 : };
149 :
150 0 : struct DistributionInZ {
151 0 : using type = std::vector<domain::CoordinateMaps::Distribution>;
152 0 : static constexpr Options::String help = {
153 : "Select the distribution of grid points along the z-axis in each "
154 : "layer. The lowermost layer must have a 'Linear' distribution, because "
155 : "both a 'Logarithmic' and 'Inverse' distribution places its "
156 : "singularity at 'LowerZBound'. The 'PartitioningInZ' determines the "
157 : "number of layers."};
158 0 : static size_t lower_bound_on_size() { return 1; }
159 : };
160 :
161 0 : struct BoundaryConditions {
162 0 : static constexpr Options::String help =
163 : "Options for the boundary conditions";
164 : };
165 :
166 : template <typename BoundaryConditionsBase>
167 0 : struct LowerZBoundaryCondition {
168 0 : using group = BoundaryConditions;
169 0 : static std::string name() { return "LowerZ"; }
170 0 : static constexpr Options::String help =
171 : "The boundary condition to be imposed on the lower base of the "
172 : "cylinder, i.e. at the `LowerZBound` in the z-direction.";
173 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
174 : };
175 :
176 : template <typename BoundaryConditionsBase>
177 0 : struct UpperZBoundaryCondition {
178 0 : using group = BoundaryConditions;
179 0 : static std::string name() { return "UpperZ"; }
180 0 : static constexpr Options::String help =
181 : "The boundary condition to be imposed on the upper base of the "
182 : "cylinder, i.e. at the `UpperZBound` in the z-direction.";
183 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
184 : };
185 :
186 : template <typename BoundaryConditionsBase>
187 0 : struct MantleBoundaryCondition {
188 0 : using group = BoundaryConditions;
189 0 : static std::string name() { return "Mantle"; }
190 0 : static constexpr Options::String help =
191 : "The boundary condition to be imposed on the mantle of the "
192 : "cylinder, i.e. at the `OuterRadius` in the radial direction.";
193 0 : using type = std::unique_ptr<BoundaryConditionsBase>;
194 : };
195 :
196 : template <typename Metavariables>
197 0 : using options = tmpl::append<
198 : tmpl::list<InnerRadius, OuterRadius, LowerZBound, UpperZBound>,
199 : tmpl::conditional_t<
200 : domain::BoundaryConditions::has_boundary_conditions_base_v<
201 : typename Metavariables::system>,
202 : tmpl::list<
203 : LowerZBoundaryCondition<
204 : domain::BoundaryConditions::get_boundary_conditions_base<
205 : typename Metavariables::system>>,
206 : UpperZBoundaryCondition<
207 : domain::BoundaryConditions::get_boundary_conditions_base<
208 : typename Metavariables::system>>,
209 : MantleBoundaryCondition<
210 : domain::BoundaryConditions::get_boundary_conditions_base<
211 : typename Metavariables::system>>>,
212 : tmpl::list<IsPeriodicInZ>>,
213 : tmpl::list<InitialRefinement, InitialGridPoints, UseEquiangularMap,
214 : RadialPartitioning, PartitioningInZ, RadialDistribution,
215 : DistributionInZ>>;
216 :
217 0 : static constexpr Options::String help{
218 : "Creates a right circular Cylinder with a square prism surrounded by \n"
219 : "wedges. \n"
220 : "The cylinder can be partitioned radially into multiple cylindrical \n"
221 : "shells as well as partitioned along the cylinder's height into \n"
222 : "multiple layers. Including this partitioning, the number of Blocks is \n"
223 : "given by (1 + 4*(1+n_s)) * (1+n_z), where n_s is the \n"
224 : "length of RadialPartitioning and n_z the length of \n"
225 : "HeightPartitioning. The block numbering starts at the inner square \n"
226 : "and goes counter-clockwise, starting with the eastern wedge \n"
227 : "(+x-direction) through consecutive shells, then repeats this pattern \n"
228 : "for all layers bottom to top. The wedges are named as follows: \n"
229 : " +x-direction: East \n"
230 : " +y-direction: North \n"
231 : " -x-direction: West \n"
232 : " -y-direction: South \n"
233 : "The circularity of the wedge changes from 0 to 1 within the first \n"
234 : "shell.\n"
235 : "Equiangular coordinates give better gridpoint spacings in the angular\n"
236 : "direction, while equidistant coordinates give better gridpoint\n"
237 : "spacings in the center block."};
238 :
239 0 : Cylinder(
240 : double inner_radius, double outer_radius, double lower_z_bound,
241 : double upper_z_bound, bool is_periodic_in_z,
242 : const typename InitialRefinement::type& initial_refinement,
243 : const typename InitialGridPoints::type& initial_number_of_grid_points,
244 : bool use_equiangular_map, std::vector<double> radial_partitioning = {},
245 : std::vector<double> partitioning_in_z = {},
246 : std::vector<domain::CoordinateMaps::Distribution> radial_distribution =
247 : {domain::CoordinateMaps::Distribution::Linear},
248 : std::vector<domain::CoordinateMaps::Distribution> distribution_in_z =
249 : {domain::CoordinateMaps::Distribution::Linear},
250 : const Options::Context& context = {});
251 :
252 0 : Cylinder(
253 : double inner_radius, double outer_radius, double lower_z_bound,
254 : double upper_z_bound,
255 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
256 : lower_z_boundary_condition,
257 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
258 : upper_z_boundary_condition,
259 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
260 : mantle_boundary_condition,
261 : const typename InitialRefinement::type& initial_refinement,
262 : const typename InitialGridPoints::type& initial_number_of_grid_points,
263 : bool use_equiangular_map, std::vector<double> radial_partitioning = {},
264 : std::vector<double> partitioning_in_z = {},
265 : std::vector<domain::CoordinateMaps::Distribution> radial_distribution =
266 : {domain::CoordinateMaps::Distribution::Linear},
267 : std::vector<domain::CoordinateMaps::Distribution> distribution_in_z =
268 : {domain::CoordinateMaps::Distribution::Linear},
269 : const Options::Context& context = {});
270 :
271 0 : Cylinder() = default;
272 0 : Cylinder(const Cylinder&) = delete;
273 0 : Cylinder(Cylinder&&) = default;
274 0 : Cylinder& operator=(const Cylinder&) = delete;
275 0 : Cylinder& operator=(Cylinder&&) = default;
276 0 : ~Cylinder() override = default;
277 :
278 0 : Domain<3> create_domain() const override;
279 :
280 : std::vector<DirectionMap<
281 : 3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
282 1 : external_boundary_conditions() const override;
283 :
284 1 : std::vector<std::array<size_t, 3>> initial_extents() const override;
285 :
286 1 : std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
287 :
288 1 : std::vector<std::string> block_names() const override { return block_names_; }
289 :
290 : std::unordered_map<std::string, std::unordered_set<std::string>>
291 1 : block_groups() const override {
292 : return block_groups_;
293 : }
294 :
295 : private:
296 0 : double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
297 0 : double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
298 0 : double lower_z_bound_{std::numeric_limits<double>::signaling_NaN()};
299 0 : double upper_z_bound_{std::numeric_limits<double>::signaling_NaN()};
300 0 : bool is_periodic_in_z_{true};
301 0 : std::vector<std::array<size_t, 3>> initial_refinement_{};
302 0 : std::vector<std::array<size_t, 3>> initial_number_of_grid_points_{};
303 0 : bool use_equiangular_map_{false};
304 0 : std::vector<double> radial_partitioning_{};
305 0 : std::vector<double> partitioning_in_z_{};
306 0 : std::vector<domain::CoordinateMaps::Distribution> radial_distribution_{};
307 0 : std::vector<domain::CoordinateMaps::Distribution> distribution_in_z_{};
308 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
309 0 : lower_z_boundary_condition_{};
310 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
311 0 : upper_z_boundary_condition_{};
312 : std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
313 0 : mantle_boundary_condition_{};
314 0 : std::vector<std::string> block_names_{};
315 : std::unordered_map<std::string, std::unordered_set<std::string>>
316 0 : block_groups_{};
317 : };
318 : } // namespace domain::creators
|
