Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : #pragma once
       4             : 
       5             : #include <cstddef>
       6             : #include <cstdint>
       7             : 
       8             : /// \cond
       9             : class DataVector;
      10             : template <size_t Dim>
      11             : class Index;
      12             : class Matrix;
      13             : template <size_t Dim>
      14             : class Mesh;
      15             : enum class Side : uint8_t;
      16             : namespace Spectral {
      17             : enum class Parity : uint8_t;
      18             : enum class Quadrature : uint8_t;
      19             : }  // namespace Spectral
      20             : /// \endcond
      21             : 
      22             : namespace evolution::dg::subcell::fd {
      23             : /*!
      24             :  * \ingroup DgSubcellGroup
      25             :  * \brief Computes the projection matrix in 1 dimension going from a DG
      26             :  * mesh to a conservative finite difference subcell mesh.
      27             :  *
      28             :  * The parity parameter is required for ZernikeB1 bases, and is ignored for
      29             :  * others.
      30             :  */
      31           1 : const Matrix& projection_matrix(const Mesh<1>& dg_mesh, size_t subcell_extents,
      32             :                                 const Spectral::Quadrature& subcell_quadrature,
      33             :                                 Spectral::Parity parity);
      34             : 
      35             : /*!
      36             :  * \ingroup DgSubcellGroup
      37             :  * \brief Computes the matrix needed for reconstructing the DG solution from
      38             :  * the subcell solution.
      39             :  *
      40             :  * Reconstructing the DG solution from the FD solution is a bit more
      41             :  * involved than projecting the DG solution to the FD subcells. Denoting the
      42             :  * projection operator by \f$\mathcal{P}\f$ and the reconstruction operator by
      43             :  * \f$\mathcal{R}\f$, we desire the property
      44             :  *
      45             :  * \f{align*}{
      46             :  *   \mathcal{R}(\mathcal{P}(u_{\breve{\imath}}
      47             :  *   J_{\breve{\imath}}))=u_{\breve{\imath}} J_{\breve{\imath}},
      48             :  * \f}
      49             :  *
      50             :  * where \f$\breve{\imath}\f$ denotes a grid point on the DG grid, \f$u\f$ is
      51             :  * the solution on the DG grid, and \f$J\f$ is the determinant of the Jacobian
      52             :  * on the DG grid. We also require that the integral of the conserved variables
      53             :  * over the subcells is equal to the integral over the DG element. That is,
      54             :  *
      55             :  * \f{align*}{
      56             :  *   \int_{\Omega}u \,d^3x =\int_{\Omega} \underline{u} \,d^3x \Longrightarrow
      57             :  *   \int_{\Omega}u J \,d^3\xi=\int_{\Omega} \underline{u} J \,d^3\xi,
      58             :  * \f}
      59             :  *
      60             :  * where \f$\underline{u}\f$ is the solution on the subcells. Because the number
      61             :  * of subcell points is larger than the number of DG points, we need to solve a
      62             :  * constrained linear least squares problem to reconstruct the DG solution from
      63             :  * the subcells.
      64             :  *
      65             :  * The final reconstruction matrix is given by
      66             :  *
      67             :  * \f{align*}{
      68             :  *   R_{\breve{\jmath}\underline{i}}
      69             :  *   &=\left\{(2 \mathcal{P}\otimes\mathcal{P})^{-1}2\mathcal{P} - (2
      70             :  *   \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\left[\mathbf{w}(2
      71             :  *   \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\right]^{-1}\mathbf{w}(2
      72             :  *   \mathcal{P}\otimes\mathcal{P})^{-1}2\mathcal{P}
      73             :  *   + (2 \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\left[\mathbf{w}(2
      74             :  *   \mathcal{P}\otimes\mathcal{P})^{-1}\vec{w}\right]^{-1}\vec{\underline{w}}
      75             :  *   \right\}_{\breve{\jmath}\underline{i}},
      76             :  * \f}
      77             :  *
      78             :  * where \f$\vec{w}\f$ is the vector of integration weights on the DG element,
      79             :  * \f$\mathbf{w}=w_{\breve{l}}\delta_{\breve{l}\breve{\jmath}}\f$, and
      80             :  * \f$\vec{\underline{w}}\f$ is the vector of integration weights over the
      81             :  * subcells. The integration weights \f$\vec{\underline{w}}\f$ on the subcells
      82             :  * are those for 6th-order integration on a uniform mesh.
      83             :  */
      84             : template <size_t Dim>
      85           1 : const Matrix& reconstruction_matrix(const Mesh<Dim>& dg_mesh,
      86             :                                     const Index<Dim>& subcell_extents);
      87             : 
      88             : /*!
      89             :  * \ingroup DgSubcellGroup
      90             :  * \brief Computes the 1D reconstruction matrix for a ZernikeB1 DG mesh.
      91             :  *
      92             :  * Uses the same constrained least-squares formula as the Legendre overload,
      93             :  * but with the parity-aware ZernikeB1 interpolation matrix and GaussRadauUpper
      94             :  * quadrature weights, notably only for even parity. For odd parity, the
      95             :  * collocation points are inconsistent with the constraint equation, so we opt
      96             :  * to directly use the pseduo-inverse. Only `DimByDim` reconstruction is
      97             :  * supported for ZernikeB1  meshes; this overload is called per-dimension for
      98             :  * the radial direction.
      99             :  *
     100             :  * The `parity` parameter must be `Even` or `Odd` (not `Uninitialized`).
     101             :  */
     102           1 : const Matrix& reconstruction_matrix(const Mesh<1>& dg_mesh,
     103             :                                     size_t subcell_extents,
     104             :                                     Spectral::Parity parity);
     105             : 
     106             : /*!
     107             :  * \ingroup DgSubcellGroup
     108             :  * \brief Computes the projection matrix in 1 dimension going from a DG
     109             :  * mesh to a conservative finite difference subcell mesh for only the ghost
     110             :  * zones.
     111             :  *
     112             :  * This is used when a neighbor sends DG volume data and we need to switch to
     113             :  * FD. In this case we need to project the DG volume data onto the ghost zone
     114             :  * cells.
     115             :  *
     116             :  * \note Currently assumes a max ghost zone size of `5` and a minimum ghost zone
     117             :  * size of 2.
     118             :  */
     119           1 : const Matrix& projection_matrix(const Mesh<1>& dg_mesh, size_t subcell_extents,
     120             :                                 size_t ghost_zone_size, Side side);
     121             : }  // namespace evolution::dg::subcell::fd