Test various reconstruction alternatives [ci-skip]

This is committed to ease discussions

src/scheme/test_developed_reconstruction.cpp

@@ -914,29 +914,38 @@ test_developed_reconstruction3(const std::shared_ptr<const MeshType>& p_mesh,
 
   compute_HP(rho_bar_cons, rhou_bar, u_h, HP_rhou);
 
-  auto u_bar_R = [&](const CellId cell_id, const Rd& x) {
+  [[maybe_unused]] auto u_bar_R = [&](const CellId cell_id, const Rd& x) {
     return u_h[cell_id] + (1. / rho_bar_cons[cell_id](x)) * HP_rhou[cell_id](x);
   };
 
-  if constexpr (Dimension == 1) {
-    std::cout << "done\n";
-
   if constexpr (Dimension == 1) {
     auto xr                  = p_mesh->xr();
     auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
 
+    auto rho_exact = [](const Rd& x) { return 2 + 3 * x[0] + 2 * (x[0] * x[0]) + 0.2 * (x[0] * x[0] * x[0]); };
+    auto u_exect   = [](const Rd& x) -> Rd { return Rd{-1 + 2 * x[0] + 3 * (x[0] * x[0]) - 3 * (x[0] * x[0] * x[0])}; };
+
+    auto rho_u_exact = [&](const Rd& x) { return rho_exact(x) * u_exect(x); };
+
     {
-        double sum = 0;
+      double L1    = 0;
+      double maxL1 = 0;
+      double Linf  = 0;
       auto qf      = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{10});
       for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
         LineTransformation<Dimension> t{xr[cell_to_node_matrix[cell_id][0]], xr[cell_to_node_matrix[cell_id][1]]};
+        double localL1 = 0;
         for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
           auto x = t(qf.point(i_q));
-            auto diff = (rho_u_bar_cons[cell_id](x) - rhou_bar[cell_id](x));
-            sum += std::abs(qf.weight(i_q) * std::abs(diff[0]));
+          // auto diff = (rho_u_bar_cons[cell_id](x) - rhou_bar[cell_id](x));
+          auto diff = (rho_u_exact(x) - rhou_bar[cell_id](x));
+          localL1 += std::abs(qf.weight(i_q) * std::abs(diff[0]));
+          Linf = std::max(Linf, std::abs(diff[0]));
         }
+        L1 += localL1;
+        maxL1 = std::max(maxL1, localL1);
       }
-        std::cout << "rho u: L1 error = " << sum << '\n';
+      std::cout << "errors: [L1, maxL1, Linf] = [" << L1 << ' ' << maxL1 << ' ' << Linf << "]\n";
     }
 
     // {
@@ -954,6 +963,145 @@ test_developed_reconstruction3(const std::shared_ptr<const MeshType>& p_mesh,
     // }
   }
 }
