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