diff --git a/src/scheme/test_developed_reconstruction.cpp b/src/scheme/test_developed_reconstruction.cpp
index 6a4e4e6a030724dc671e90e534acb61c4d948b71..291885bcd6771936a218fc954a7fb7ad1e883615 100644
--- a/src/scheme/test_developed_reconstruction.cpp
+++ b/src/scheme/test_developed_reconstruction.cpp
@@ -914,44 +914,192 @@ 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";
+ auto xr = p_mesh->xr();
+ auto cell_to_node_matrix = p_mesh->connectivity().cellToNodeMatrix();
- 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])}; };
- {
- double sum = 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]]};
- 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 rho_u_exact = [&](const Rd& x) { return rho_exact(x) * u_exect(x); };
+
+ {
+ 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[0]));
+ Linf = std::max(Linf, std::abs(diff[0]));
}
- std::cout << "rho u: L1 error = " << sum << '\n';
+ L1 += localL1;
+ maxL1 = std::max(maxL1, localL1);
}
+ std::cout << "errors: [L1, maxL1, Linf] = [" << L1 << ' ' << maxL1 << ' ' << Linf << "]\n";
+ }
+
+ // {
+ // double sum = 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]]};
+ // 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) - rho_bar_cons[cell_id](x) * (u_bar_R(cell_id, x)));
+ // sum += std::abs(qf.weight(i_q) * std::abs(diff[0]));
+ // }
+ // }
+ // std::cout << "rho u: L1 error = " << sum << '\n';
+ // }
+ }
+}
+
+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;
+ }
+ });
+ };
- // {
- // double sum = 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]]};
- // 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) - rho_bar_cons[cell_id](x) * (u_bar_R(cell_id, x)));
- // sum += std::abs(qf.weight(i_q) * std::abs(diff[0]));
- // }
- // }
- // std::cout << "rho u: L1 error = " << sum << '\n';
- // }
+ 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";
}
}
}
@@ -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());
}