Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <limits>
       8             : #include <optional>
       9             : #include <string>
      10             : #include <vector>
      11             : 
      12             : #include "Evolution/DgSubcell/ReconstructionMethod.hpp"
      13             : #include "NumericalAlgorithms/FiniteDifference/DerivativeOrder.hpp"
      14             : #include "Options/Auto.hpp"
      15             : #include "Options/Context.hpp"
      16             : #include "Options/String.hpp"
      17             : #include "Utilities/ErrorHandling/Assert.hpp"
      18             : #include "Utilities/Literals.hpp"
      19             : #include "Utilities/TMPL.hpp"
      20             : 
      21             : /// \cond
      22             : template <size_t VolumeDim>
      23             : class DomainCreator;
      24             : namespace PUP {
      25             : class er;
      26             : }  // namespace PUP
      27             : /// \endcond
      28             : 
      29             : namespace evolution::dg::subcell {
      30             : /*!
      31             :  * \brief Holds the system-agnostic subcell parameters, such as numbers
      32             :  * controlling when to switch between DG and subcell.
      33             :  */
      34           1 : class SubcellOptions {
      35             :  public:
      36             :   /// Parameters related to the troubled cell indicator (TCI) that determines
      37             :   /// when to switch between DG and FD.
      38           1 :   struct TroubledCellIndicator {
      39           0 :     static constexpr Options::String help =
      40             :         "Parameters related to the troubled cell indicator (TCI) that "
      41             :         "determines when to switch between DG and FD.";
      42             :   };
      43             : 
      44           0 :   struct PerssonTci {
      45           0 :     static constexpr Options::String help =
      46             :         "Parameters related to the Persson TCI";
      47           0 :     using group = TroubledCellIndicator;
      48             :   };
      49             :   /// The exponent \f$\alpha\f$ passed to the Persson troubled-cell indicator
      50           1 :   struct PerssonExponent {
      51           0 :     static std::string name() { return "Exponent"; }
      52           0 :     static constexpr Options::String help{
      53             :         "The exponent at which the error should decrease with (N+1-M)"};
      54           0 :     using type = double;
      55           0 :     static constexpr type lower_bound() { return 1.0; }
      56           0 :     static constexpr type upper_bound() { return 10.0; }
      57           0 :     using group = PerssonTci;
      58             :   };
      59             :   /// The number of highest modes the Persson troubled-cell indicator monitors
      60           1 :   struct PerssonNumHighestModes {
      61           0 :     static std::string name() { return "NumHighestModes"; }
      62           0 :     static constexpr Options::String help{
      63             :         "The number of highest modes M the Persson TCI monitors."};
      64           0 :     using type = size_t;
      65           0 :     static constexpr type lower_bound() { return 1_st; }
      66           0 :     static constexpr type upper_bound() { return 10_st; }
      67           0 :     using group = PerssonTci;
      68             :   };
      69             : 
      70           0 :   struct RdmpTci {
      71           0 :     static constexpr Options::String help =
      72             :         "Parameters related to the relaxed discrete maximum principle TCI";
      73           0 :     using group = TroubledCellIndicator;
      74             :   };
      75             :   /// The \f$\delta_0\f$ parameter in the relaxed discrete maximum principle
      76             :   /// troubled-cell indicator
      77           1 :   struct RdmpDelta0 {
      78           0 :     static std::string name() { return "Delta0"; }
      79           0 :     static constexpr Options::String help{"Absolute jump tolerance parameter."};
      80           0 :     using type = double;
      81           0 :     static type lower_bound() { return 0.0; }
      82           0 :     using group = RdmpTci;
      83             :   };
      84             :   /// The \f$\epsilon\f$ parameter in the relaxed discrete maximum principle
      85             :   /// troubled-cell indicator
      86           1 :   struct RdmpEpsilon {
      87           0 :     static std::string name() { return "Epsilon"; }
      88           0 :     static constexpr Options::String help{
      89             :         "The jump-dependent relaxation constant."};
      90           0 :     using type = double;
      91           0 :     static type lower_bound() { return 0.0; }
      92           0 :     static type upper_bound() { return 1.0; }
      93           0 :     using group = RdmpTci;
      94             :   };
      95             :   /// If true, then we always use the subcell method, not DG.
      96           1 :   struct AlwaysUseSubcells {
      97           0 :     static constexpr Options::String help{
      98             :         "If true, then always use the subcell method (e.g. finite-difference) "
      99             :         "instead of DG."};
     100           0 :     using type = bool;
     101           0 :     using group = TroubledCellIndicator;
     102             :   };
     103             :   /// If true, then we allow extension directions to be used for
     104             :   /// computing the ghost points in the FD subcell method to avoid
     105             :   /// extrapolation.
     106           1 :   struct EnableExtensionDirections {
     107           0 :     static constexpr Options::String help{
     108             :         "If true, then we allow extension directions to be used for "
     109             :         "computing the ghost points in the FD subcell method to avoid "
     110             :         "extrapolation."};
     111           0 :     using type = bool;
     112           0 :     using group = TroubledCellIndicator;
     113             :   };
     114             :   /// Method to use for reconstructing the DG solution from the subcell
     115             :   /// solution.
     116           1 :   struct SubcellToDgReconstructionMethod {
     117           0 :     static constexpr Options::String help{
     118             :         "Method to use for reconstructing the DG solution from the subcell "
     119             :         "solution."};
     120           0 :     using type = fd::ReconstructionMethod;
     121             :   };
     122             :   /// \brief Use a width-one halo of FD elements around any troubled element.
     123             :   ///
     124             :   /// This provides a buffer of FD subcells so that as a discontinuity moves
     125             :   /// from one element to another we do not get any Gibbs phenomenon. In the
     126             :   /// case where we evolve the spacetime metric (e.g. GH+GRMHD) a halo region
     127             :   /// provides a buffer in case the stellar surface is near an element
     128             :   /// boundary. Since the GH variables are interpolated using high-order
     129             :   /// unlimited reconstruction, they can run into issues with Gibbs phenomenon.
     130           1 :   struct UseHalo {
     131           0 :     using type = bool;
     132           0 :     static constexpr Options::String help = {
     133             :         "Use a width-one halo of FD elements around any troubled element."
     134             :         "\n"
     135             :         "This provides a buffer of FD subcells so that as a discontinuity "
     136             :         "moves from one element to another we do not get any Gibbs "
     137             :         "phenomenon."};
     138           0 :     using group = TroubledCellIndicator;
     139             :   };
     140             : 
     141             :   /// \brief A list of block names on which to never do subcell.
     142             :   ///
     143             :   /// Set to `None` to allow subcell in all blocks.
     144           1 :   struct OnlyDgBlocksAndGroups {
     145           0 :     using type =
     146             :         Options::Auto<std::vector<std::string>, Options::AutoLabel::None>;
     147           0 :     static constexpr Options::String help = {
     148             :         "A list of block and group names on which to never do subcell.\n"
     149             :         "Set to 'None' to not restrict where FD can be used."};
     150           0 :     using group = TroubledCellIndicator;
     151             :   };
     152             : 
     153             :   /// \brief The order of the FD derivative used.
     154             :   ///
     155             :   /// Must be one of 2, 4, 6, 8, or 10. If `Auto` then the derivative order is
     156             :   /// determined for you, typically the next-lowest even order compared with
     157             :   /// the reconstruction scheme order. E.g. for a 5th-order reconstruction we
     158             :   /// would use 4th order derivatives.
     159           1 :   struct FiniteDifferenceDerivativeOrder {
     160           0 :     using type = ::fd::DerivativeOrder;
     161           0 :     static constexpr Options::String help = {
     162             :         "The finite difference derivative order to use. If computed from the "
     163             :         "reconstruction, then the reconstruction method must support returning "
     164             :         "its reconstruction order."};
     165             :   };
     166             : 
     167           0 :   struct FdToDgTci {
     168           0 :     static constexpr Options::String help =
     169             :         "Options related to how quickly we switch from FD to DG.";
     170           0 :     using group = TroubledCellIndicator;
     171             :   };
     172             :   /// The number of time steps taken between calls to the TCI to check if we
     173             :   /// can go back to the DG grid. A value of `1` means every time step, while
     174             :   /// `2` means every other time step.
     175           1 :   struct NumberOfStepsBetweenTciCalls {
     176           0 :     static constexpr Options::String help{
     177             :         "The number of time steps taken between calls to the TCI to check if "
     178             :         "we can go back to the DG grid. A value of `1` means every time step, "
     179             :         "while `2` means every other time step."};
     180           0 :     using type = size_t;
     181           0 :     static constexpr type lower_bound() { return 1; }
     182           0 :     using group = FdToDgTci;
     183             :   };
     184             :   /// The number of time steps/TCI calls after a switch from DG to FD before we
     185             :   /// allow switching back to DG.
     186           1 :   struct MinTciCallsAfterRollback {
     187           0 :     static constexpr Options::String help{
     188             :         "The number of time steps/TCI calls after a switch from DG to FD "
     189             :         "before we allow switching back to DG."};
     190           0 :     using type = size_t;
     191           0 :     static constexpr type lower_bound() { return 1; }
     192           0 :     using group = FdToDgTci;
     193             :   };
     194             :   /// The number of time steps/TCI calls that the TCI needs to have decided
     195             :   /// switching to DG is fine before we actually do the switch.
     196           1 :   struct MinimumClearTcis {
     197           0 :     static constexpr Options::String help{
     198             :         "The number of time steps/TCI calls that the TCI needs to have decided "
     199             :         "switching to DG is fine before we actually do the switch."};
     200           0 :     using type = size_t;
     201           0 :     static constexpr type lower_bound() { return 1; }
     202           0 :     using group = FdToDgTci;
     203             :   };
     204             : 
     205             :   /// \brief FD to FD interpolation order.
     206             :   ///
     207           1 :   struct FdInterpolationOrder {
     208           0 :     using type = size_t;
     209           0 :     static constexpr Options::String help = {
     210             :         "The interpolation order to use when interpolating from FD block to FD "
     211             :         "block."};
     212           0 :     static constexpr type lower_bound() { return 1; }
     213           0 :     static constexpr type upper_bound() { return 3; }
     214             :   };
     215             : 
     216             :   /// \brief Subcell steps per slab in regions allowed to do subcell
     217             :   /// in LTS mode.  Must be a power of 2.
     218           1 :   struct LtsStepsPerSlab {
     219           0 :     using type = size_t;
     220           0 :     static constexpr Options::String help = {
     221             :         "Steps per slab in regions allowed to do subcell.  Must be a power of "
     222             :         "2."};
     223           0 :     static constexpr type lower_bound() { return 1; }
     224           0 :     static constexpr type upper_bound() { return 1_st << 31; }
     225             :   };
     226             : 
     227           0 :   using gts_options =
     228             :       tmpl::list<PerssonExponent, PerssonNumHighestModes, RdmpDelta0,
     229             :                  RdmpEpsilon, AlwaysUseSubcells, EnableExtensionDirections,
     230             :                  SubcellToDgReconstructionMethod, UseHalo,
     231             :                  OnlyDgBlocksAndGroups, FiniteDifferenceDerivativeOrder,
     232             :                  NumberOfStepsBetweenTciCalls, MinTciCallsAfterRollback,
     233             :                  MinimumClearTcis, FdInterpolationOrder>;
     234             : 
     235             :   template <typename Metavars>
     236           0 :   using options =
     237             :       tmpl::conditional_t<Metavars::local_time_stepping,
     238             :                           tmpl::push_back<gts_options, LtsStepsPerSlab>,
     239             :                           gts_options>;
     240             : 
     241           0 :   static constexpr Options::String help{
     242             :       "System-agnostic options for the DG-subcell method."};
     243             : 
     244           0 :   SubcellOptions() = default;
     245           0 :   SubcellOptions(
     246             :       double persson_exponent, size_t persson_num_highest_modes,
     247             :       double rdmp_delta0, double rdmp_epsilon, bool always_use_subcells,
     248             :       bool enable_extension_directions, fd::ReconstructionMethod recons_method,
     249             :       bool use_halo,
     250             :       std::optional<std::vector<std::string>> only_dg_block_and_group_names,
     251             :       ::fd::DerivativeOrder finite_difference_derivative_order,
     252             :       size_t number_of_steps_between_tci_calls,
     253             :       size_t min_tci_calls_after_rollback, size_t min_clear_tci_before_dg,
     254             :       size_t fd_to_fd_interp_order = 1_st,
     255             :       const std::optional<size_t>& lts_steps_per_slab = std::nullopt,
     256             :       const Options::Context& context = {});
     257             : 
     258             :   /// \brief Given an existing SubcellOptions that was created from block and
     259             :   /// group names, create one that stores block IDs.
     260             :   ///
     261             :   /// The `DomainCreator` is used to convert block and group names into IDs
     262             :   /// and also to check that all listed block names and groups are in the
     263             :   /// domain.
     264             :   ///
     265             :   /// \note This is a workaround since our option parser does not allow us to
     266             :   /// retrieve options specified somewhere completely different in the input
     267             :   /// file.
     268             :   template <size_t Dim>
     269           1 :   SubcellOptions(const SubcellOptions& subcell_options_with_block_names,
     270             :                  const DomainCreator<Dim>& domain_creator);
     271             : 
     272           0 :   void pup(PUP::er& p);
     273             : 
     274           0 :   double persson_exponent() const { return persson_exponent_; }
     275             : 
     276           0 :   size_t persson_num_highest_modes() const {
     277             :     return persson_num_highest_modes_;
     278             :   }
     279             : 
     280           0 :   double rdmp_delta0() const { return rdmp_delta0_; }
     281             : 
     282           0 :   double rdmp_epsilon() const { return rdmp_epsilon_; }
     283             : 
     284           0 :   bool always_use_subcells() const { return always_use_subcells_; }
     285             : 
     286           0 :   bool enable_extension_directions() const {
     287             :     return enable_extension_directions_;
     288             :   }
     289             : 
     290           0 :   fd::ReconstructionMethod reconstruction_method() const {
     291             :     return reconstruction_method_;
     292             :   }
     293             : 
     294           0 :   bool use_halo() const { return use_halo_; }
     295             : 
     296           0 :   const std::vector<size_t>& only_dg_block_ids() const {
     297             :     ASSERT(only_dg_block_ids_.has_value(),
     298             :            "The block IDs on which we are only allowed to do DG have not been "
     299             :            "set.");
     300             :     return only_dg_block_ids_.value();
     301             :   }
     302             : 
     303           0 :   size_t get_fd_to_fd_interp_order() const { return fd_to_fd_interp_order_; }
     304             : 
     305           0 :   ::fd::DerivativeOrder finite_difference_derivative_order() const {
     306             :     return finite_difference_derivative_order_;
     307             :   }
     308             : 
     309             :   /// The number of time steps between when we check the TCI to see if we can
     310             :   /// switch from FD back to DG.
     311           1 :   size_t number_of_steps_between_tci_calls() const {
     312             :     return number_of_steps_between_tci_calls_;
     313             :   }
     314             : 
     315             :   /// The number of time steps after a rollback before we check the TCI to see
     316             :   /// if we can switch from FD back to DG.
     317           1 :   size_t min_tci_calls_after_rollback() const {
     318             :     return min_tci_calls_after_rollback_;
     319             :   }
     320             : 
     321             :   /// The number of times the TCI must have flagged that we can switch back to
     322             :   /// DG before doing so.
     323             :   ///
     324             :   /// Note that if we only check the TCI every
     325             :   /// `number_of_steps_between_tci_calls()` then it takes
     326             :   /// `number_of_steps_between_tci_calls * min_clear_tci_before_dg` time steps
     327             :   /// before switching from FD back to DG.
     328             :   ///
     329             :   /// `0 means
     330           1 :   size_t min_clear_tci_before_dg() const { return min_clear_tci_before_dg_; }
     331             : 
     332             :   /// The number of steps per slab in regions allowed to do subcell.
     333             :   /// Only valid in LTS-subcell evolutions.
     334           1 :   size_t lts_steps_per_slab() const {
     335             :     ASSERT(lts_steps_per_slab_.has_value(),
     336             :            "Attempted to access lts_steps_per_slab in GTS evolution.");
     337             :     return lts_steps_per_slab_.value();
     338             :   }
     339             : 
     340             :  private:
     341           0 :   friend bool operator==(const SubcellOptions& lhs, const SubcellOptions& rhs);
     342             : 
     343           0 :   double persson_exponent_ = std::numeric_limits<double>::signaling_NaN();
     344           0 :   size_t persson_num_highest_modes_ =
     345             :       std::numeric_limits<size_t>::signaling_NaN();
     346           0 :   double rdmp_delta0_ = std::numeric_limits<double>::signaling_NaN();
     347           0 :   double rdmp_epsilon_ = std::numeric_limits<double>::signaling_NaN();
     348           0 :   bool always_use_subcells_ = false;
     349           0 :   bool enable_extension_directions_ = false;
     350           0 :   fd::ReconstructionMethod reconstruction_method_ =
     351             :       fd::ReconstructionMethod::AllDimsAtOnce;
     352           0 :   bool use_halo_{false};
     353           0 :   std::optional<std::vector<std::string>> only_dg_block_and_group_names_{};
     354           0 :   std::optional<std::vector<size_t>> only_dg_block_ids_{};
     355           0 :   ::fd::DerivativeOrder finite_difference_derivative_order_{};
     356           0 :   size_t number_of_steps_between_tci_calls_{1};
     357           0 :   size_t min_tci_calls_after_rollback_{1};
     358           0 :   size_t min_clear_tci_before_dg_{0};
     359           0 :   size_t fd_to_fd_interp_order_{};
     360           0 :   std::optional<size_t> lts_steps_per_slab_{};
     361             : };
     362             : 
     363           0 : bool operator!=(const SubcellOptions& lhs, const SubcellOptions& rhs);
     364             : }  // namespace evolution::dg::subcell