Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : 
       8             : #include "DataStructures/Variables.hpp"
       9             : #include "Utilities/ErrorHandling/Assert.hpp"
      10             : #include "Utilities/TMPL.hpp"
      11             : 
      12             : namespace evolution::dg::subcell {
      13             : /*!
      14             :  * \brief Troubled cell indicator using a relaxed discrete maximum principle,
      15             :  * comparing the solution on two grids at the same point in time.
      16             :  *
      17             :  * Checks that the subcell solution \f$\underline{u}\f$ and the DG solution
      18             :  * \f$u\f$ satisfy
      19             :  *
      20             :  * \f{align*}{ \min(u)-\delta \le \underline{u} \le \max(u)+\delta \f}
      21             :  *
      22             :  * where
      23             :  *
      24             :  * \f{align*}{ \delta = \max\left[\delta_0, \epsilon(\max(u) - \min(u))\right]
      25             :  * \f}
      26             :  *
      27             :  * where \f$\delta_0\f$ and \f$\epsilon\f$ are constants controlling the maximum
      28             :  * absolute and relative change allowed when projecting the DG solution to the
      29             :  * subcell grid. We currently specify one value of \f$\delta_0\f$ and
      30             :  * \f$\epsilon\f$ for all variables, but this could be generalized to choosing
      31             :  * the allowed variation in a variable-specific manner.
      32             :  *
      33             :  * If all checks are passed and cell is not troubled, returns an integer `0`.
      34             :  * Otherwise returns 1-based index of the tag in the input Variables that fails
      35             :  * the check. For instance, if we have
      36             :  *
      37             :  *  - `Variables<tmpl::list<DgVar1, DgVar2, DgVar3>>` for `dg_evolved_vars`
      38             :  *  - `Variables<tmpl::list<SubVar1, SubVar2, SubVar3>>` for
      39             :  *    `subcell_evolved_vars`
      40             :  *
      41             :  * as inputs and TCI flags the second pair `DgVar2` and `SubVar2` not satisfying
      42             :  * two-mesh RDMP criteria, returned value is `2` since the second pair of tags
      43             :  * failed the check.
      44             :  *
      45             :  * \note Once a single pair of tags fails to satisfy the check, checks for the
      46             :  * remaining part of the input variables are skipped. In the example above, for
      47             :  * instance if the second pair (`DgVar2`,`SubVar2`) is flagged, the third pair
      48             :  * (`DgVar3`,`SubVar3`) is ignored and not checked.
      49             :  *
      50             :  */
      51             : template <typename... DgEvolvedVarsTags, typename... SubcellEvolvedVarsTags>
      52           1 : int two_mesh_rdmp_tci(
      53             :     const Variables<tmpl::list<DgEvolvedVarsTags...>>& dg_evolved_vars,
      54             :     const Variables<tmpl::list<SubcellEvolvedVarsTags...>>&
      55             :         subcell_evolved_vars,
      56             :     const double rdmp_delta0, const double rdmp_epsilon) {
      57             :   static_assert(sizeof...(DgEvolvedVarsTags) ==
      58             :                 sizeof...(SubcellEvolvedVarsTags));
      59             :   ASSERT(rdmp_delta0 > 0.0, "The RDMP delta0 parameter must be positive.");
      60             :   ASSERT(rdmp_epsilon > 0.0, "The RDMP epsilon parameter must be positive.");
      61             : 
      62             :   bool cell_is_troubled = false;
      63             :   int tci_status = 0;
      64             :   size_t tag_index = 0;
      65             : 
      66             :   tmpl::for_each<
      67             :       tmpl::list<tmpl::list<DgEvolvedVarsTags, SubcellEvolvedVarsTags>...>>(
      68             :       [&cell_is_troubled, &tag_index, &dg_evolved_vars, rdmp_delta0,
      69             :        rdmp_epsilon, &subcell_evolved_vars, &tci_status](auto tag_v) {
      70             :         if (cell_is_troubled) {
      71             :           return;
      72             :         }
      73             : 
      74             :         using tags_list = tmpl::type_from<decltype(tag_v)>;
      75             :         using tag = tmpl::front<tags_list>;
      76             :         using inactive_tag = tmpl::back<tags_list>;
      77             :         const auto& dg_var = get<tag>(dg_evolved_vars);
      78             :         const auto& subcell_var = get<inactive_tag>(subcell_evolved_vars);
      79             : 
      80             :         for (auto dg_it = dg_var.begin(), subcell_it = subcell_var.begin();
      81             :              dg_it != dg_var.end() and subcell_it != subcell_var.end();
      82             :              (void)++dg_it, (void)++subcell_it) {
      83             :           ASSERT(not cell_is_troubled,
      84             :                  "If a cell has already been marked as troubled during the "
      85             :                  "two mesh RDMP TCI, we should not be continuing to check "
      86             :                  "other variables.");
      87             :           using std::max;
      88             :           using std::min;
      89             : 
      90             :           const double max_dg = max(*dg_it);
      91             :           const double min_dg = min(*dg_it);
      92             :           const double max_subcell = max(*subcell_it);
      93             :           const double min_subcell = min(*subcell_it);
      94             :           const double delta =
      95             :               max(rdmp_delta0, rdmp_epsilon * (max_dg - min_dg));
      96             :           cell_is_troubled =
      97             :               max_subcell > max_dg + delta or min_subcell < min_dg - delta;
      98             :           if (cell_is_troubled) {
      99             :             tci_status = static_cast<int>(tag_index + 1);
     100             :             return;
     101             :           }
     102             :           ++tag_index;
     103             :         }
     104             :       });
     105             :   return tci_status;
     106             : }
     107             : }  // namespace evolution::dg::subcell