Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : 
       9             : #include "NumericalAlgorithms/Spectral/Mesh.hpp"
      10             : #include "Utilities/Gsl.hpp"
      11             : #include "Utilities/TMPL.hpp"
      12             : 
      13             : /// \cond
      14             : class DataVector;
      15             : /// \endcond
      16             : 
      17             : /*!
      18             :  * \brief Items for assessing truncation error in spectral methods.
      19             :  */
      20           1 : namespace PowerMonitors {
      21             : 
      22             : /// @{
      23             : /*!
      24             :  * \ingroup SpectralGroup
      25             :  * \brief Returns array of power monitors in each spatial dimension.
      26             :  *
      27             :  * Computed following Sec. 5.1 of Ref. \cite Szilagyi2014fna.
      28             :  * For example, in the x dimension (indexed by \f$ k_0 \f$), we compute
      29             :  *
      30             :  * \f{align*}{
      31             :  *  P_{k_0}[\psi] = \sqrt{ \frac{1}{N_1 N_2}
      32             :  *   \sum_{k_1,k_2} \left| C_{k_0,k_1,k_2} \right|^2} ,
      33             :  * \f}
      34             :  *
      35             :  * where \f$ C_{k_0,k_1,k_2}\f$ are the modal coefficients
      36             :  * of variable \f$ \psi \f$.
      37             :  *
      38             :  */
      39             : template <size_t Dim>
      40           1 : void power_monitors(gsl::not_null<std::array<DataVector, Dim>*> result,
      41             :                     const DataVector& u, const Mesh<Dim>& mesh);
      42             : 
      43             : template <size_t Dim>
      44           1 : std::array<DataVector, Dim> power_monitors(const DataVector& u,
      45             :                                            const Mesh<Dim>& mesh);
      46             : /// @}
      47             : 
      48             : /// @{
      49             : /*!
      50             :  * \ingroup SpectralGroup
      51             :  * \brief Compute the negative log10 of the relative truncation error.
      52             :  *
      53             :  * Truncation error according to Eqs. (57) and (58) of Ref.
      54             :  * \cite Szilagyi2014fna, i.e.,
      55             :  *
      56             :  * \f{align*}{
      57             :  *  \mathcal{T}\left[P_k\right] = \log_{10} \max \left(P_0, P_1\right)
      58             :  *   - \dfrac{\sum_{j=0}^{j_{\text{max}, k}} \log_{10} \left(P_j\right) w_j}
      59             :  *   {\sum_{j=0}^{j_{\text{max}, k}} w_j} , \f}
      60             :  *
      61             :  * with weights
      62             :  *
      63             :  * \f{align*}{
      64             :  *  w_j = \exp\left[ - \left(j - j_{\text{max}, k}
      65             :  *          + \dfrac{1}{2}\right)^2 \right] .
      66             :  * \f}
      67             :  *
      68             :  * where \f$ j_{\text{max}, k}  = N_k - 1 \f$ and  \f$ N_k \f$ is the number of
      69             :  * modes or gridpoints in dimension k. Here the second term is a weighted
      70             :  * average with larger weights toward the highest modes. This number should
      71             :  * correspond to the number of digits resolved by the spectral expansion.
      72             :  *
      73             :  * \note Modes below a cutoff of $100 \epsilon \mathrm{max}_k(P_k)$ are ignored
      74             :  * in the weighted average, where $\epsilon$ is the machine epsilon. This
      75             :  * ensures that we don't underestimate the truncation error if some modes are
      76             :  * zero (e.g. by symmetry). Furthermore, if the last two or more modes are zero,
      77             :  * we assume that the function is represented exactly and return a relative
      78             :  * truncation error of zero.
      79             :  *
      80             :  * \details The number of modes (`num_modes_to_use`) argument needs to be less
      81             :  * or equal than the total number of power monitors (`power_monitor.size()`).
      82             :  * In contrast with Ref. \cite Szilagyi2014fna, here we index the modes starting
      83             :  * from zero.
      84             :  *
      85             :  */
      86           1 : double relative_truncation_error(const DataVector& power_monitor,
      87             :                                  const size_t num_modes_to_use);
      88             : /// @}
      89             : 
      90             : /*!
      91             :  * \brief The relative truncation error in each logical direction of the grid
      92             :  *
      93             :  * This overload is intended for visualization purposes only. It takes a tensor
      94             :  * component as input, so it can be used as a kernel to post-process volume data
      95             :  * with Python bindings (see `TransformVolumeData.py`).
      96             :  *
      97             :  * This function returns the relative truncation error directly, as opposed to
      98             :  * the other overload that returns the negative log10 of the relative truncation
      99             :  * error.
     100             :  */
     101             : template <size_t Dim>
     102           1 : std::array<double, Dim> relative_truncation_error(
     103             :     const DataVector& tensor_component, const Mesh<Dim>& mesh);
     104             : 
     105             : /// @{
     106             : /*!
     107             :  * \ingroup SpectralGroup
     108             :  * \brief Returns an estimate of the absolute truncation error in each
     109             :  * dimension.
     110             :  *
     111             :  * The estimate of the numerical error is given by
     112             :  *
     113             :  * \f{align*}{
     114             :  *  \mathcal{E}\left[P_k\right] = u_\mathrm{max} \times 10^{- \mathcal{T}[P_k]},
     115             :  * \f}
     116             :  *
     117             :  * where \f$ u_\mathrm{max} = \mathrm{max} |u|\f$ in the corresponding element
     118             :  * and \f$ \mathcal{T}[P_k] \f$ is the relative error estimate
     119             :  * computed from the power monitors \f$ P_k \f$.
     120             :  *
     121             :  * \warning This estimate is intended for visualization purposes only.
     122             :  */
     123             : template <size_t Dim>
     124           1 : std::array<double, Dim> absolute_truncation_error(
     125             :     const DataVector& tensor_component, const Mesh<Dim>& mesh);
     126             : /// @}
     127             : }  // namespace PowerMonitors