SpECTRE Documentation Coverage Report
Current view: top level - Domain/Creators - Sphere.hpp Hit Total Coverage
Commit: 9a905b0737f373631c1b8e8389b8f26e67fa5313 Lines: 9 99 9.1 %
Date: 2024-03-28 09:03:18
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14