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             :   /// Method to use for reconstructing the DG solution from the subcell
     102             :   /// solution.
     103           1 :   struct SubcellToDgReconstructionMethod {
     104           0 :     static constexpr Options::String help{
     105             :         "Method to use for reconstructing the DG solution from the subcell "
     106             :         "solution."};
     107           0 :     using type = fd::ReconstructionMethod;
     108             :   };
     109             :   /// \brief Use a width-one halo of FD elements around any troubled element.
     110             :   ///
     111             :   /// This provides a buffer of FD subcells so that as a discontinuity moves
     112             :   /// from one element to another we do not get any Gibbs phenomenon. In the
     113             :   /// case where we evolve the spacetime metric (e.g. GH+GRMHD) a halo region
     114             :   /// provides a buffer in case the stellar surface is near an element
     115             :   /// boundary. Since the GH variables are interpolated using high-order
     116             :   /// unlimited reconstruction, they can run into issues with Gibbs phenomenon.
     117           1 :   struct UseHalo {
     118           0 :     using type = bool;
     119           0 :     static constexpr Options::String help = {
     120             :         "Use a width-one halo of FD elements around any troubled element."
     121             :         "\n"
     122             :         "This provides a buffer of FD subcells so that as a discontinuity "
     123             :         "moves from one element to another we do not get any Gibbs "
     124             :         "phenomenon."};
     125           0 :     using group = TroubledCellIndicator;
     126             :   };
     127             : 
     128             :   /// \brief A list of block names on which to never do subcell.
     129             :   ///
     130             :   /// Set to `None` to allow subcell in all blocks.
     131           1 :   struct OnlyDgBlocksAndGroups {
     132           0 :     using type =
     133             :         Options::Auto<std::vector<std::string>, Options::AutoLabel::None>;
     134           0 :     static constexpr Options::String help = {
     135             :         "A list of block and group names on which to never do subcell.\n"
     136             :         "Set to 'None' to not restrict where FD can be used."};
     137           0 :     using group = TroubledCellIndicator;
     138             :   };
     139             : 
     140             :   /// \brief The order of the FD derivative used.
     141             :   ///
     142             :   /// Must be one of 2, 4, 6, 8, or 10. If `Auto` then the derivative order is
     143             :   /// determined for you, typically the next-lowest even order compared with
     144             :   /// the reconstruction scheme order. E.g. for a 5th-order reconstruction we
     145             :   /// would use 4th order derivatives.
     146           1 :   struct FiniteDifferenceDerivativeOrder {
     147           0 :     using type = ::fd::DerivativeOrder;
     148           0 :     static constexpr Options::String help = {
     149             :         "The finite difference derivative order to use. If computed from the "
     150             :         "reconstruction, then the reconstruction method must support returning "
     151             :         "its reconstruction order."};
     152             :   };
     153             : 
     154           0 :   struct FdToDgTci {
     155           0 :     static constexpr Options::String help =
     156             :         "Options related to how quickly we switch from FD to DG.";
     157           0 :     using group = TroubledCellIndicator;
     158             :   };
     159             :   /// The number of time steps taken between calls to the TCI to check if we
     160             :   /// can go back to the DG grid. A value of `1` means every time step, while
     161             :   /// `2` means every other time step.
     162           1 :   struct NumberOfStepsBetweenTciCalls {
     163           0 :     static constexpr Options::String help{
     164             :         "The number of time steps taken between calls to the TCI to check if "
     165             :         "we can go back to the DG grid. A value of `1` means every time step, "
     166             :         "while `2` means every other time step."};
     167           0 :     using type = size_t;
     168           0 :     static constexpr type lower_bound() { return 1; }
     169           0 :     using group = FdToDgTci;
     170             :   };
     171             :   /// The number of time steps/TCI calls after a switch from DG to FD before we
     172             :   /// allow switching back to DG.
     173           1 :   struct MinTciCallsAfterRollback {
     174           0 :     static constexpr Options::String help{
     175             :         "The number of time steps/TCI calls after a switch from DG to FD "
     176             :         "before we allow switching back to DG."};
     177           0 :     using type = size_t;
     178           0 :     static constexpr type lower_bound() { return 1; }
     179           0 :     using group = FdToDgTci;
     180             :   };
     181             :   /// The number of time steps/TCI calls that the TCI needs to have decided
     182             :   /// switching to DG is fine before we actually do the switch.
     183           1 :   struct MinimumClearTcis {
     184           0 :     static constexpr Options::String help{
     185             :         "The number of time steps/TCI calls that the TCI needs to have decided "
     186             :         "switching to DG is fine before we actually do the switch."};
     187           0 :     using type = size_t;
     188           0 :     static constexpr type lower_bound() { return 1; }
     189           0 :     using group = FdToDgTci;
     190             :   };
     191             : 
     192           0 :   using options = tmpl::list<
     193             :       PerssonExponent, PerssonNumHighestModes, RdmpDelta0, RdmpEpsilon,
     194             :       AlwaysUseSubcells, SubcellToDgReconstructionMethod, UseHalo,
     195             :       OnlyDgBlocksAndGroups, FiniteDifferenceDerivativeOrder,
     196             :       NumberOfStepsBetweenTciCalls, MinTciCallsAfterRollback, MinimumClearTcis>;
     197             : 
     198           0 :   static constexpr Options::String help{
     199             :       "System-agnostic options for the DG-subcell method."};
     200             : 
     201           0 :   SubcellOptions() = default;
     202           0 :   SubcellOptions(
     203             :       double persson_exponent, size_t persson_num_highest_modes,
     204             :       double rdmp_delta0, double rdmp_epsilon, bool always_use_subcells,
     205             :       fd::ReconstructionMethod recons_method, bool use_halo,
     206             :       std::optional<std::vector<std::string>> only_dg_block_and_group_names,
     207             :       ::fd::DerivativeOrder finite_difference_derivative_order,
     208             :       size_t number_of_steps_between_tci_calls,
     209             :       size_t min_tci_calls_after_rollback, size_t min_clear_tci_before_dg);
     210             : 
     211             :   /// \brief Given an existing SubcellOptions that was created from block and
     212             :   /// group names, create one that stores block IDs.
     213             :   ///
     214             :   /// The `DomainCreator` is used to convert block and group names into IDs
     215             :   /// and also to check that all listed block names and groups are in the
     216             :   /// domain.
     217             :   ///
     218             :   /// \note This is a workaround since our option parser does not allow us to
     219             :   /// retrieve options specified somewhere completely different in the input
     220             :   /// file.
     221             :   template <size_t Dim>
     222           1 :   SubcellOptions(const SubcellOptions& subcell_options_with_block_names,
     223             :                  const DomainCreator<Dim>& domain_creator);
     224             : 
     225           0 :   void pup(PUP::er& p);
     226             : 
     227           0 :   double persson_exponent() const { return persson_exponent_; }
     228             : 
     229           0 :   size_t persson_num_highest_modes() const {
     230             :     return persson_num_highest_modes_;
     231             :   }
     232             : 
     233           0 :   double rdmp_delta0() const { return rdmp_delta0_; }
     234             : 
     235           0 :   double rdmp_epsilon() const { return rdmp_epsilon_; }
     236             : 
     237           0 :   bool always_use_subcells() const { return always_use_subcells_; }
     238             : 
     239           0 :   fd::ReconstructionMethod reconstruction_method() const {
     240             :     return reconstruction_method_;
     241             :   }
     242             : 
     243           0 :   bool use_halo() const { return use_halo_; }
     244             : 
     245           0 :   const std::vector<size_t>& only_dg_block_ids() const {
     246             :     ASSERT(only_dg_block_ids_.has_value(),
     247             :            "The block IDs on which we are only allowed to do DG have not been "
     248             :            "set.");
     249             :     return only_dg_block_ids_.value();
     250             :   }
     251             : 
     252           0 :   ::fd::DerivativeOrder finite_difference_derivative_order() const {
     253             :     return finite_difference_derivative_order_;
     254             :   }
     255             : 
     256             :   /// The number of time steps between when we check the TCI to see if we can
     257             :   /// switch from FD back to DG.
     258           1 :   size_t number_of_steps_between_tci_calls() const {
     259             :     return number_of_steps_between_tci_calls_;
     260             :   }
     261             : 
     262             :   /// The number of time steps after a rollback before we check the TCI to see
     263             :   /// if we can switch from FD back to DG.
     264           1 :   size_t min_tci_calls_after_rollback() const {
     265             :     return min_tci_calls_after_rollback_;
     266             :   }
     267             : 
     268             :   /// The number of times the TCI must have flagged that we can switch back to
     269             :   /// DG before doing so.
     270             :   ///
     271             :   /// Note that if we only check the TCI every
     272             :   /// `number_of_steps_between_tci_calls()` then it takes
     273             :   /// `number_of_steps_between_tci_calls * min_clear_tci_before_dg` time steps
     274             :   /// before switching from FD back to DG.
     275             :   ///
     276             :   /// `0 means
     277           1 :   size_t min_clear_tci_before_dg() const { return min_clear_tci_before_dg_; }
     278             : 
     279             :  private:
     280           0 :   friend bool operator==(const SubcellOptions& lhs, const SubcellOptions& rhs);
     281             : 
     282           0 :   double persson_exponent_ = std::numeric_limits<double>::signaling_NaN();
     283           0 :   size_t persson_num_highest_modes_ =
     284             :       std::numeric_limits<size_t>::signaling_NaN();
     285           0 :   double rdmp_delta0_ = std::numeric_limits<double>::signaling_NaN();
     286           0 :   double rdmp_epsilon_ = std::numeric_limits<double>::signaling_NaN();
     287           0 :   bool always_use_subcells_ = false;
     288           0 :   fd::ReconstructionMethod reconstruction_method_ =
     289             :       fd::ReconstructionMethod::AllDimsAtOnce;
     290           0 :   bool use_halo_{false};
     291           0 :   std::optional<std::vector<std::string>> only_dg_block_and_group_names_{};
     292           0 :   std::optional<std::vector<size_t>> only_dg_block_ids_{};
     293           0 :   ::fd::DerivativeOrder finite_difference_derivative_order_{};
     294           0 :   size_t number_of_steps_between_tci_calls_{1};
     295           0 :   size_t min_tci_calls_after_rollback_{1};
     296           0 :   size_t min_clear_tci_before_dg_{0};
     297             : };
     298             : 
     299           0 : bool operator!=(const SubcellOptions& lhs, const SubcellOptions& rhs);
     300             : }  // namespace evolution::dg::subcell