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 : * The parity parameter is required for ZernikeB1 bases, and is ignored for
117 : * others.
118 : *
119 : * \note Currently assumes a max ghost zone size of `5` and a minimum ghost zone
120 : * size of 2.
121 : */
122 1 : const Matrix& projection_matrix(const Mesh<1>& dg_mesh, size_t subcell_extents,
123 : size_t ghost_zone_size, Side side,
124 : Spectral::Parity parity);
125 : } // namespace evolution::dg::subcell::fd
|
