diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 0f3dadaadcb2bb955ab11b59d25dd00f77f0e71f..212afc5523875f4b25c8af5d03155e136629d299 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -682,6 +682,18 @@ SchemeModule::SchemeModule()
));
+#warning REMOVE AFTER TESTS FINISHED
+ this->_addBuiltinFunction("test_reconstruction",
+ std::function(
+
+ [](const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>&
+ discrete_function_variant_list,
+ uint64_t degree) -> void {
+ test_reconstruction(discrete_function_variant_list, degree);
+ }
+
+ ));
+
MathFunctionRegisterForVh{*this};
}
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index c20fbc897c060fa182eb30683b109160938d68c6..0d0c829c7ce0d7f91c635dd59ae7636717ab7e2b 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -270,17 +270,15 @@ _computeEjkMean(const QuadratureFormula<3>& quadrature,
}
template <typename NodeListT, size_t Dimension>
-void
-_computeEjkMean(const TinyVector<Dimension>& Xj,
- const NodeValue<const TinyVector<Dimension>>& xr,
- const NodeListT& node_list,
- const CellType& cell_type,
- const double Vi,
- const size_t degree,
- const size_t basis_dimension,
- SmallArray<double>& inv_Vj_wq_detJ_ek,
- SmallArray<double>& mean_of_ejk)
-{}
+void _computeEjkMean(const TinyVector<Dimension>& Xj,
+ const NodeValue<const TinyVector<Dimension>>& xr,
+ const NodeListT& node_list,
+ const CellType& cell_type,
+ const double Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_wq_detJ_ek,
+ SmallArray<double>& mean_of_ejk);
template <typename NodeListT>
void
@@ -306,6 +304,86 @@ _computeEjkMean(const TinyVector<1>& Xj,
}
}
+template <typename ConformTransformationT>
+void
+_computeEjkBoundaryMean(const QuadratureFormula<1>& quadrature,
+ const ConformTransformationT& T,
+ const TinyVector<2>& Xj,
+ const double Vi,
+ const size_t degree,
+ const size_t dimension,
+ SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ using Rd = TinyVector<2>;
+
+ const double velocity_perp_e1 = T.velocity()[1];
+
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const TinyVector<1> xi_q = quadrature.point(i_q);
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+
+ {
+ size_t k = 0;
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1 / Vi;
+ for (; k <= degree; ++k) {
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
+ x_xj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * ((1. * k) / (k + 1));
+ }
+
+ for (size_t i_y = 1; i_y <= degree; ++i_y) {
+ const size_t begin_i_y_1 = ((i_y - 1) * (2 * degree - i_y + 4)) / 2;
+ for (size_t l = 0; l <= degree - i_y; ++l) {
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
+ }
+ }
+ }
+
+ for (size_t k = 1; k < dimension; ++k) {
+ mean_of_ejk[k - 1] += inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
+ }
+ }
+}
+
+void
+_computeEjkMeanByBoundary(const TinyVector<2>& Xj,
+ const CellId& cell_id,
+ const NodeValue<const TinyVector<2>>& xr,
+ const auto& cell_to_face_matrix,
+ const auto& face_to_node_matrix,
+ const auto& cell_face_is_reversed,
+ const CellValue<const double>& Vi,
+ const size_t degree,
+ const size_t basis_dimension,
+ SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek,
+ SmallArray<double>& mean_of_ejk)
+{
+ const auto& quadrature = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(degree + 1));
+
+ mean_of_ejk.fill(0);
+ auto cell_face_list = cell_to_face_matrix[cell_id];
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ bool is_reversed = cell_face_is_reversed[cell_id][i_face];
+
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ if (is_reversed) {
+ const LineTransformation<2> T{xr[face_node_list[1]], xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{xr[face_node_list[0]], xr[face_node_list[1]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, Vi[cell_id], degree, basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_of_ejk);
+ }
+ }
+}
+
template <typename NodeListT>
void
_computeEjkMean(const TinyVector<2>& Xj,
@@ -481,6 +559,7 @@ PolynomialReconstruction::_build(
auto cell_is_owned = mesh.connectivity().cellIsOwned();
auto cell_type = mesh.connectivity().cellType();
auto cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+ auto cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
const size_t max_stencil_size = [&]() {
size_t max_size = 0;
@@ -509,6 +588,9 @@ PolynomialReconstruction::_build(
SmallArray<SmallArray<double>> mean_j_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
SmallArray<SmallArray<double>> mean_i_of_ejk_pool(Kokkos::DefaultExecutionSpace::concurrency());
+ SmallArray<SmallArray<double>> inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool(Kokkos::DefaultExecutionSpace::concurrency());
+ SmallArray<SmallArray<double>> mean_l_of_ejk_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);
@@ -518,15 +600,17 @@ PolynomialReconstruction::_build(
inv_Vj_wq_detJ_ek_pool[i] = SmallArray<double>(basis_dimension);
mean_j_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
mean_i_of_ejk_pool[i] = SmallArray<double>(basis_dimension - 1);
+
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[i] = SmallArray<double>(basis_dimension);
}
#warning remove after rework
- const bool FORCE_INTEGRATE = []() {
+ [[maybe_unused]] const int INTEGRATE = []() {
auto value = std::getenv("INTEGRATE");
if (value == nullptr) {
- return false;
+ return 0;
} else {
- return true;
+ return std::atoi(value);
}
}();
@@ -596,7 +680,7 @@ PolynomialReconstruction::_build(
ShrinkMatrixView A(A_pool[t], stencil_cell_list.size());
- if ((m_descriptor.degree() == 1) and (not FORCE_INTEGRATE)) {
+ if ((m_descriptor.degree() == 1) and (INTEGRATE == 0)) {
const Rd& Xj = xj[cell_j_id];
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
const CellId cell_i_id = stencil_cell_list[i];
@@ -605,26 +689,59 @@ PolynomialReconstruction::_build(
A(i, l) = Xi_Xj[l];
}
}
- } else {
- SmallArray<double>& inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_ek_pool[t];
- SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
- SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
- const Rd& Xj = xj[cell_j_id];
+ } else if ((INTEGRATE == 2) or (INTEGRATE == 1)) {
+ if ((INTEGRATE == 2) and (MeshType::Dimension == 2)) {
+ if constexpr (MeshType::Dimension == 2) {
+ SmallArray<double>& inv_Vj_alpha_p_1_wq_X_prime_orth_ek = inv_Vj_alpha_p_1_wq_X_prime_orth_ek_pool[t];
+ SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
- _computeEjkMean(Xj, xr, cell_to_node_matrix[cell_j_id], cell_type[cell_j_id], Vj[cell_j_id],
- m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_j_of_ejk);
+ auto face_to_node_matrix = p_mesh->connectivity().faceToNodeMatrix();
+ auto cell_face_is_reversed = p_mesh->connectivity().cellFaceIsReversed();
+ const Rd& Xj = xj[cell_j_id];
- for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
- const CellId cell_i_id = stencil_cell_list[i];
+ _computeEjkMeanByBoundary(Xj, cell_j_id, xr, cell_to_face_matrix, face_to_node_matrix,
+ cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_j_of_ejk);
- _computeEjkMean(Xj, xr, cell_to_node_matrix[cell_i_id], cell_type[cell_i_id], Vj[cell_i_id],
- m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_i_of_ejk);
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
- for (size_t l = 0; l < basis_dimension - 1; ++l) {
- A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ _computeEjkMeanByBoundary(Xj, cell_i_id, xr, cell_to_face_matrix, face_to_node_matrix,
+ cell_face_is_reversed, Vj, m_descriptor.degree(), basis_dimension,
+ inv_Vj_alpha_p_1_wq_X_prime_orth_ek, mean_i_of_ejk);
+
+ for (size_t l = 0; l < basis_dimension - 1; ++l) {
+ A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ }
+ }
+ } else {
+ throw UnexpectedError("invalid mesh dimension");
+ }
+ } else {
+ SmallArray<double>& inv_Vj_wq_detJ_ek = inv_Vj_wq_detJ_ek_pool[t];
+ SmallArray<double>& mean_j_of_ejk = mean_j_of_ejk_pool[t];
+ SmallArray<double>& mean_i_of_ejk = mean_i_of_ejk_pool[t];
+
+ const Rd& Xj = xj[cell_j_id];
+
+ _computeEjkMean(Xj, xr, cell_to_node_matrix[cell_j_id], cell_type[cell_j_id], Vj[cell_j_id],
+ m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_j_of_ejk);
+
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+
+ _computeEjkMean(Xj, xr, cell_to_node_matrix[cell_i_id], cell_type[cell_i_id], Vj[cell_i_id],
+ m_descriptor.degree(), basis_dimension, inv_Vj_wq_detJ_ek, mean_i_of_ejk);
+
+ for (size_t l = 0; l < basis_dimension - 1; ++l) {
+ A(i, l) = mean_i_of_ejk[l] - mean_j_of_ejk[l];
+ }
}
}
+ } else {
+ throw UnexpectedError("invalid integration strategy");
}
if (m_descriptor.rowWeighting()) {
@@ -679,12 +796,6 @@ PolynomialReconstruction::_build(
Givens::solveCollectionInPlace(A, X, B);
}
- // Results are now written from the {ejk} basis to the {ek} basis
- if (m_descriptor.degree() > 1) {
- throw NotImplementedError(
- "write the result from the {ejk} basis to the {ek} basis to get correct polynomial values");
- }
-
column_begin = 0;
for (size_t i_dpk_variant = 0; i_dpk_variant < mutable_discrete_function_dpk_variant_list.size();
++i_dpk_variant) {
diff --git a/src/scheme/test_reconstruction.cpp b/src/scheme/test_reconstruction.cpp
index 3d255a5474bf9319a9079340fc7c20ab7ae7f834..4ba54bda2dc6a586e14e16fb01afee0417bb2476 100644
--- a/src/scheme/test_reconstruction.cpp
+++ b/src/scheme/test_reconstruction.cpp
@@ -4,42 +4,17 @@
#include <utils/Timer.hpp>
void
-test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list)
+test_reconstruction(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+ uint64_t degree)
{
- PolynomialReconstructionDescriptor descriptor{1};
- // {
- // std::cout << "** variable list one by one (50 times)\n";
- // Timer t;
- // for (auto discrete_function_v : discrete_function_variant_list) {
- // std::visit(
- // [&](auto&& discrete_function) {
- // auto mesh_v = discrete_function.meshVariant();
- // std::visit(
- // [&](auto&& p_mesh) {
- // using DiscreteFunctionT = std::decay_t<decltype(discrete_function)>;
- // PolynomialReconstruction reconstruction{descriptor};
- // if constexpr (is_discrete_function_P0_v<DiscreteFunctionT>) {
- // for (size_t i = 0; i < 50; ++i) {
- // auto res = reconstruction.build(*p_mesh, discrete_function);
- // }
- // } else if constexpr (is_discrete_function_P0_vector_v<DiscreteFunctionT>) {
- // std::cout << "Omitting discrete function p0 vector!\n";
- // }
- // },
- // mesh_v->variant());
- // },
- // discrete_function_v->discreteFunction());
- // }
- // t.pause();
- // std::cout << "t = " << t << '\n';
- // }
-
+ PolynomialReconstructionDescriptor descriptor{degree};
{
- std::cout << "** variable list at once (50 times)\n";
+ const size_t nb_loops = 1;
+ std::cout << "** variable list at once (" << nb_loops << " times)\n";
PolynomialReconstruction reconstruction{descriptor};
Timer t;
- for (size_t i = 0; i < 50; ++i) {
+ for (size_t i = 0; i < nb_loops; ++i) {
auto res = reconstruction.build(discrete_function_variant_list);
}
t.pause();
diff --git a/src/scheme/test_reconstruction.hpp b/src/scheme/test_reconstruction.hpp
index 6620acd44ca5e199ec7f4dcedde3a2d9afadda12..d2fc5bd031989eecdf4d837f566c137b3f874888 100644
--- a/src/scheme/test_reconstruction.hpp
+++ b/src/scheme/test_reconstruction.hpp
@@ -8,6 +8,7 @@
#include <vector>
void test_reconstruction(
- const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list);
+ const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list,
+ uint64_t degree = 1);
#endif // TEST_RECONSTRUCTION_HPP