diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp index 788c6df649c37090a22efcfbd8a202be22d45cd9..e3fe0193092fa5daccda3009c5c65d97b5afb31e 100644 --- a/src/scheme/PolynomialReconstruction.cpp +++ b/src/scheme/PolynomialReconstruction.cpp @@ -103,6 +103,8 @@ template <MeshConcept MeshType> [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> PolynomialReconstruction::_build( const size_t degree, + const bool preconditioning, + const bool row_weighting, const std::shared_ptr<const MeshType>& p_mesh, const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const { @@ -194,11 +196,13 @@ PolynomialReconstruction::_build( SmallArray<SmallMatrix<double>> A_pool(Kokkos::DefaultExecutionSpace::concurrency()); SmallArray<SmallMatrix<double>> B_pool(Kokkos::DefaultExecutionSpace::concurrency()); + SmallArray<SmallArray<double>> G_pool(Kokkos::DefaultExecutionSpace::concurrency()); SmallArray<SmallMatrix<double>> X_pool(Kokkos::DefaultExecutionSpace::concurrency()); for (size_t i = 0; i < A_pool.size(); ++i) { A_pool[i] = SmallMatrix<double>(max_stencil_size, basis_dimension - 1); B_pool[i] = SmallMatrix<double>(max_stencil_size, number_of_columns); + G_pool[i] = SmallArray<double>(basis_dimension - 1); X_pool[i] = SmallMatrix<double>(basis_dimension - 1, number_of_columns); } @@ -278,19 +282,55 @@ PolynomialReconstruction::_build( } } - for (size_t i = 0; i < stencil_cell_list.size(); ++i) { - const CellId cell_i_id = stencil_cell_list[i]; - const double weight = 1. / l2Norm(xj[cell_i_id] - Xj); - for (size_t l = 0; l < basis_dimension - 1; ++l) { - A(i, l) *= weight; - } - for (size_t l = 0; l < number_of_columns; ++l) { - B(i, l) *= weight; + if (row_weighting) { + // Add row weighting (give more importance to cells that are + // closer to j) + for (size_t i = 0; i < stencil_cell_list.size(); ++i) { + const CellId cell_i_id = stencil_cell_list[i]; + const double weight = 1. / l2Norm(xj[cell_i_id] - Xj); + for (size_t l = 0; l < basis_dimension - 1; ++l) { + A(i, l) *= weight; + } + for (size_t l = 0; l < number_of_columns; ++l) { + B(i, l) *= weight; + } } } const SmallMatrix<double>& X = X_pool[t]; - Givens::solveCollectionInPlace(A, X, B); + + if (preconditioning) { + // Add column weighting preconditioning (increase the presition) + SmallArray<double>& G = G_pool[t]; + + for (size_t l = 0; l < basis_dimension - 1; ++l) { + double g = 0; + for (size_t i = 0; i < stencil_cell_list.size(); ++i) { + const double Ail = A(i, l); + + g += Ail * Ail; + } + G[l] = std::sqrt(g); + } + + for (size_t l = 0; l < basis_dimension - 1; ++l) { + const double Gl = G[l]; + for (size_t i = 0; i < stencil_cell_list.size(); ++i) { + A(i, l) *= Gl; + } + } + + Givens::solveCollectionInPlace(A, X, B); + + for (size_t l = 0; l < basis_dimension - 1; ++l) { + const double Gl = G[l]; + for (size_t i = 0; i < number_of_columns; ++i) { + X(l, i) *= Gl; + } + } + } else { + Givens::solveCollectionInPlace(A, X, B); + } column_begin = 0; for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size(); @@ -408,6 +448,8 @@ PolynomialReconstruction::_build( std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> PolynomialReconstruction::build( const size_t degree, + const bool preconditioning, + const bool row_weighting, const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const { if (not hasSameMesh(discrete_function_variant_list)) { @@ -416,6 +458,9 @@ PolynomialReconstruction::build( auto mesh_v = getCommonMesh(discrete_function_variant_list); - return std::visit([&](auto&& p_mesh) { return this->_build(degree, p_mesh, discrete_function_variant_list); }, - mesh_v->variant()); + return std::visit( + [&](auto&& p_mesh) { + return this->_build(degree, preconditioning, row_weighting, p_mesh, discrete_function_variant_list); + }, + mesh_v->variant()); } diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp index 78f7bc5bad29c8da74249c21647d62763e9e6ee9..74bb827854facd86fc7e51191468908b72d6b481 100644 --- a/src/scheme/PolynomialReconstruction.hpp +++ b/src/scheme/PolynomialReconstruction.hpp @@ -29,6 +29,8 @@ class PolynomialReconstruction template <MeshConcept MeshType> [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> _build( const size_t degree, + const bool preconditioning, + const bool row_weighting, const std::shared_ptr<const MeshType>& p_mesh, const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const; @@ -265,11 +267,14 @@ class PolynomialReconstruction [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> build( const size_t degree, + const bool preconditioning, + const bool row_weighting, const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const; +#warning add reconstruction descriptor/option to handle the options to be used for reconstruction template <typename... DiscreteFunctionT> [[nodiscard]] std::vector<std::shared_ptr<const DiscreteFunctionDPkVariant>> - build(const size_t degree, DiscreteFunctionT... input) const + build(const size_t degree, const bool preconditioning, const bool row_weighting, DiscreteFunctionT... input) const { std::vector<std::shared_ptr<const DiscreteFunctionVariant>> variant_vector; auto convert = [&variant_vector](auto&& df) { @@ -289,7 +294,7 @@ class PolynomialReconstruction }; (convert(std::forward<DiscreteFunctionT>(input)), ...); - return this->build(degree, variant_vector); + return this->build(degree, preconditioning, row_weighting, variant_vector); } PolynomialReconstruction() {} diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp index 87b9bcc16193ed00a3d505f859c12c7ee14da214..1267aabecc233fc3b042b026e2dd6897afe2b1c6 100644 --- a/src/scheme/test_reconstruction.cpp +++ b/src/scheme/test_reconstruction.cpp @@ -87,7 +87,7 @@ test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVari PolynomialReconstruction reconstruction; Timer t; for (size_t i = 0; i < 50; ++i) { - auto res = reconstruction.build(1, discrete_function_variant_list); + auto res = reconstruction.build(1, true, true, discrete_function_variant_list); } t.pause(); std::cout << "t = " << t << '\n'; diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp index e019cc0162e0ccf52a583020f8d8a8c9346809b9..ada21453a8266a90ad0749613fc39dcc6629ca17 100644 --- a/tests/test_PolynomialReconstruction.cpp +++ b/tests/test_PolynomialReconstruction.cpp @@ -43,7 +43,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, fh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, fh); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>(); @@ -53,7 +53,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -64,7 +64,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7)); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14)); } } } @@ -92,7 +92,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, uh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3>>(); @@ -102,7 +102,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -113,7 +113,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1})); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14)); } } } @@ -143,7 +143,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R3x3_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, Ah); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, Ah); auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>(); @@ -153,7 +153,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -169,7 +169,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, // frobeniusNorm(reconstructed_slope - slops)); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13)); } } } @@ -200,7 +200,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") } }); - auto reconstructions = PolynomialReconstruction{}.build(1, Vh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>(); @@ -212,7 +212,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i])); } } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -226,8 +226,8 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i])); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13)); } + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13)); } } } @@ -287,9 +287,10 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") } }); - auto reconstructions = PolynomialReconstruction{}.build(1, std::make_shared<DiscreteFunctionVariant>(fh), - uh, std::make_shared<DiscreteFunctionP0<R3x3>>(Ah), - DiscreteFunctionVariant(Vh)); + auto reconstructions = + PolynomialReconstruction{}.build(1, true, true, std::make_shared<DiscreteFunctionVariant>(fh), uh, + std::make_shared<DiscreteFunctionP0<R3x3>>(Ah), + DiscreteFunctionVariant(Vh)); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<1, const double>>(); @@ -299,7 +300,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -310,7 +311,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7)); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14)); } auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>(); @@ -321,7 +322,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -332,7 +333,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1})); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-14)); } auto dpk_Ah = reconstructions[2]->get<DiscreteFunctionDPk<1, const R3x3>>(); @@ -343,7 +344,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R3x3_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -359,7 +360,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, // frobeniusNorm(reconstructed_slope - slops)); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13)); } auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>(); @@ -372,7 +373,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i])); } } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -386,8 +387,8 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - slope[i])); } - REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-13)); } + REQUIRE(parallel::allReduceMax(max_slope_error) == Catch::Approx(0).margin(1E-13)); } } } @@ -414,7 +415,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, fh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, fh); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<2, const double>>(); @@ -424,7 +425,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -435,7 +436,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7)); } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-14)); } { @@ -446,7 +447,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3))); } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-14)); } } } @@ -474,7 +475,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, uh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<2, const R3>>(); @@ -484,7 +485,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -495,7 +496,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1})); } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-14)); } { @@ -506,7 +507,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1})); } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-14)); } } } @@ -533,7 +534,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, Ah); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, Ah); auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<2, const R2x2>>(); @@ -543,7 +544,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -556,7 +557,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, +2.1, // -0.6, -2.3})); } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12)); } { @@ -569,7 +570,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, // -2.1, +1.3})); } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12)); } } } @@ -601,7 +602,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") } }); - auto reconstructions = PolynomialReconstruction{}.build(1, Vh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<2, const double>>(); @@ -613,7 +614,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i])); } } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -627,7 +628,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i])); } } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12)); } { @@ -641,7 +642,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i])); } } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12)); } } } @@ -668,7 +669,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = R_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, fh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, fh); auto dpk_fh = reconstructions[0]->get<DiscreteFunctionDPk<3, const double>>(); @@ -678,7 +679,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_fh[cell_id](xj[cell_id]) - R_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -689,7 +690,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - 1.7)); } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-14)); } { @@ -700,7 +701,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - (-1.3))); } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-14)); } { @@ -711,7 +712,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - 2.1)); } - REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-14)); } } } @@ -737,7 +738,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { uh[cell_id] = R3_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, uh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, uh); auto dpk_uh = reconstructions[0]->get<DiscreteFunctionDPk<3, const R3>>(); @@ -747,7 +748,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, l2Norm(dpk_uh[cell_id](xj[cell_id]) - R3_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-14)); } { @@ -758,7 +759,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_x_slope_error = std::max(max_x_slope_error, l2Norm(reconstructed_slope - R3{1.7, -0.6, 3.1})); } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-13)); } { @@ -769,7 +770,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1})); } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-13)); } { @@ -780,7 +781,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_z_slope_error = std::max(max_z_slope_error, l2Norm(reconstructed_slope - R3{1.8, -3.7, 1.9})); } - REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-13)); } } } @@ -808,7 +809,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") parallel_for( mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { Ah[cell_id] = R2x2_affine(xj[cell_id]); }); - auto reconstructions = PolynomialReconstruction{}.build(1, Ah); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, Ah); auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<3, const R2x2>>(); @@ -818,7 +819,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, frobeniusNorm(dpk_Ah[cell_id](xj[cell_id]) - R2x2_affine(xj[cell_id]))); } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -830,7 +831,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_x_slope_error = std::max(max_x_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.7, 2.1, // -2.3, +3.1})); } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12)); } { @@ -843,7 +844,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") std::max(max_y_slope_error, frobeniusNorm(reconstructed_slope - R2x2{+1.2, -2.2, // 1.3, +0.8})); } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12)); } { @@ -856,7 +857,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") std::max(max_z_slope_error, frobeniusNorm(reconstructed_slope - R2x2{-1.3, -2.4, // +1.4, -1.8})); } - REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12)); } } } @@ -889,7 +890,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") } }); - auto reconstructions = PolynomialReconstruction{}.build(1, Vh); + auto reconstructions = PolynomialReconstruction{}.build(1, true, true, Vh); auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<3, const double>>(); @@ -901,7 +902,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_mean_error = std::max(max_mean_error, std::abs(dpk_Vh(cell_id, i)(xj[cell_id]) - vector[i])); } } - REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-13)); + REQUIRE(parallel::allReduceMax(max_mean_error) == Catch::Approx(0).margin(1E-13)); } { @@ -915,7 +916,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_x_slope_error = std::max(max_x_slope_error, std::abs(reconstructed_slope - slope[i])); } } - REQUIRE(max_x_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_x_slope_error) == Catch::Approx(0).margin(1E-12)); } { @@ -929,7 +930,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_y_slope_error = std::max(max_y_slope_error, std::abs(reconstructed_slope - slope[i])); } } - REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_y_slope_error) == Catch::Approx(0).margin(1E-12)); } { @@ -943,7 +944,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") max_z_slope_error = std::max(max_z_slope_error, std::abs(reconstructed_slope - slope[i])); } } - REQUIRE(max_z_slope_error == Catch::Approx(0).margin(1E-12)); + REQUIRE(parallel::allReduceMax(max_z_slope_error) == Catch::Approx(0).margin(1E-12)); } } } @@ -958,7 +959,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]") auto p_mesh2 = MeshDataBaseForTests::get().cartesian1DMesh()->get<Mesh<1>>(); DiscreteFunctionP0<double> f2{p_mesh2}; - REQUIRE_THROWS_WITH(PolynomialReconstruction{}.build(1, f1, f2), + REQUIRE_THROWS_WITH(PolynomialReconstruction{}.build(1, true, true, f1, f2), "error: cannot reconstruct functions living of different meshes simultaneously"); } }