+
+template <MeshConcept MeshType>
+void
+test_developed_reconstruction4(const std::shared_ptr<const MeshType>& p_mesh,
+                               const std::shared_ptr<const DiscreteFunctionVariant>& rho_v,
+                               const std::shared_ptr<const DiscreteFunctionVariant>& u_v,
+                               const std::shared_ptr<const DiscreteFunctionVariant>& rhou_v)
+{
+  static constexpr size_t Dimension = MeshType::Dimension;
+  using Rd                          = TinyVector<Dimension>;
+  const size_t degree               = 3;
+
+  PolynomialReconstructionDescriptor descriptor{IntegrationMethodType::element, degree};
+
+  auto rho_h  = rho_v->get<DiscreteFunctionP0<const double>>();
+  auto u_h    = u_v->get<DiscreteFunctionP0<const double>>();
+  auto rhou_h = rhou_v->get<DiscreteFunctionP0<const double>>();
+
+  auto reconstructions = PolynomialReconstruction{descriptor}.build(rho_v, u_v);
+
+  PolynomialReconstructionDescriptor descriptor_deg_6{IntegrationMethodType::boundary, 6};
+
+  auto reconstructions_deg_6 = PolynomialReconstruction{descriptor_deg_6}.build(rhou_h);
+
+  auto rho_u_bar_cons = reconstructions_deg_6[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+  auto polynomial_basis_means = PolynomialBasisDataManager::instance()
+                                  .getPolynomialBasisData(*p_mesh)
+                                  .getPolynomialBasisMeans(descriptor.integrationMethodType(), descriptor.degree());
+
+  auto remove_mean = [&](const auto& f_h, const auto& g_h, auto fg_bar) {
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+        auto basis_means   = polynomial_basis_means[cell_id];
+        using DataType_fg  = std::decay_t<decltype(fg_bar_coeffs[0])>;
+
+        DataType_fg fg_bar_mean;
+        if constexpr (std::is_arithmetic_v<DataType_fg>) {
+          fg_bar_mean = 0;
+        } else {
+          fg_bar_mean = zero;
+        }
+
+        for (size_t i = 0; i < basis_means.size(); ++i) {
+          fg_bar_mean += basis_means[i] * fg_bar_coeffs[i];
+        }
+
+        fg_bar_coeffs[0] += f_h[cell_id] * g_h[cell_id] - fg_bar_mean;
+      });
+  };
+
+  auto compute_HP = [&](const auto& f_bar, const auto& fg_bar, const auto& g_h, auto& HP_fg) {
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        auto fg_bar_coeffs = fg_bar.coefficients(cell_id);
+        auto f_bar_coeffs  = f_bar.coefficients(cell_id);
+
+        auto HP_fg_coeffs = HP_fg.coefficients(cell_id);
+
+        const auto& gj = g_h[cell_id];
+
+        for (size_t i = 0; i < HP_fg_coeffs.size(); ++i) {
+          HP_fg_coeffs[i] = fg_bar_coeffs[i] - f_bar_coeffs[i] * gj;
+        }
+      });
+  };
+
+  auto rho_bar_cons = reconstructions[0]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+  auto u_bar_cons   = reconstructions[1]->template get<DiscreteFunctionDPk<Dimension, const double>>();
+
+  DiscreteFunctionDPk<Dimension, const double> rho_bar = [&] {
+    auto v = copy(rho_bar_cons);
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        auto v_coeffs = v.coefficients(cell_id);
+        v_coeffs[0]   = rho_h[cell_id];
+      });
+    return v;
+  }();
+
+  DiscreteFunctionDPk<Dimension, const double> u_bar = [&] {
+    auto v = copy(u_bar_cons);
+    parallel_for(
+      p_mesh->numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+        auto u_bar_coeffs = v.coefficients(cell_id);
+        u_bar_coeffs[0]   = u_h[cell_id];
+      });
+    return v;
+  }();
+
+  DiscreteFunctionDPk<Dimension, double> rhou_bar{p_mesh, rho_bar.degree()};
+
+  computeFGBar(p_mesh, rho_bar, u_bar, rhou_bar);
+  remove_mean(rho_h, u_h, rhou_bar);
+
+  DiscreteFunctionDPk<Dimension, double> HP_rhou{p_mesh, rho_bar.degree()};
+
+  compute_HP(rho_bar_cons, rhou_bar, u_h, HP_rhou);
+
+  [[maybe_unused]] auto u_bar_R = [&](const CellId cell_id, const Rd& x) {
+    return u_h[cell_id] + (1. / rho_bar_cons[cell_id](x)) * HP_rhou[cell_id](x);
+  };
+
+  if constexpr (Dimension == 1) {
+    auto xr                  = p_mesh->xr();
+    auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
+
+    auto rho_exact = [](const Rd& x) { return 2 + 3 * x[0] + 2 * (x[0] * x[0]) + 0.2 * (x[0] * x[0] * x[0]); };
+    auto u_exact   = [](const Rd& x) { return -1 + 2 * x[0] + 3 * (x[0] * x[0]) - 3 * (x[0] * x[0] * x[0]); };
+
+    auto rho_u_exact = [&](const Rd& x) { return rho_exact(x) * u_exact(x); };
+
+    {
+      double L1_deg_6 = 0;
+
+      double L1    = 0;
+      double maxL1 = 0;
+      double Linf  = 0;
+      auto qf      = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor{10});
+      for (CellId cell_id = 0; cell_id < p_mesh->numberOfCells(); ++cell_id) {
+        LineTransformation<Dimension> t{xr[cell_to_node_matrix[cell_id][0]], xr[cell_to_node_matrix[cell_id][1]]};
+        double localL1 = 0;
+        for (size_t i_q = 0; i_q < qf.numberOfPoints(); ++i_q) {
+          auto x = t(qf.point(i_q));
+          // auto diff = (rho_u_bar_cons[cell_id](x) - rhou_bar[cell_id](x));
+          auto diff = (rho_u_exact(x) - rhou_bar[cell_id](x));
+          localL1 += std::abs(qf.weight(i_q) * std::abs(diff));
+          Linf = std::max(Linf, std::abs(diff));
+
+          L1_deg_6 += std::abs(qf.weight(i_q) * (rho_u_exact(x) - rho_u_bar_cons[cell_id](x)));
+        }
+        L1 += localL1;
+        maxL1 = std::max(maxL1, localL1);
+      }
+      std::cout << "errors: [L1, maxL1, Linf] = [" << L1 << ' ' << maxL1 << ' ' << Linf << "]\n";
+      std::cout << "errors: [L1_deg_6] = [" << L1_deg_6 << "]\n";
+    }
+  }
 }
 
 void
@@ -963,5 +1111,5 @@ test_developed_reconstruction(const std::shared_ptr<const DiscreteFunctionVarian
 {
   auto mesh_v = getCommonMesh({rho_v, u_v, E_v});
 
-  std::visit([&](auto&& p_mesh) { test_developed_reconstruction3(p_mesh, rho_v, u_v, E_v); }, mesh_v->variant());
+  std::visit([&](auto&& p_mesh) { test_developed_reconstruction4(p_mesh, rho_v, u_v, E_v); }, mesh_v->variant());
 }