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             :   /// \brief FD to FD interpolation order.
     204             :   ///
     205           1 :   struct FdInterpolationOrder {
     206           0 :     using type = size_t;
     207           0 :     static constexpr Options::String help = {
     208             :         "The interpolation order to use when interpolating from FD block to FD "
     209             :         "block."};
     210           0 :     static constexpr type lower_bound() { return 1; }
     211           0 :     static constexpr type upper_bound() { return 3; }
     212             :   };
     213             : 
     214           0 :   using options =
     215             :       tmpl::list<PerssonExponent, PerssonNumHighestModes, RdmpDelta0,
     216             :                  RdmpEpsilon, AlwaysUseSubcells, EnableExtensionDirections,
     217             :                  SubcellToDgReconstructionMethod, UseHalo,
     218             :                  OnlyDgBlocksAndGroups, FiniteDifferenceDerivativeOrder,
     219             :                  NumberOfStepsBetweenTciCalls, MinTciCallsAfterRollback,
     220             :                  MinimumClearTcis, FdInterpolationOrder>;
     221             : 
     222           0 :   static constexpr Options::String help{
     223             :       "System-agnostic options for the DG-subcell method."};
     224             : 
     225           0 :   SubcellOptions() = default;
     226           0 :   SubcellOptions(
     227             :       double persson_exponent, size_t persson_num_highest_modes,
     228             :       double rdmp_delta0, double rdmp_epsilon, bool always_use_subcells,
     229             :       bool enable_extension_directions, fd::ReconstructionMethod recons_method,
     230             :       bool use_halo,
     231             :       std::optional<std::vector<std::string>> only_dg_block_and_group_names,
     232             :       ::fd::DerivativeOrder finite_difference_derivative_order,
     233             :       size_t number_of_steps_between_tci_calls,
     234             :       size_t min_tci_calls_after_rollback, size_t min_clear_tci_before_dg,
     235             :       size_t fd_to_fd_interp_order = 1_st);
     236             : 
     237             :   /// \brief Given an existing SubcellOptions that was created from block and
     238             :   /// group names, create one that stores block IDs.
     239             :   ///
     240             :   /// The `DomainCreator` is used to convert block and group names into IDs
     241             :   /// and also to check that all listed block names and groups are in the
     242             :   /// domain.
     243             :   ///
     244             :   /// \note This is a workaround since our option parser does not allow us to
     245             :   /// retrieve options specified somewhere completely different in the input
     246             :   /// file.
     247             :   template <size_t Dim>
     248           1 :   SubcellOptions(const SubcellOptions& subcell_options_with_block_names,
     249             :                  const DomainCreator<Dim>& domain_creator);
     250             : 
     251           0 :   void pup(PUP::er& p);
     252             : 
     253           0 :   double persson_exponent() const { return persson_exponent_; }
     254             : 
     255           0 :   size_t persson_num_highest_modes() const {
     256             :     return persson_num_highest_modes_;
     257             :   }
     258             : 
     259           0 :   double rdmp_delta0() const { return rdmp_delta0_; }
     260             : 
     261           0 :   double rdmp_epsilon() const { return rdmp_epsilon_; }
     262             : 
     263           0 :   bool always_use_subcells() const { return always_use_subcells_; }
     264             : 
     265           0 :   bool enable_extension_directions() const {
     266             :     return enable_extension_directions_;
     267             :   }
     268             : 
     269           0 :   fd::ReconstructionMethod reconstruction_method() const {
     270             :     return reconstruction_method_;
     271             :   }
     272             : 
     273           0 :   bool use_halo() const { return use_halo_; }
     274             : 
     275           0 :   const std::vector<size_t>& only_dg_block_ids() const {
     276             :     ASSERT(only_dg_block_ids_.has_value(),
     277             :            "The block IDs on which we are only allowed to do DG have not been "
     278             :            "set.");
     279             :     return only_dg_block_ids_.value();
     280             :   }
     281             : 
     282           0 :   size_t get_fd_to_fd_interp_order() const { return fd_to_fd_interp_order_; }
     283             : 
     284           0 :   ::fd::DerivativeOrder finite_difference_derivative_order() const {
     285             :     return finite_difference_derivative_order_;
     286             :   }
     287             : 
     288             :   /// The number of time steps between when we check the TCI to see if we can
     289             :   /// switch from FD back to DG.
     290           1 :   size_t number_of_steps_between_tci_calls() const {
     291             :     return number_of_steps_between_tci_calls_;
     292             :   }
     293             : 
     294             :   /// The number of time steps after a rollback before we check the TCI to see
     295             :   /// if we can switch from FD back to DG.
     296           1 :   size_t min_tci_calls_after_rollback() const {
     297             :     return min_tci_calls_after_rollback_;
     298             :   }
     299             : 
     300             :   /// The number of times the TCI must have flagged that we can switch back to
     301             :   /// DG before doing so.
     302             :   ///
     303             :   /// Note that if we only check the TCI every
     304             :   /// `number_of_steps_between_tci_calls()` then it takes
     305             :   /// `number_of_steps_between_tci_calls * min_clear_tci_before_dg` time steps
     306             :   /// before switching from FD back to DG.
     307             :   ///
     308             :   /// `0 means
     309           1 :   size_t min_clear_tci_before_dg() const { return min_clear_tci_before_dg_; }
     310             : 
     311             :  private:
     312           0 :   friend bool operator==(const SubcellOptions& lhs, const SubcellOptions& rhs);
     313             : 
     314           0 :   double persson_exponent_ = std::numeric_limits<double>::signaling_NaN();
     315           0 :   size_t persson_num_highest_modes_ =
     316             :       std::numeric_limits<size_t>::signaling_NaN();
     317           0 :   double rdmp_delta0_ = std::numeric_limits<double>::signaling_NaN();
     318           0 :   double rdmp_epsilon_ = std::numeric_limits<double>::signaling_NaN();
     319           0 :   bool always_use_subcells_ = false;
     320           0 :   bool enable_extension_directions_ = false;
     321           0 :   fd::ReconstructionMethod reconstruction_method_ =
     322             :       fd::ReconstructionMethod::AllDimsAtOnce;
     323           0 :   bool use_halo_{false};
     324           0 :   std::optional<std::vector<std::string>> only_dg_block_and_group_names_{};
     325           0 :   std::optional<std::vector<size_t>> only_dg_block_ids_{};
     326           0 :   ::fd::DerivativeOrder finite_difference_derivative_order_{};
     327           0 :   size_t number_of_steps_between_tci_calls_{1};
     328           0 :   size_t min_tci_calls_after_rollback_{1};
     329           0 :   size_t min_clear_tci_before_dg_{0};
     330           0 :   size_t fd_to_fd_interp_order_{};
     331             : };
     332             : 
     333           0 : bool operator!=(const SubcellOptions& lhs, const SubcellOptions& rhs);
     334             : }  // namespace evolution::dg::subcell