Update subcell projection to work with ZernikeB1 bases

ZernikeB1 bases, used with Cartoon, require the parity of the component to be
known as the interpolation requires this to be sufficiently accurate. When the
projection is of a DataVector, it must either be a single component with the
parity passed as an argument, or the underlying tmpl::list to types that make up
the DataVector must be passed as a template.

The GrMhd systems are updated to correctly call the projection, while all other
systems are simply updated to work with the new interface.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: michaeldmurphy1

src/Evolution/DgSubcell/Actions/ReconstructionCommunication.hpp

@@ -34,6 +34,7 @@
 #include "Evolution/DgSubcell/CombineVolumeGhostData.hpp"
 #include "Evolution/DgSubcell/GhostData.hpp"
 #include "Evolution/DgSubcell/NeighborRdmpAndVolumeData.hpp"
+#include "Evolution/DgSubcell/NeighborRdmpAndVolumeData.tpp"
 #include "Evolution/DgSubcell/Projection.hpp"
 #include "Evolution/DgSubcell/RdmpTci.hpp"
 #include "Evolution/DgSubcell/RdmpTciData.hpp"
@@ -803,7 +804,9 @@ struct ReceiveDataForReconstruction {
                 *boundary_data.ghost_cell_data, number_of_rdmp_vars,
                 directional_element_id,
                 mesh_for_ghost_data->at(directional_element_id), element,
-                subcell_mesh, ghost_zone_size, neighbor_dg_to_fd_interpolants);
+                subcell_mesh, ghost_zone_size, neighbor_dg_to_fd_interpolants,
+                tmpl::list<typename Metavariables::SubcellOptions::
+                               GhostVariables::ghost_variables_tag_list>{});
             ASSERT(neighbor_tci_decisions->contains(directional_element_id),
                    "The NeighorTciDecisions should contain the neighbor ("
                        << directional_element_id.direction() << ", "

src/Evolution/DgSubcell/Actions/TciAndRollback.hpp

@@ -27,6 +27,7 @@
 #include "Evolution/DgSubcell/ActiveGrid.hpp"
 #include "Evolution/DgSubcell/GhostData.hpp"
 #include "Evolution/DgSubcell/NeighborRdmpAndVolumeData.hpp"
+#include "Evolution/DgSubcell/NeighborRdmpAndVolumeData.tpp"
 #include "Evolution/DgSubcell/Projection.hpp"
 #include "Evolution/DgSubcell/RdmpTci.hpp"
 #include "Evolution/DgSubcell/RdmpTciData.hpp"
@@ -247,12 +248,14 @@ struct TciAndRollback {
           for (const auto& [directional_element_id, mesh_for_ghost_data] :
                meshes_for_ghost_data) {
             evolution::dg::subcell::insert_or_update_neighbor_volume_data<
-                false>(ghost_data_ptr,
-                       ghost_data_ptr->at(directional_element_id)
-                           .neighbor_ghost_data_for_reconstruction(),
-                       0, directional_element_id, mesh_for_ghost_data, element,
-                       subcell_mesh, ghost_zone_size,
-                       neighbor_dg_to_fd_interpolants);
+                false>(
+                ghost_data_ptr,
+                ghost_data_ptr->at(directional_element_id)
+                    .neighbor_ghost_data_for_reconstruction(),
+                0, directional_element_id, mesh_for_ghost_data, element,
+                subcell_mesh, ghost_zone_size, neighbor_dg_to_fd_interpolants,
+                tmpl::list<typename Metavariables::SubcellOptions::
+                               GhostVariables::ghost_variables_tag_list>{});
           }
 
           // Note: We do _not_ project the boundary history here because

src/Evolution/DgSubcell/CMakeLists.txt

@@ -27,6 +27,7 @@ spectre_target_headers(
   Matrices.hpp
   Mesh.hpp
   NeighborRdmpAndVolumeData.hpp
+  NeighborRdmpAndVolumeData.tpp
   NeighborReconstructedFaceSolution.hpp
   NeighborReconstructedFaceSolution.tpp
   NeighborTciDecision.hpp
@@ -59,7 +60,6 @@ spectre_target_sources(
   InitialTciData.cpp
   Matrices.cpp
   Mesh.cpp
-  NeighborRdmpAndVolumeData.cpp
   NeighborTciDecision.cpp
   PerssonTci.cpp
   Projection.cpp

src/Evolution/DgSubcell/CellCenteredFlux.hpp

@@ -18,6 +18,7 @@
 #include "Evolution/DgSubcell/Tags/Mesh.hpp"
 #include "Evolution/DgSubcell/Tags/SubcellOptions.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
+#include "NumericalAlgorithms/Spectral/Parity.hpp"
 #include "Utilities/Gsl.hpp"
 #include "Utilities/TMPL.hpp"
 
@@ -74,7 +75,8 @@ struct CellCenteredFlux {
             const DataVector& cell_centered_mesh_velocity =
                 evolution::dg::subcell::fd::project(
                     dg_mesh_velocity.value().get(i), dg_mesh,
-                    subcell_mesh.extents());
+                    subcell_mesh.extents(),
+                    i == 0 ? Spectral::Parity::Odd : Spectral::Parity::Even);
 
             tmpl::for_each<flux_variables>(
                 [&cell_centered_flux_vars, &cell_centered_mesh_velocity,

src/Evolution/DgSubcell/CorrectPackagedData.hpp

@@ -122,7 +122,8 @@ void correct_package_data(
           projected_data = DataVector{
               Variables<DgPackageFieldTags>::number_of_independent_components *
               subcell_face_mesh.number_of_grid_points()};
-          evolution::dg::subcell::fd::detail::project_impl(
+          evolution::dg::subcell::fd::detail::project_impl_with_tag_list<
+              DgPackageFieldTags>(
               gsl::make_span(projected_data.data(), projected_data.size()),
               gsl::make_span(
                   slice_data,

src/Evolution/DgSubcell/Matrices.cpp

@@ -19,6 +19,7 @@
 #include "NumericalAlgorithms/Spectral/MaximumNumberOfPoints.hpp"
 #include "NumericalAlgorithms/Spectral/Mesh.hpp"
 #include "NumericalAlgorithms/Spectral/MinimumNumberOfPoints.hpp"
+#include "NumericalAlgorithms/Spectral/Parity.hpp"
 #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
 #include "NumericalAlgorithms/Spectral/QuadratureWeights.hpp"
 #include "Utilities/Algorithm.hpp"
@@ -31,12 +32,15 @@
 #include "Utilities/StaticCache.hpp"
 
 namespace evolution::dg::subcell::fd {
-const Matrix& projection_matrix(
-    const Mesh<1>& dg_mesh, const size_t subcell_extents,
-    const Spectral::Quadrature& subcell_quadrature) {
+const Matrix& projection_matrix(const Mesh<1>& dg_mesh,
+                                const size_t subcell_extents,
+                                const Spectral::Quadrature& subcell_quadrature,
+                                const Spectral::Parity parity) {
   ASSERT(dg_mesh.basis(0) == Spectral::Basis::Legendre or
-             dg_mesh.basis(0) == Spectral::Basis::Cartoon,
-         "FD Subcell projection only supports Legendre or Cartoon bases right "
+             dg_mesh.basis(0) == Spectral::Basis::Cartoon or
+             dg_mesh.basis(0) == Spectral::Basis::ZernikeB1,
+         "FD Subcell projection only supports Legendre, Cartoon, or ZernikeB1 "
+         "bases right "
          "now but got basis "
              << dg_mesh.basis(0));
   ASSERT(subcell_quadrature == Spectral::Quadrature::FaceCentered or
@@ -46,7 +50,39 @@ const Matrix& projection_matrix(
          "subcell_quadrature option in projection_matrix should be "
          "FaceCentered, CellCentered, or a Cartoon quadrature, but got "
              << subcell_quadrature);
-  static const auto cache = make_static_cache<
+  ASSERT(dg_mesh.basis(0) != Spectral::Basis::ZernikeB1 or
+             parity != Spectral::Parity::Uninitialized,
+         "Parity must be set when using ZernikeB1");
+
+  const static auto zernike_b1_cache = make_static_cache<
+      CacheRange<
+          Spectral::minimum_number_of_points<
+              Spectral::Basis::ZernikeB1,
+              Spectral::Quadrature::GaussRadauUpper>,
+          Spectral::maximum_number_of_points<Spectral::Basis::ZernikeB1> + 1>,
+      CacheRange<Spectral::minimum_number_of_points<
+                     Spectral::Basis::FiniteDifference,
+                     Spectral::Quadrature::CellCentered>,
+                 Spectral::maximum_number_of_points<
+                     Spectral::Basis::FiniteDifference> +
+                     1>,
+      CacheEnumeration<Spectral::Quadrature, Spectral::Quadrature::CellCentered,
+                       Spectral::Quadrature::FaceCentered>,
+      CacheEnumeration<Spectral::Parity, Spectral::Parity::Even,
+                       Spectral::Parity::Odd>>(
+      [](const size_t local_num_dg_points, const size_t local_num_fd_points,
+         const Spectral::Quadrature local_subcell_quadrature,
+         const Spectral::Parity local_parity) {
+        return Spectral::interpolation_matrix<
+            Spectral::Basis::ZernikeB1, Spectral::Quadrature::GaussRadauUpper>(
+            local_num_dg_points,
+            Spectral::collocation_points(
+                Mesh<1>{local_num_fd_points, Spectral::Basis::FiniteDifference,
+                        local_subcell_quadrature}),
+            local_parity);
+      });
+
+  const static auto legendre_cache = make_static_cache<
       CacheRange<
           Spectral::minimum_number_of_points<
               Spectral::Basis::Legendre, Spectral::Quadrature::GaussLobatto>,
@@ -74,9 +110,13 @@ const Matrix& projection_matrix(
   if (dg_mesh.basis(0) == Spectral::Basis::Cartoon) {
     static const Matrix cartoon_matrix{1, 1, 1.0};
     return cartoon_matrix;
+  } else if (dg_mesh.basis(0) == Spectral::Basis::ZernikeB1) {
+    return zernike_b1_cache(dg_mesh.extents(0), subcell_extents,
+                            subcell_quadrature, parity);
+  } else {
+    return legendre_cache(dg_mesh.extents(0), subcell_extents,
+                          dg_mesh.quadrature(0), subcell_quadrature);
   }
-  return cache(dg_mesh.extents(0), subcell_extents, dg_mesh.quadrature(0),
-               subcell_quadrature);
 }
 
 namespace {

src/Evolution/DgSubcell/Matrices.hpp

@@ -14,6 +14,7 @@ template <size_t Dim>
 class Mesh;
 enum class Side : uint8_t;
 namespace Spectral {
+enum class Parity : uint8_t;
 enum class Quadrature : uint8_t;
 }  // namespace Spectral
 /// \endcond
@@ -23,9 +24,12 @@ namespace evolution::dg::subcell::fd {
  * \ingroup DgSubcellGroup
  * \brief Computes the projection matrix in 1 dimension going from a DG
  * mesh to a conservative finite difference subcell mesh.
+ *
+ * The parity parameter is required for ZernikeB1 bases.
  */
 const Matrix& projection_matrix(const Mesh<1>& dg_mesh, size_t subcell_extents,
-                                const Spectral::Quadrature& subcell_quadrature);
+                                const Spectral::Quadrature& subcell_quadrature,
+                                Spectral::Parity parity);
 
 /*!
  * \ingroup DgSubcellGroup
@@ -95,4 +99,4 @@ const Matrix& reconstruction_matrix(const Mesh<Dim>& dg_mesh,
  */
 const Matrix& projection_matrix(const Mesh<1>& dg_mesh, size_t subcell_extents,
                                 size_t ghost_zone_size, Side side);
-}  // namespace evolution::dg::subcellfd
+}  // namespace evolution::dg::subcell::fd

src/Evolution/DgSubcell/NeighborRdmpAndVolumeData.hpp

@@ -15,6 +15,7 @@
 #include "Evolution/DgSubcell/RdmpTciData.hpp"
 #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
 #include "Utilities/Gsl.hpp"
+#include "Utilities/TMPL.hpp"
 
 /// \cond
 template <size_t Dim>
@@ -31,25 +32,34 @@ namespace evolution::dg::subcell {
  *
  * This is intended to be used during a rollback from DG to make sure neighbor
  * data is projected to the FD grid.
+ *
+ * The `tmpl::list<VolumeFields>` meta parameter encodes the tag list of the
+ * packed ghost DataVector so that per-component parity information can be
+ * deduced for ZernikeB1 meshes.
  */
-template <bool InsertIntoMap, size_t Dim>
+template <bool InsertIntoMap, size_t Dim, typename VolumeFields>
 void insert_or_update_neighbor_volume_data(
     gsl::not_null<DirectionalIdMap<Dim, GhostData>*> ghost_data_ptr,
     const DataVector& neighbor_subcell_data,
-    const size_t number_of_rdmp_vars_in_buffer,
+    size_t number_of_rdmp_vars_in_buffer,
     const DirectionalId<Dim>& directional_element_id,
     const Mesh<Dim>& neighbor_mesh, const Element<Dim>& element,
     const Mesh<Dim>& subcell_mesh, size_t number_of_ghost_zones,
     const DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>&
-        Neighbor_dg_to_fd_interpolants);
+        neighbor_dg_to_fd_interpolants,
+    tmpl::list<VolumeFields> /*meta*/);
 
 /*!
  * \brief Check whether the neighbor sent is DG volume or FD ghost data, and
  * orient project DG volume data if necessary.
  *
  * This is intended to be used by the `ReceiveDataForReconstruction` action.
+ *
+ * The `tmpl::list<VolumeFields>` meta parameter encodes the tag list of the
+ * packed ghost DataVector so that per-component parity information can be
+ * deduced for ZernikeB1 meshes.
  */
-template <size_t Dim>
+template <size_t Dim, typename VolumeFields>
 void insert_neighbor_rdmp_and_volume_data(
     gsl::not_null<RdmpTciData*> rdmp_tci_data_ptr,
     gsl::not_null<DirectionalIdMap<Dim, GhostData>*> ghost_data_ptr,
@@ -59,5 +69,6 @@ void insert_neighbor_rdmp_and_volume_data(
     const Mesh<Dim>& neighbor_mesh, const Element<Dim>& element,
     const Mesh<Dim>& subcell_mesh, size_t number_of_ghost_zones,
     const DirectionalIdMap<Dim, std::optional<intrp::Irregular<Dim>>>&
-        neighbor_dg_to_fd_interpolants);
+        neighbor_dg_to_fd_interpolants,
+    tmpl::list<VolumeFields> /*meta*/);
 }  // namespace evolution::dg::subcell