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 <string>
      10             : #include <unordered_map>
      11             : #include <unordered_set>
      12             : #include <variant>
      13             : #include <vector>
      14             : 
      15             : #include "Domain/BoundaryConditions/BoundaryCondition.hpp"
      16             : #include "Domain/BoundaryConditions/Cartoon.hpp"
      17             : #include "Domain/BoundaryConditions/GetBoundaryConditionsBase.hpp"
      18             : #include "Domain/CoordinateMaps/Distribution.hpp"
      19             : #include "Domain/Creators/DomainCreator.hpp"
      20             : #include "Domain/Creators/Sphere.hpp"
      21             : #include "Domain/Creators/TimeDependence/TimeDependence.hpp"
      22             : #include "Domain/Domain.hpp"
      23             : #include "Domain/Structure/DirectionMap.hpp"
      24             : #include "Options/Context.hpp"
      25             : #include "Options/String.hpp"
      26             : #include "Utilities/TMPL.hpp"
      27             : 
      28             : /// \cond
      29             : namespace domain {
      30             : namespace CoordinateMaps {
      31             : class Affine;
      32             : class Equiangular;
      33             : template <size_t Dim>
      34             : class Identity;
      35             : template <typename Map1, typename Map2>
      36             : class ProductOf2Maps;
      37             : template <typename Map1, typename Map2, typename Map3>
      38             : class ProductOf3Maps;
      39             : template <size_t Dim>
      40             : class Wedge;
      41             : template <size_t Dim>
      42             : class DiscreteRotation;
      43             : }  // namespace CoordinateMaps
      44             : 
      45             : template <typename SourceFrame, typename TargetFrame, typename... Maps>
      46             : class CoordinateMap;
      47             : }  // namespace domain
      48             : /// \endcond
      49             : 
      50             : namespace domain::creators::detail {
      51             : /// Options for filling the interior of the disk with a square
      52             : struct InnerSquare {
      53             :   static constexpr Options::String help = {
      54             :       "Fill the interior of the 2D sphere with a half-square."};
      55             :   struct Sphericity {
      56             :     static std::string name() { return "FillWithSphericity"; }
      57             :     using type = double;
      58             :     static constexpr Options::String help = {
      59             :         "Sphericity of the inner half-square. Only a sphericity of "
      60             :         "0.0 is currently implemented."};
      61             :     static double lower_bound() { return 0.0; }
      62             :     static double upper_bound() { return 0.0; }
      63             :   };
      64             :   using options = tmpl::list<Sphericity>;
      65             :   InnerSquare() = default;
      66             :   explicit InnerSquare(double sphericity_in) : sphericity(sphericity_in) {}
      67             :   double sphericity = std::numeric_limits<double>::signaling_NaN();
      68             : };
      69             : }  // namespace domain::creators::detail
      70             : 
      71             : namespace domain::creators {
      72             : /// Create a 3D Domain with a half-disk computational domain employing axial
      73             : /// symmetry. The third dimension uses a Cartoon basis with Killing vector
      74             : /// along the $\phi$ direction.
      75           1 : class CartoonSphere2D : public DomainCreator<3> {
      76             :  public:
      77           0 :   using maps_list = tmpl::list<
      78             :       domain::CoordinateMap<
      79             :           Frame::BlockLogical, Frame::Inertial,
      80             :           CoordinateMaps::ProductOf2Maps<
      81             :               CoordinateMaps::ProductOf2Maps<CoordinateMaps::Affine,
      82             :                                              CoordinateMaps::Affine>,
      83             :               CoordinateMaps::Identity<1>>>,
      84             :       domain::CoordinateMap<
      85             :           Frame::BlockLogical, Frame::Inertial,
      86             :           CoordinateMaps::ProductOf2Maps<
      87             :               CoordinateMaps::ProductOf2Maps<CoordinateMaps::Equiangular,
      88             :                                              CoordinateMaps::Equiangular>,
      89             :               CoordinateMaps::Identity<1>>>,
      90             :       domain::CoordinateMap<
      91             :           Frame::BlockLogical, Frame::Inertial,
      92             :           CoordinateMaps::DiscreteRotation<3>,
      93             :           CoordinateMaps::ProductOf2Maps<CoordinateMaps::Wedge<2>,
      94             :                                          CoordinateMaps::Identity<1>>>>;
      95             : 
      96           0 :   struct InnerRadius {
      97           0 :     using type = double;
      98           0 :     static constexpr Options::String help = {
      99             :         "Radius of the circle circumscribing the inner half-square, or the "
     100             :         "radius of the inner boundary if ExciseCenter = true."};
     101             :   };
     102             : 
     103           0 :   struct OuterRadius {
     104           0 :     using type = double;
     105           0 :     static constexpr Options::String help = {"Outer radius of the half-disk."};
     106             :   };
     107             : 
     108           0 :   struct InitialAngularRefinement {
     109           0 :     using type = size_t;
     110           0 :     static constexpr Options::String help = {
     111             :         "Initial refinement level in theta. It will be applied to all "
     112             :         "blocks because h-refinement is not supported for the ZernikeB1 "
     113             :         "basis.\n"
     114             :         "Note: the half-wedges will have their theta values decremented by "
     115             :         "one.\n"
     116             :         "Note: the inner square, if included, will have refinement set to "
     117             :         "the theta value for the center half-circle for both dimensions "
     118             :         "(decremented by one in the halved dimension)."};
     119             :   };
     120             : 
     121           0 :   struct InitialRadialRefinement {
     122           0 :     using type = std::variant<size_t, std::vector<size_t>>;
     123           0 :     static constexpr Options::String help = {
     124             :         "Initial refinement level in r. If one value is given, it will be "
     125             :         "applied to all blocks, otherwise each shell can be specificed, "
     126             :         "from innermost to outermost.\n"};
     127             :   };
     128             : 
     129           0 :   struct InitialGridPoints {
     130           0 :     using type = std::array<size_t, 2>;
     131           0 :     static constexpr Options::String help = {
     132             :         "Initial number of grid points in [r,theta]. It will be applied to all "
     133             :         "blocks because p-refinement is not supported for the ZernikeB1 "
     134             :         "basis.\n"
     135             :         "Note: if included, the inner square will have both dimensions'"
     136             :         "number of grid points set to the theta value of the surrounding "
     137             :         "wedges."};
     138             :   };
     139             : 
     140           0 :   struct RadialPartitioning {
     141           0 :     using type = std::vector<double>;
     142           0 :     static constexpr Options::String help = {
     143             :         "Radial coordinates of the boundaries splitting the radial shells. "
     144             :         "Leave empty for only one layer."};
     145             :   };
     146             : 
     147           0 :   struct UseEquiangularMap {
     148           0 :     using type = bool;
     149           0 :     static constexpr Options::String help = {
     150             :         "Use equiangular instead of equidistant coordinates in wedges."};
     151             :   };
     152             : 
     153           0 :   using Excision = detail::Excision;
     154           0 :   using InnerSquare = detail::InnerSquare;
     155             : 
     156           0 :   struct Interior {
     157           0 :     using type = std::variant<Excision, InnerSquare>;
     158           0 :     static constexpr Options::String help = {
     159             :         "Specify 'ExciseWithBoundaryCondition' and a boundary condition to "
     160             :         "excise the interior of the sphere, leaving a spherical shell "
     161             :         "(or just 'Excise' if boundary conditions are disabled). "
     162             :         "Or specify 'FillWithSphericity' to fill the interior."};
     163             :   };
     164             : 
     165             :   template <typename BoundaryConditionsBase>
     166           0 :   struct OuterBoundaryCondition {
     167           0 :     using type = std::unique_ptr<BoundaryConditionsBase>;
     168           0 :     static constexpr Options::String help = {
     169             :         "The boundary condition to impose at the outer boundary of the "
     170             :         "domain."};
     171             :   };
     172             : 
     173           0 :   struct RadialDistribution {
     174           0 :     using type = std::variant<CoordinateMaps::Distribution,
     175             :                               std::vector<CoordinateMaps::Distribution>>;
     176           0 :     static constexpr Options::String help = {
     177             :         "Distribution of grid points for the radial direction of the wedges. "
     178             :         "A single value will be applied to all shells, or every shell can be "
     179             :         "specified individually, in which case for N partitions there must be "
     180             :         "N+1 distributions. For anything but linear in the innermost shell, "
     181             :         "the center must be excised."};
     182             :   };
     183             : 
     184           0 :   struct TimeDependence {
     185           0 :     using type =
     186             :         std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>;
     187           0 :     static constexpr Options::String help = {
     188             :         "The time dependence of the moving mesh domain. Specify `None` for no "
     189             :         "time dependant maps."};
     190             :   };
     191             : 
     192           0 :   using basic_options =
     193             :       tmpl::list<InnerRadius, OuterRadius, InitialAngularRefinement,
     194             :                  InitialRadialRefinement, InitialGridPoints, RadialPartitioning,
     195             :                  UseEquiangularMap, Interior, RadialDistribution,
     196             :                  TimeDependence>;
     197             : 
     198             :   template <typename Metavariables>
     199           0 :   using options = tmpl::conditional_t<
     200             :       domain::BoundaryConditions::has_boundary_conditions_base_v<
     201             :           typename Metavariables::system>,
     202             :       tmpl::push_back<
     203             :           basic_options,
     204             :           OuterBoundaryCondition<
     205             :               domain::BoundaryConditions::get_boundary_conditions_base<
     206             :                   typename Metavariables::system>>>,
     207             :       basic_options>;
     208             : 
     209           0 :   static constexpr Options::String help{
     210             :       "Creates a sphere that requires/enforces axial-symmetry.\n"
     211             :       "The computational domain is a 2D half-disk, consisting of an inner\n"
     212             :       "half-square with two half-wedges above and below and a full wedge to\n"
     213             :       "the right. This can be extended to have mulitple shells of wedges by\n"
     214             :       "specifying a radial partition. The inner half-square can be excised,\n"
     215             :       "in which case the inner sphericity is set to 1.\n"
     216             :       "Equiangular coordinates give better gridpoint spacings in the angular\n"
     217             :       "direction, while equidistant coordinates give better gridpoint\n"
     218             :       "spacings in the center half-square.\n"
     219             :       "Elements touching the x=0 axis use ZernikeB1 bases in the x direction\n"
     220             :       "and automatically apply a system's Cartoon-type boundary condition.\n"
     221             :       "Angular refinement must be globally set as the required mortar\n"
     222             :       "projection has not been implemented."};
     223             : 
     224           0 :   CartoonSphere2D(
     225             :       double inner_radius, double outer_radius,
     226             :       size_t initial_angular_refinement,
     227             :       const typename InitialRadialRefinement::type& initial_radial_refinement,
     228             :       std::array<size_t, 2> initial_number_of_grid_points,
     229             :       std::vector<double> radial_partitioning, bool use_equiangular_map,
     230             :       std::variant<Excision, InnerSquare> interior,
     231             :       const typename RadialDistribution::type& radial_distribution,
     232             :       std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
     233             :           time_dependence,
     234             :       std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
     235             :           outer_boundary_condition,
     236             :       std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
     237             :           cartoon_boundary_condition,
     238             :       const Options::Context& context);
     239             : 
     240           0 :   CartoonSphere2D() = default;
     241           0 :   CartoonSphere2D(const CartoonSphere2D&) = delete;
     242           0 :   CartoonSphere2D(CartoonSphere2D&&) = default;
     243           0 :   CartoonSphere2D& operator=(const CartoonSphere2D&) = delete;
     244           0 :   CartoonSphere2D& operator=(CartoonSphere2D&&) = default;
     245           0 :   ~CartoonSphere2D() override = default;
     246             : 
     247           0 :   Domain<3> create_domain() const override;
     248             : 
     249             :   std::vector<DirectionMap<
     250             :       3, std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>>>
     251           1 :   external_boundary_conditions() const override;
     252             : 
     253           1 :   std::vector<std::array<size_t, 3>> initial_extents() const override;
     254             : 
     255           1 :   std::vector<std::array<size_t, 3>> initial_refinement_levels() const override;
     256             : 
     257             :   /// The block names are Shell0_LowerY, Shell0_UpperX, Shell0_UpperY,
     258             :   /// Shell1_LowerY, and so on. The naming and numbering goes from outer-most
     259             :   /// shell, starting from bottom and going counterclockwise, going to inside
     260             :   /// neighboring shell. If the center is included, the "center" half-circle
     261             :   /// follows same numbering with the _HalfSquare being the last block.
     262           1 :   std::vector<std::string> block_names() const override { return block_names_; }
     263             : 
     264             :   /// The block groups are Shell0, Shell1, ... starting from the outermost
     265             :   /// partition and working in.
     266             :   std::unordered_map<std::string, std::unordered_set<std::string>>
     267           1 :   block_groups() const override {
     268             :     return block_groups_;
     269             :   }
     270             : 
     271           1 :   auto functions_of_time(const std::unordered_map<std::string, double>&
     272             :                              initial_expiration_times = {}) const
     273             :       -> std::unordered_map<
     274             :           std::string,
     275             :           std::unique_ptr<domain::FunctionsOfTime::FunctionOfTime>> override;
     276             : 
     277             :  private:
     278           0 :   double inner_radius_{};
     279           0 :   double outer_radius_{};
     280           0 :   size_t initial_angular_refinement_{};
     281           0 :   std::vector<size_t> initial_radial_refinement_{};
     282           0 :   std::array<size_t, 2> initial_number_of_grid_points_{};
     283           0 :   std::vector<double> radial_partitioning_{};
     284           0 :   bool use_equiangular_map_{false};
     285           0 :   std::variant<Excision, InnerSquare> interior_{};
     286           0 :   bool fill_interior_{false};
     287             :   std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
     288           0 :       time_dependence_ = nullptr;
     289           0 :   std::vector<CoordinateMaps::Distribution> radial_distributions_{};
     290             :   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
     291           0 :       outer_boundary_condition_ = nullptr;
     292             :   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
     293           0 :       cartoon_boundary_condition_ = nullptr;
     294           0 :   std::vector<std::string> block_names_;
     295             :   std::unordered_map<std::string, std::unordered_set<std::string>>
     296           0 :       block_groups_;
     297           0 :   size_t num_shells_{};
     298           0 :   size_t num_blocks_{};
     299             : };
     300             : }  // namespace domain::creators
     301             : 
     302             : namespace domain::creators::detail {
     303             : /// \brief Helper struct for CartoonSphere2D options parsing so the internal
     304             : /// cartoon boundary condition at $x = 0$ is not a required argument.
     305             : ///
     306             : /// \details To get the cartoon-type boundary condition from a system we need
     307             : /// access to the system's metavariables, which is only present when parsing
     308             : /// options. The point of this design is so that the input file does not require
     309             : /// a dummy value for an InnerBoundaryCondition that is always the same for a
     310             : /// given system.
     311             : ///
     312             : /// This helper's constructor does not take an InnerBoundaryCondition, which
     313             : /// means if we parse a CartoonSphere2D as a CartonSphere2DOptionsHelper
     314             : /// by having an options parsing specialization (so the input file values are
     315             : /// the helper's construtor, not CartoonSphere2D's), we can detect the
     316             : /// cartoon-type boundary condition for the given system and only then call the
     317             : /// CartoonSphere2D constructor with the extra information.
     318             : struct CartoonSphere2DOptionsHelper {
     319             :   // Inherit the options template from CartoonSphere2D
     320             :   template <typename Metavariables>
     321             :   using options = typename domain::creators::CartoonSphere2D::template options<
     322             :       Metavariables>;
     323             : 
     324             :   using InitialRadialRefinement =
     325             :       domain::creators::CartoonSphere2D::InitialRadialRefinement;
     326             :   using Interior = domain::creators::CartoonSphere2D::Interior;
     327             : 
     328             :   template <typename BoundaryConditionsBase>
     329             :   using OuterBoundaryCondition =
     330             :       domain::creators::CartoonSphere2D::OuterBoundaryCondition<
     331             :           BoundaryConditionsBase>;
     332             : 
     333             :   static constexpr Options::String help = {"CartoonSphere2DOptionsHelper."};
     334             : 
     335             :   // Default constructor required by Options system
     336             :   CartoonSphere2DOptionsHelper() = default;
     337             : 
     338             :   // Does not take cartoon BC
     339             :   CartoonSphere2DOptionsHelper(
     340             :       double inner_radius, double outer_radius,
     341             :       size_t initial_angular_refinement,
     342             :       typename InitialRadialRefinement::type&& initial_radial_refinement,
     343             :       std::array<size_t, 2> initial_number_of_grid_points,
     344             :       std::vector<double> radial_partitioning, bool use_equiangular_map,
     345             :       std::variant<Excision, InnerSquare> interior,
     346             :       typename domain::creators::CartoonSphere2D::RadialDistribution::type
     347             :           radial_distribution = CoordinateMaps::Distribution::Linear,
     348             :       std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
     349             :           time_dependence = nullptr,
     350             :       std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
     351             :           outer_boundary_condition = nullptr,
     352             :       Options::Context&& context = {});
     353             : 
     354             :   // Same members as in CartoonSphere2D; public, to be extracted
     355             :   double inner_radius_{std::numeric_limits<double>::signaling_NaN()};
     356             :   double outer_radius_{std::numeric_limits<double>::signaling_NaN()};
     357             :   size_t initial_angular_refinement_{0};
     358             :   typename InitialRadialRefinement::type initial_radial_refinement_{};
     359             :   std::array<size_t, 2> initial_number_of_grid_points_{};
     360             :   std::vector<double> radial_partitioning_{};
     361             :   bool use_equiangular_map_{false};
     362             :   typename Interior::type interior_{};
     363             :   typename domain::creators::CartoonSphere2D::RadialDistribution::type
     364             :       radial_distribution_{CoordinateMaps::Distribution::Linear};
     365             :   std::unique_ptr<domain::creators::time_dependence::TimeDependence<3>>
     366             :       time_dependence_{nullptr};
     367             :   std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
     368             :       outer_boundary_condition_{nullptr};
     369             :   Options::Context context_{};
     370             : };
     371             : }  // namespace domain::creators::detail
     372             : 
     373             : // Options parsing specialization to automate CartoonSphere2D's cartoon
     374             : // boundary condition
     375             : template <>
     376           0 : struct Options::create_from_yaml<domain::creators::CartoonSphere2D> {
     377             :   template <typename Metavariables>
     378           0 :   static domain::creators::CartoonSphere2D create(
     379             :       const Options::Option& options) {
     380             :     auto helper =
     381             :         options.parse_as<domain::creators::detail::CartoonSphere2DOptionsHelper,
     382             :                          Metavariables>();
     383             : 
     384             :     // Create cartoon BC if system supports it, if not will throw parse error in
     385             :     // real constructor
     386             :     std::unique_ptr<domain::BoundaryConditions::BoundaryCondition>
     387             :         cartoon_boundary_condition = nullptr;
     388             :     if constexpr (domain::BoundaryConditions::has_boundary_conditions_base_v<
     389             :                       typename Metavariables::system>) {
     390             :       if constexpr (domain::BoundaryConditions::system_has_cartoon_bc_v<
     391             :                         Metavariables>) {
     392             :         cartoon_boundary_condition =
     393             :             domain::BoundaryConditions::make_cartoon_boundary_condition<
     394             :                 Metavariables>();
     395             :       }
     396             :     }
     397             : 
     398             :     // Construct CartoonSphere2D using the helper's parsed data + cartoon BC
     399             :     return domain::creators::CartoonSphere2D(
     400             :         helper.inner_radius_, helper.outer_radius_,
     401             :         helper.initial_angular_refinement_,
     402             :         std::move(helper.initial_radial_refinement_),
     403             :         std::move(helper.initial_number_of_grid_points_),
     404             :         std::move(helper.radial_partitioning_), helper.use_equiangular_map_,
     405             :         std::move(helper.interior_), helper.radial_distribution_,
     406             :         std::move(helper.time_dependence_),
     407             :         std::move(helper.outer_boundary_condition_),
     408             :         std::move(cartoon_boundary_condition), helper.context_);
     409             :   }
     410             : };