Merge pull request #4718 from pajkosmi/inner_cube_surrounded_by_spher…

…es_domain

Add domain: inner cube surrounded by sphere with radial partitions

src/Domain/Creators/Python/Sphere.cpp

@@ -17,7 +17,7 @@ namespace py = pybind11;
 namespace domain::creators::py_bindings {
 void bind_sphere(py::module& m) {
   py::class_<Sphere, DomainCreator<3>>(m, "Sphere")
-      .def(py::init<double, double, double, size_t, std::array<size_t, 2>,
+      .def(py::init<double, double, double, size_t, std::array<size_t, 3>,
                     bool>(),
            py::arg("inner_radius"), py::arg("outer_radius"),
            py::arg("inner_cube_sphericity"), py::arg("initial_refinement"),

src/Domain/Creators/Sphere.cpp

@@ -7,7 +7,9 @@
 #include <cmath>
 #include <cstddef>
 #include <memory>
+#include <unordered_map>
 #include <utility>
+#include <variant>
 #include <vector>
 
 #include "Domain/Block.hpp"
@@ -22,6 +24,7 @@
 #include "Domain/CoordinateMaps/ProductMaps.tpp"
 #include "Domain/CoordinateMaps/Wedge.hpp"
 #include "Domain/Creators/DomainCreator.hpp"
+#include "Domain/Creators/ExpandOverBlocks.hpp"
 #include "Domain/Creators/TimeDependence/None.hpp"
 #include "Domain/Domain.hpp"
 #include "Domain/DomainHelpers.hpp"
@@ -39,33 +42,121 @@ struct BlockLogical;
 
 namespace domain::creators {
 Sphere::Sphere(
-    typename InnerRadius::type inner_radius,
-    typename OuterRadius::type outer_radius, const double inner_cube_sphericity,
-    typename InitialRefinement::type initial_refinement,
-    typename InitialGridPoints::type initial_number_of_grid_points,
-    typename UseEquiangularMap::type use_equiangular_map,
+    double inner_radius, double outer_radius, double inner_cube_sphericity,
+    const typename InitialRefinement::type& initial_refinement,
+    const typename InitialGridPoints::type& initial_number_of_grid_points,
+    bool use_equiangular_map, std::vector<double> radial_partitioning,
+    std::vector<domain::CoordinateMaps::Distribution> radial_distribution,
     std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
         time_dependence,
     std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
         boundary_condition,
     const Options::Context& context)
-    // clang-tidy: trivially copyable
-    : inner_radius_(std::move(inner_radius)),  // NOLINT
-      outer_radius_(std::move(outer_radius)),  // NOLINT
+    : inner_radius_(inner_radius),
+      outer_radius_(outer_radius),
       inner_cube_sphericity_(inner_cube_sphericity),
-      initial_refinement_(                                   // NOLINT
-          std::move(initial_refinement)),                    // NOLINT
-      initial_number_of_grid_points_(                        // NOLINT
-          std::move(initial_number_of_grid_points)),         // NOLINT
-      use_equiangular_map_(std::move(use_equiangular_map)),  // NOLINT
-      time_dependence_(std::move(time_dependence)),          // NOLINT
+      use_equiangular_map_(use_equiangular_map),
+      radial_partitioning_(std::move(radial_partitioning)),
+      radial_distribution_(std::move(radial_distribution)),
+      time_dependence_(std::move(time_dependence)),
+
       boundary_condition_(std::move(boundary_condition)) {
   if (inner_cube_sphericity_ < 0.0 or inner_cube_sphericity_ >= 1.0) {
     PARSE_ERROR(
         context,
         "Inner cube sphericity must be >= 0.0 and strictly < 1.0, not " +
             get_output(inner_cube_sphericity_));
   }
+
+  if (inner_radius_ > outer_radius_) {
+    PARSE_ERROR(context,
+                "Inner radius must be smaller than outer radius, but inner "
+                "radius is " +
+                    std::to_string(inner_radius_) + " and outer radius is " +
+                    std::to_string(outer_radius_) + ".");
+  }
+  if (not std::is_sorted(radial_partitioning_.begin(),
+                         radial_partitioning_.end())) {
+    PARSE_ERROR(context,
+                "Specify radial partitioning in ascending order. Specified "
+                "radial partitioning is: " +
+                    get_output(radial_partitioning_));
+  }
+  if (not radial_partitioning_.empty()) {
+    if (radial_partitioning_.front() <= inner_radius_) {
+      PARSE_ERROR(context,
+                  "First radial partition must be larger than inner "
+                  "radius, but is: " +
+                      std::to_string(inner_radius_));
+    }
+    if (radial_partitioning_.back() >= outer_radius_) {
+      PARSE_ERROR(context,
+                  "Last radial partition must be smaller than outer "
+                  "radius, but is: " +
+                      std::to_string(outer_radius_));
+    }
+  }
+
+  const size_t num_shells = 1 + radial_partitioning_.size();
+  if (radial_distribution_.size() != num_shells) {
+    PARSE_ERROR(context,
+                "Specify a 'RadialDistribution' for every spherical shell. You "
+                "specified "
+                    << radial_distribution_.size()
+                    << " items, but the domain has " << num_shells
+                    << " shells.");
+  }
+  if (radial_distribution_.front() !=
+      domain::CoordinateMaps::Distribution::Linear) {
+    PARSE_ERROR(context,
+                "The 'RadialDistribution' must be 'Linear' for the innermost "
+                "shell because it changes in sphericity. Add entries to "
+                "'RadialPartitioning' to add outer shells for which you can "
+                "select different radial distributions.");
+  }
+
+  // Create block names and groups
+  static std::array<std::string, 6> wedge_directions{
+      "UpperZ", "LowerZ", "UpperY", "LowerY", "UpperX", "LowerX"};
+  for (size_t shell = 0; shell < num_shells; ++shell) {
+    std::string shell_prefix = "Shell" + std::to_string(shell);
+    for (size_t direction = 0; direction < 6; ++direction) {
+      const std::string wedge_name =
+          shell_prefix + gsl::at(wedge_directions, direction);
+      block_names_.emplace_back(wedge_name);
+      if (num_shells > 1) {
+        block_groups_[shell_prefix].insert(wedge_name);
+      }
+      block_groups_["Wedges"].insert(wedge_name);
+    }
+  }
+  block_names_.emplace_back("InnerCube");
+
+  // Expand initial refinement and number of grid points over all blocks
+  const ExpandOverBlocks<size_t, 3> expand_over_blocks{block_names_,
+                                                       block_groups_};
+  try {
+    initial_refinement_ = std::visit(expand_over_blocks, initial_refinement);
+  } catch (const std::exception& error) {
+    PARSE_ERROR(context, "Invalid 'InitialRefinement': " << error.what());
+  }
+  try {
+    initial_number_of_grid_points_ =
+        std::visit(expand_over_blocks, initial_number_of_grid_points);
+  } catch (const std::exception& error) {
+    PARSE_ERROR(context, "Invalid 'InitialGridPoints': " << error.what());
+  }
+
+  // The central cube has no notion of a "radial" direction, so we set
+  // refinement and number of grid points of the central cube z direction to
+  // its y value, which corresponds to the azimuthal direction of the
+  // wedges. This keeps the boundaries conforming when the radial direction
+  // is chosen differently to the angular directions.
+  auto& central_cube_refinement = initial_refinement_.back();
+  auto& central_cube_grid_points = initial_number_of_grid_points_.back();
+  central_cube_refinement[2] = central_cube_refinement[1];
+  central_cube_grid_points[2] = central_cube_grid_points[1];
+
   using domain::BoundaryConditions::is_none;
   if (is_none(boundary_condition_)) {
     PARSE_ERROR(
@@ -85,14 +176,14 @@ Sphere::Sphere(
 }
 
 Domain<3> Sphere::create_domain() const {
+  const size_t num_shells = 1 + radial_partitioning_.size();
   std::vector<std::array<size_t, 8>> corners =
-      corners_for_radially_layered_domains(1, true);
+      corners_for_radially_layered_domains(num_shells, true);
 
-  auto coord_maps = domain::make_vector_coordinate_map_base<Frame::BlockLogical,
-                                                            Frame::Inertial, 3>(
-      sph_wedge_coordinate_maps(inner_radius_, outer_radius_,
-                                inner_cube_sphericity_, 1.0,
-                                use_equiangular_map_));
+  auto coord_maps = domain::make_vector_coordinate_map_base<
+      Frame::BlockLogical, Frame::Inertial, 3>(sph_wedge_coordinate_maps(
+      inner_radius_, outer_radius_, inner_cube_sphericity_, 1.0,
+      use_equiangular_map_, false, radial_partitioning_, radial_distribution_));
   if (inner_cube_sphericity_ == 0.0) {
     if (use_equiangular_map_) {
       coord_maps.emplace_back(
@@ -121,7 +212,8 @@ Domain<3> Sphere::create_domain() const {
                        use_equiangular_map_}));
   }
 
-  Domain<3> domain{std::move(coord_maps), corners};
+  Domain<3> domain{std::move(coord_maps), corners,      {}, {},
+                   block_names_,          block_groups_};
 
   if (not time_dependence_->is_none()) {
     const size_t number_of_blocks = domain.blocks().size();
@@ -148,31 +240,19 @@ Sphere::external_boundary_conditions() const {
   if (boundary_condition_ == nullptr) {
     return {};
   }
+
+  // number of blocks = 1 inner_block + 6 * (number of shells)
+  size_t number_of_blocks = 1 + 6 * (radial_partitioning_.size() + 1);
+
   std::vector<DirectionMap<
       3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
-      boundary_conditions{7};
-  for (size_t i = 0; i < 6; ++i) {
+      boundary_conditions{number_of_blocks};
+  for (size_t i = number_of_blocks - 7; i < number_of_blocks - 1; ++i) {
     boundary_conditions[i][Direction<3>::upper_zeta()] =
         boundary_condition_->get_clone();
   }
   return boundary_conditions;
 }
-
-std::vector<std::array<size_t, 3>> Sphere::initial_extents() const {
-  std::vector<std::array<size_t, 3>> extents{
-      6,
-      {{initial_number_of_grid_points_[1], initial_number_of_grid_points_[1],
-        initial_number_of_grid_points_[0]}}};
-  extents.push_back(
-      {{initial_number_of_grid_points_[1], initial_number_of_grid_points_[1],
-        initial_number_of_grid_points_[1]}});
-  return extents;
-}
-
-std::vector<std::array<size_t, 3>> Sphere::initial_refinement_levels() const {
-  return {7, make_array<3>(initial_refinement_)};
-}
-
 std::unordered_map<std::string,
                    std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>>
 Sphere::functions_of_time(const std::unordered_map<std::string, double>&

src/Domain/Creators/Sphere.hpp

@@ -6,10 +6,14 @@
 #include <array>
 #include <cstddef>
 #include <memory>
+#include <string>
+#include <unordered_map>
+#include <variant>
 #include <vector>
 
 #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
 #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
+#include "Domain/CoordinateMaps/Distribution.hpp"
 #include "Domain/Creators/DomainCreator.hpp"
 #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
 #include "Domain/Domain.hpp"
@@ -83,21 +87,53 @@ class Sphere : public DomainCreator<3> {
   };
 
   struct InitialRefinement {
-    using type = size_t;
+    using type =
+        std::variant<size_t, std::array<size_t, 3>,
+                     std::vector<std::array<size_t, 3>>,
+                     std::unordered_map<std::string, std::array<size_t, 3>>>;
     static constexpr Options::String help = {
-        "Initial refinement level in each dimension."};
+        "Initial refinement level. Specify one of: a single number, a "
+        "list representing [phi, theta, r], or such a list for every block "
+        "in the domain. The central cube always uses the value for 'theta' "
+        "in both y- and z-direction."};
   };
 
   struct InitialGridPoints {
-    using type = std::array<size_t, 2>;
+    using type =
+        std::variant<size_t, std::array<size_t, 3>,
+                     std::vector<std::array<size_t, 3>>,
+                     std::unordered_map<std::string, std::array<size_t, 3>>>;
     static constexpr Options::String help = {
-        "Initial number of grid points in [r,angular]."};
+        "Initial number of grid points. Specify one of: a single number, a "
+        "list representing [phi, theta, r], or such a list for every block "
+        "in the domain. The central cube always uses the value for 'theta' "
+        "in both y- and z-direction."};
   };
 
   struct UseEquiangularMap {
     using type = bool;
     static constexpr Options::String help = {
-        "Use equiangular instead of equidistant coordinates."};
+        "Use equiangular instead of equidistant coordinates. Equiangular "
+        "coordinates give better gridpoint spacings in the angular "
+        "directions, while equidistant coordinates give better gridpoint "
+        "spacings in the center block."};
+  };
+
+  struct RadialPartitioning {
+    using type = std::vector<double>;
+    static constexpr Options::String help = {
+        "Radial coordinates of the boundaries splitting the spherical shell "
+        "between InnerRadius and OuterRadius."};
+  };
+
+  struct RadialDistribution {
+    using type = std::vector<domain::CoordinateMaps::Distribution>;
+    static constexpr Options::String help = {
+        "Select the radial distribution of grid points in each spherical "
+        "shell. The innermost shell must have a 'Linear' distribution because "
+        "it changes in sphericity. The 'RadialPartitioning' determines the "
+        "number of shells."};
+    static size_t lower_bound_on_size() { return 1; }
   };
 
   struct TimeDependence {
@@ -118,7 +154,7 @@ class Sphere : public DomainCreator<3> {
   using basic_options =
       tmpl::list<InnerRadius, OuterRadius, InnerCubeSphericity,
                  InitialRefinement, InitialGridPoints, UseEquiangularMap,
-                 TimeDependence>;
+                 RadialPartitioning, RadialDistribution, TimeDependence>;
 
   template <typename Metavariables>
   using options = tmpl::conditional_t<
@@ -132,22 +168,16 @@ class Sphere : public DomainCreator<3> {
       basic_options>;
 
   static constexpr Options::String help{
-      "Creates a 3D Sphere with seven Blocks.\n"
-      "Only one refinement level for all dimensions is currently supported.\n"
-      "The number of gridpoints in the radial direction can be set\n"
-      "independently of the number of gridpoints in the angular directions.\n"
-      "The number of gridpoints along the dimensions of the cube is equal\n"
-      "to the number of gridpoints along the angular dimensions of the "
-      "wedges.\n"
-      "Equiangular coordinates give better gridpoint spacings in the angular\n"
-      "directions, while equidistant coordinates give better gridpoint\n"
-      "spacings in the center block."};
-
-  Sphere(typename InnerRadius::type inner_radius,
-         typename OuterRadius::type outer_radius, double inner_cube_sphericity,
-         typename InitialRefinement::type initial_refinement,
-         typename InitialGridPoints::type initial_number_of_grid_points,
-         typename UseEquiangularMap::type use_equiangular_map,
+      "A 3D cubed sphere. Six wedges surround a central cube. Additional "
+      "spherical shells, each composed of six wedges, can be added with the "
+      "'RadialPartitioning' option."};
+
+  Sphere(double inner_radius, double outer_radius, double inner_cube_sphericity,
+         const typename InitialRefinement::type& initial_refinement,
+         const typename InitialGridPoints::type& initial_number_of_grid_points,
+         bool use_equiangular_map, std::vector<double> radial_partitioning = {},
+         std::vector<domain::CoordinateMaps::Distribution> radial_distribution =
+             {domain::CoordinateMaps::Distribution::Linear},
          std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
              time_dependence = nullptr,
          std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
@@ -167,9 +197,21 @@ class Sphere : public DomainCreator<3> {
       3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
   external_boundary_conditions() const override;
 
-  std::vector<std::array<size_t, 3>> initial_extents() const override;
+  std::vector<std::array<size_t, 3>> initial_extents() const override {
+    return initial_number_of_grid_points_;
+  }
+
+  std::vector<std::array<size_t, 3>> initial_refinement_levels()
+      const override {
+    return initial_refinement_;
+  }
+
+  std::vector<std::string> block_names() const override { return block_names_; }
 
-  std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
+  std::unordered_map<std::string, std::unordered_set<std::string>>
+  block_groups() const override {
+    return block_groups_;
+  }
 
   auto functions_of_time(const std::unordered_map<std::string, double>&
                              initial_expiration_times = {}) const
@@ -178,16 +220,21 @@ class Sphere : public DomainCreator<3> {
           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
 
  private:
-  typename InnerRadius::type inner_radius_{};
-  typename OuterRadius::type outer_radius_{};
+  double inner_radius_{};
+  double outer_radius_{};
   double inner_cube_sphericity_{};
-  typename InitialRefinement::type initial_refinement_{};
-  typename InitialGridPoints::type initial_number_of_grid_points_{};
-  typename UseEquiangularMap::type use_equiangular_map_ = false;
+  std::vector<std::array<size_t, 3>> initial_refinement_{};
+  std::vector<std::array<size_t, 3>> initial_number_of_grid_points_{};
+  bool use_equiangular_map_ = false;
+  std::vector<double> radial_partitioning_{};
+  std::vector<domain::CoordinateMaps::Distribution> radial_distribution_{};
   std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
       time_dependence_;
   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
       boundary_condition_;
+  std::vector<std::string> block_names_{};
+  std::unordered_map<std::string, std::unordered_set<std::string>>
+      block_groups_{};
 };
 }  // namespace creators
 }  // namespace domain