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