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