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