Skip to content
Snippets Groups Projects
Commit 647f38ba authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Improve tests in view of testing reconstructions of degree 2+ (WIP)

parent 0c8cd0a1
No related branches found
No related tags found
1 merge request!205High-order polynomial reconstruction
Show whitespace changes
Inline Side-by-side
tests/test_PolynomialReconstruction_degree_1.cpp
+ 137
− 531
View file @ 647f38ba
...@@ -12,7 +12,10 @@ ...@@ -12,7 +12,10 @@
#include <analysis/GaussQuadratureDescriptor.hpp> #include <analysis/GaussQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp> #include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp> #include <analysis/QuadratureManager.hpp>
#include <geometry/CubeTransformation.hpp>
#include <geometry/LineTransformation.hpp> #include <geometry/LineTransformation.hpp>
#include <geometry/PrismTransformation.hpp>
#include <geometry/PyramidTransformation.hpp>
#include <geometry/SquareTransformation.hpp> #include <geometry/SquareTransformation.hpp>
#include <geometry/TetrahedronTransformation.hpp> #include <geometry/TetrahedronTransformation.hpp>
#include <geometry/TriangleTransformation.hpp> #include <geometry/TriangleTransformation.hpp>
...@@ -104,16 +107,36 @@ max_reconstruction_error(const MeshType& mesh, ...@@ -104,16 +107,36 @@ max_reconstruction_error(const MeshType& mesh,
} }
break; break;
} }
// case CellType::Pyramid: { case CellType::Pyramid: {
// SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
// xr[cell_nodes[3]]}; xr[cell_nodes[4]]};
// auto qf = QuadratureManager::instance().getSquareFormula(GaussLobattoQuadratureDescriptor{dpk_f.degree() + 1}); auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
// for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) { for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
// auto x = T(qf.point(i_quadrarture)); auto x = T(qf.point(i_quadrarture));
// max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x))); max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
// } }
// break; break;
// } }
case CellType::Prism: {
PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
auto x = T(qf.point(i_quadrarture));
max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
}
break;
}
case CellType::Hexahedron: {
CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1});
for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
auto x = T(qf.point(i_quadrarture));
max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x)));
}
break;
}
default: { default: {
throw UnexpectedError("unexpected cell type"); throw UnexpectedError("unexpected cell type");
} }
...@@ -194,16 +217,45 @@ max_reconstruction_error( ...@@ -194,16 +217,45 @@ max_reconstruction_error(
} }
break; break;
} }
// case CellType::Pyramid: { case CellType::Pyramid: {
// SquareTransformation<MeshType::Dimension> T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], PyramidTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
// xr[cell_nodes[3]]}; xr[cell_nodes[4]]};
// auto qf = QuadratureManager::instance().getSquareFormula(GaussQuadratureDescriptor{dpk_f.degree() + 1}); auto qf = QuadratureManager::instance().getPyramidFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
// for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) { for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
// auto x = T(qf.point(i_quadrarture)); auto x = T(qf.point(i_quadrarture));
// max_error = std::max(max_error, get_max_error(dpk_f[cell_id](x), exact(x))); for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
// } max_error =
// break; std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
// } }
}
break;
}
case CellType::Prism: {
PrismTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]],
xr[cell_nodes[3]], xr[cell_nodes[4]], xr[cell_nodes[5]]};
auto qf = QuadratureManager::instance().getPrismFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
auto x = T(qf.point(i_quadrarture));
for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
max_error =
std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
}
}
break;
}
case CellType::Hexahedron: {
CubeTransformation T{xr[cell_nodes[0]], xr[cell_nodes[1]], xr[cell_nodes[2]], xr[cell_nodes[3]],
xr[cell_nodes[4]], xr[cell_nodes[5]], xr[cell_nodes[6]], xr[cell_nodes[7]]};
auto qf = QuadratureManager::instance().getCubeFormula(GaussQuadratureDescriptor{dpk_v.degree() + 1});
for (size_t i_quadrarture = 0; i_quadrarture < qf.numberOfPoints(); ++i_quadrarture) {
auto x = T(qf.point(i_quadrarture));
for (size_t i_component = 0; i_component < NbComponents; ++i_component) {
max_error =
std::max(max_error, get_max_error(dpk_v(cell_id, i_component)(x), vector_exact[i_component](x)));
}
}
break;
}
default: { default: {
throw UnexpectedError("unexpected cell type"); throw UnexpectedError("unexpected cell type");
} }
...@@ -366,12 +418,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -366,12 +418,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
std::array<std::function<R3(const R1&)>, 2> vector_exact = std::array<std::function<R3(const R1&)>, 2> vector_exact =
{[](const R1& x) -> R3 { {[](const R1& x) -> R3 { return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]}; },
return R3{+2.3 + 1.7 * x[0], -1.7 + 2.1 * x[0], +1.4 - 0.6 * x[0]}; [](const R1& x) -> R3 { return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]}; }};
},
[](const R1& x) -> R3 {
return R3{+1.6 + 0.7 * x[0], -2.1 + 1.2 * x[0], +1.1 - 0.3 * x[0]};
}};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
...@@ -619,62 +667,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -619,62 +667,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<3>>(); auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto R_affine = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; }; auto R_exact = [](const R3& x) { return 2.3 + 1.7 * x[0] - 1.3 * x[1] + 2.1 * x[2]; };
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0<double> fh{p_mesh}; DiscreteFunctionP0<double> fh{p_mesh};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); }); mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>(); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R_exact));
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error =
std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id])));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
(dpk_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0})) /
0.2;
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7));
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
(dpk_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0})) /
0.2;
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3)));
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
(dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1})) /
0.2;
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1));
}
REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
} }
} }
} }
...@@ -687,7 +693,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -687,7 +693,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<3>>(); auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto R3_affine = [](const R3& x) -> R3 { auto R3_exact = [](const R3& x) -> R3 {
return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2], // return R3{+2.3 + 1.7 * x[0] - 2.2 * x[1] + 1.8 * x[2], //
+1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2], // +1.4 - 0.6 * x[0] + 1.3 * x[1] - 3.7 * x[2], //
-0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]}; -0.2 + 3.1 * x[0] - 1.1 * x[1] + 1.9 * x[2]};
...@@ -697,53 +703,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -697,53 +703,14 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
DiscreteFunctionP0<R3> uh{p_mesh}; DiscreteFunctionP0<R3> uh{p_mesh};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); }); mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(uh); auto reconstructions = PolynomialReconstruction{descriptor}.build(uh);
auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>(); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error =
std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
dpk_uh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1}));
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
dpk_uh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope = (1 / 0.2) * (dpk_uh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
dpk_uh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9}));
}
REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
} }
} }
} }
...@@ -758,7 +725,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -758,7 +725,7 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<3>>(); auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto R2x2_affine = [](const R3& x) -> R2x2 { auto R2x2_exact = [](const R3& x) -> R2x2 {
return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2], return R2x2{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2],
// //
+2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]}; +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]};
...@@ -768,88 +735,40 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -768,88 +735,40 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
DiscreteFunctionP0<R2x2> Ah{p_mesh}; DiscreteFunctionP0<R2x2> Ah{p_mesh};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); }); mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_exact(xj[cell_id]); });
descriptor.setRowWeighting(false); descriptor.setRowWeighting(false);
auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah); auto reconstructions = PolynomialReconstruction{descriptor}.build(Ah);
auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>(); auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R2x2_exact));
{ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
double max_mean_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error =
std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id])));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0.1, 0, 0} + xj[cell_id]) -
dpk_Ah[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
max_x_slope_error =
std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1, //
-2.3, +3.1}));
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0.1, 0} + xj[cell_id]) -
dpk_Ah[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
max_y_slope_error =
std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, //
1.3, +0.8}));
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R2x2 reconstructed_slope = (1 / 0.2) * (dpk_Ah[cell_id](R3{0, 0, 0.1} + xj[cell_id]) -
dpk_Ah[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
max_z_slope_error =
std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4, //
+1.4, -1.8}));
}
REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
} }
} }
} }
SECTION("vector data") SECTION("vector data")
{ {
using R4 = TinyVector<4>;
for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) { for (auto named_mesh : MeshDataBaseForTests::get().all3DMeshes()) {
SECTION(named_mesh.name()) SECTION(named_mesh.name())
{ {
auto p_mesh = named_mesh.mesh()->get<Mesh<3>>(); auto p_mesh = named_mesh.mesh()->get<Mesh<3>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto vector_affine = [](const R3& x) -> R4 { std::array<std::function<double(const R3&)>, 4> vector_exact =
return R4{+2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2], -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2], {[](const R3& x) -> double { return +2.3 + 1.7 * x[0] + 1.2 * x[1] - 1.3 * x[2]; },
// [](const R3& x) -> double { return -1.7 + 2.1 * x[0] - 2.2 * x[1] - 2.4 * x[2]; },
+2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2], -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]}; [](const R3& x) -> double { return +2.4 - 2.3 * x[0] + 1.3 * x[1] + 1.4 * x[2]; },
}; [](const R3& x) -> double { return -0.2 + 3.1 * x[0] + 0.8 * x[1] - 1.8 * x[2]; }};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0Vector<double> Vh{p_mesh, 4}; DiscreteFunctionP0Vector<double> Vh{p_mesh, 4};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
auto vector = vector_affine(xj[cell_id]); for (size_t i = 0; i < vector_exact.size(); ++i) {
for (size_t i = 0; i < vector.dimension(); ++i) { Vh[cell_id][i] = vector_exact[i](xj[cell_id]);
Vh[cell_id][i] = vector[i];
} }
}); });
...@@ -858,58 +777,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -858,58 +777,8 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>(); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
auto vector = vector_affine(xj[cell_id]);
for (size_t i = 0; i < vector.dimension(); ++i) {
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i]));
}
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
const R4 slope{+1.7, 2.1, -2.3, +3.1};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
for (size_t i = 0; i < slope.dimension(); ++i) {
const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0.1, 0, 0} + xj[cell_id]) -
dpk_Vh(cell_id, i)(xj[cell_id] - R3{0.1, 0, 0}));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i]));
}
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_y_slope_error = 0;
const R4 slope{+1.2, -2.2, 1.3, +0.8};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
for (size_t i = 0; i < slope.dimension(); ++i) {
const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0.1, 0} + xj[cell_id]) -
dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0.1, 0}));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i]));
}
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
{
double max_z_slope_error = 0;
const R4 slope{-1.3, -2.4, +1.4, -1.8};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
for (size_t i = 0; i < slope.dimension(); ++i) {
const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, i)(R3{0, 0, 0.1} + xj[cell_id]) -
dpk_Vh(cell_id, i)(xj[cell_id] - R3{0, 0, 0.1}));
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i]));
}
}
REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
}
} }
} }
} }
...@@ -1070,37 +939,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -1070,37 +939,20 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto R1_affine = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; }; auto R1_exact = [](const R1& x) { return R1{1.7 * (x[0] + 1)}; };
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0<R1> fh{p_mesh}; DiscreteFunctionP0<R1> fh{p_mesh};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R1_affine(xj[cell_id]); }); mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R1_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1>>(); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R1>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R1_exact));
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error =
std::max(max_mean_error, std::abs((dpk_fh[cell_id](xj[cell_id]) - R1_affine(xj[cell_id]))[0]));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
(dpk_fh[cell_id](R1{0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R1{0.1}))[0] / 0.2;
max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
}
REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
}
} }
SECTION("R1 vector data") SECTION("R1 vector data")
...@@ -1111,9 +963,9 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -1111,9 +963,9 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>(); auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto vector_affine0 = [](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; }; std::array<std::function<R1(const R1&)>, 2> vector_exact //
= {[](const R1& x) -> R1 { return R1{+1.7 * (x[0] + 1)}; },
auto vector_affine1 = [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }; [](const R1& x) -> R1 { return R1{-0.3 * (x[0] + 1)}; }};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
...@@ -1121,54 +973,16 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -1121,54 +973,16 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
Vh[cell_id][0] = vector_affine0(xj[cell_id]); Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
Vh[cell_id][1] = vector_affine1(xj[cell_id]); Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
}); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1>>(); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R1>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error =
std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 0)(xj[cell_id]) - vector_affine0(xj[cell_id])));
max_mean_error =
std::max(max_mean_error, l2Norm(dpk_Vh(cell_id, 1)(xj[cell_id]) - vector_affine1(xj[cell_id])));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_slope_error = 0;
{
const TinyVector<1> slope0{+1.7};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
for (size_t i = 0; i < R1::Dimension; ++i) {
const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R1{0.1} + xj[cell_id])[i] -
dpk_Vh(cell_id, 0)(xj[cell_id] - R1{0.1})[i]);
max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope0[i]));
}
}
}
{
const TinyVector<1> slope1{-0.3};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
for (size_t i = 0; i < R1::Dimension; ++i) {
const double reconstructed_slope = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R1{0.1} + xj[cell_id])[i] -
dpk_Vh(cell_id, 1)(xj[cell_id] - R1{0.1})[i]);
max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope1[i]));
}
}
}
REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14));
}
} }
} }
} }
...@@ -1211,128 +1025,47 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -1211,128 +1025,47 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto R2_affine = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; }; auto R2_exact = [](const R2& x) -> R2 { return R2{2.3 * (x[0] - 2), -1.3 * (x[1] - 1)}; };
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0<R2> fh{p_mesh}; DiscreteFunctionP0<R2> fh{p_mesh};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R2_affine(xj[cell_id]); }); mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R2_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>(); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R2_exact));
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error =
std::max(max_mean_error, l2Norm(dpk_fh[cell_id](xj[cell_id]) - R2_affine(xj[cell_id])));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R2 reconstructed_slope =
1. / 0.2 * (dpk_fh[cell_id](R2{0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0.1, 0}));
max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R2{2.3, 0}));
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R2 reconstructed_slope =
1 / 0.2 * (dpk_fh[cell_id](R2{0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R2{0, 0.1}));
max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - (R2{0, -1.3})));
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
}
} }
SECTION("vector of R2") SECTION("vector of R2")
{ {
using R4 = TinyVector<4>;
auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>(); auto p_mesh = MeshDataBaseForTests::get().hybrid2DMesh()->get<Mesh<2>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto vector_affine = [](const R2& x) -> R4 { std::array<std::function<R2(const R2&)>, 2> vector_exact //
return R4{+1.7 * (x[0] - 2), // = {[](const R2& x) -> R2 { return R2{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1)}; },
-0.6 * (x[1] - 1), // [](const R2& x) -> R2 { return R2{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1)}; }};
-2.3 * (x[0] - 2), //
+1.1 * (x[1] - 1)};
};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0Vector<R2> Vh{p_mesh, 2}; DiscreteFunctionP0Vector<R2> Vh{p_mesh, 2};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
auto vector = vector_affine(xj[cell_id]); Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
Vh[cell_id][0][0] = vector[0]; Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
Vh[cell_id][0][1] = vector[1];
Vh[cell_id][1][0] = vector[2];
Vh[cell_id][1][1] = vector[3];
}); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>(); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const R2>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
auto vector = vector_affine(xj[cell_id]);
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[0] - vector[0]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[1] - vector[1]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[0] - vector[2]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[1] - vector[3]));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
const R4 slope{+1.7, 0, -2.3, 0};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R2 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R2{0.1, 0} + xj[cell_id]) -
dpk_Vh(cell_id, 0)(xj[cell_id] - R2{0.1, 0}));
const R2 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R2{0.1, 0} + xj[cell_id]) -
dpk_Vh(cell_id, 1)(xj[cell_id] - R2{0.1, 0}));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[0] - slope[2]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[1] - slope[3]));
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_y_slope_error = 0;
const R4 slope{0, -0.6, 0, 1.1};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R2 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R2{0, 0.1} + xj[cell_id]) -
dpk_Vh(cell_id, 0)(xj[cell_id] - R2{0, 0.1}));
const R2 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R2{0, 0.1} + xj[cell_id]) -
dpk_Vh(cell_id, 1)(xj[cell_id] - R2{0, 0.1}));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[0] - slope[2]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[1] - slope[3]));
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
} }
} }
} }
...@@ -1379,174 +1112,47 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]") ...@@ -1379,174 +1112,47 @@ TEST_CASE("PolynomialReconstruction_degree_1", "[scheme]")
auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto R3_affine = [](const R3& x) -> R3 { auto R3_exact = [](const R3& x) -> R3 { return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)}; };
return R3{2.3 * (x[0] - 2), -1.3 * (x[1] - 1), 1.4 * (x[2] - 1)};
};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0<R3> fh{p_mesh}; DiscreteFunctionP0<R3> fh{p_mesh};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R3_affine(xj[cell_id]); }); mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R3_exact(xj[cell_id]); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(fh); auto reconstructions = PolynomialReconstruction{descriptor}.build(fh);
auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>(); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_fh, std::function(R3_exact));
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error =
std::max(max_mean_error, l2Norm(dpk_fh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id])));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope =
1. / 0.2 *
(dpk_fh[cell_id](R3{0.1, 0, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0.1, 0, 0}));
max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{2.3, 0, 0}));
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_y_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope =
1 / 0.2 *
(dpk_fh[cell_id](R3{0, 0.1, 0} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0.1, 0}));
max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - (R3{0, -1.3, 0})));
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12));
} }
SECTION("vector of R3")
{ {
double max_z_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope =
1 / 0.2 *
(dpk_fh[cell_id](R3{0, 0, 0.1} + xj[cell_id]) - dpk_fh[cell_id](xj[cell_id] - R3{0, 0, 0.1}));
max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - (R3{0, 0, 1.4})));
}
REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12));
}
}
SECTION("vector of R2")
{
using R6 = TinyVector<6>;
auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>(); auto p_mesh = MeshDataBaseForTests::get().hybrid3DMesh()->get<Mesh<3>>();
auto& mesh = *p_mesh; auto& mesh = *p_mesh;
auto vector_affine = [](const R3& x) -> R6 { std::array<std::function<R3(const R3&)>, 2> vector_exact //
return R6{+1.7 * (x[0] - 2), // = {[](const R3& x) -> R3 { return R3{+1.7 * (x[0] - 2), -0.6 * (x[1] - 1), +1.2 * (x[2] - 1)}; },
-0.6 * (x[1] - 1), // [](const R3& x) -> R3 { return R3{-2.3 * (x[0] - 2), +1.1 * (x[1] - 1), -0.3 * (x[2] - 1)}; }};
+1.2 * (x[2] - 1), //
-2.3 * (x[0] - 2), //
+1.1 * (x[1] - 1), //
-0.3 * (x[2] - 1)};
};
auto xj = MeshDataManager::instance().getMeshData(mesh).xj(); auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2}; DiscreteFunctionP0Vector<R3> Vh{p_mesh, 2};
parallel_for( parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
auto vector = vector_affine(xj[cell_id]); Vh[cell_id][0] = vector_exact[0](xj[cell_id]);
Vh[cell_id][0][0] = vector[0]; Vh[cell_id][1] = vector_exact[1](xj[cell_id]);
Vh[cell_id][0][1] = vector[1];
Vh[cell_id][0][2] = vector[2];
Vh[cell_id][1][0] = vector[3];
Vh[cell_id][1][1] = vector[4];
Vh[cell_id][1][2] = vector[5];
}); });
auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh); auto reconstructions = PolynomialReconstruction{descriptor}.build(Vh);
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>(); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const R3>>();
{ double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
double max_mean_error = 0; REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
auto vector = vector_affine(xj[cell_id]);
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[0] - vector[0]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[1] - vector[1]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 0)(xj[cell_id])[2] - vector[2]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[0] - vector[3]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[1] - vector[4]));
max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, 1)(xj[cell_id])[2] - vector[5]));
}
REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14));
}
{
double max_x_slope_error = 0;
const R6 slope{+1.7, 0, 0, -2.3, 0, 0};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0.1, 0, 0} + xj[cell_id]) -
dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0.1, 0, 0}));
const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0.1, 0, 0} + xj[cell_id]) -
dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0.1, 0, 0}));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
}
REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_y_slope_error = 0;
const R6 slope{0, -0.6, 0, 0, 1.1, 0};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0, 0.1, 0} + xj[cell_id]) -
dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0, 0.1, 0}));
const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0, 0.1, 0} + xj[cell_id]) -
dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0, 0.1, 0}));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
}
REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13));
}
{
double max_z_slope_error = 0;
const R6 slope{0, 0, 1.2, 0, 0, -0.3};
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const R3 reconstructed_slope0 = (1 / 0.2) * (dpk_Vh(cell_id, 0)(R3{0, 0, 0.1} + xj[cell_id]) -
dpk_Vh(cell_id, 0)(xj[cell_id] - R3{0, 0, 0.1}));
const R3 reconstructed_slope1 = (1 / 0.2) * (dpk_Vh(cell_id, 1)(R3{0, 0, 0.1} + xj[cell_id]) -
dpk_Vh(cell_id, 1)(xj[cell_id] - R3{0, 0, 0.1}));
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[0] - slope[0]));
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[1] - slope[1]));
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope0[2] - slope[2]));
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[0] - slope[3]));
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[1] - slope[4]));
max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope1[2] - slope[5]));
}
REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13));
}
} }
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment