diff --git a/CMakeLists.txt b/CMakeLists.txt
index bb8832a50565bde6f37ab589737a15e08f116030..46126f9121cfafeb028c9d6f32a2ddf826c9fcf1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,6 +3,10 @@ cmake_minimum_required (VERSION 3.19)
# CMake utils
list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
+set(CMAKE_INSTALL_LIBDIR "lib")
+set(CMAKE_INSTALL_INCLUDEDIR "include")
+set(CMAKE_INSTALL_BINDIR "bin")
+
# Forbids in-source builds
include(CheckNotInSources)
diff --git a/src/geometry/SquareTransformation.hpp b/src/geometry/SquareTransformation.hpp
index d6474f33e4ecbaa4bcfba525dadc9fd38e95962d..b0e970c6eda673b48c5cf5c0d0a1ad7e73f8ef74 100644
--- a/src/geometry/SquareTransformation.hpp
+++ b/src/geometry/SquareTransformation.hpp
@@ -84,6 +84,25 @@ class SquareTransformation
return m_coefficients * X + m_shift;
}
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation(const TinyVector<2>& X) const
+ {
+ const auto& T = m_coefficients;
+ const double& x = X[0];
+ const double& y = X[1];
+
+ const TinyVector<Dimension> dxT{T(0, 0) + T(0, 2) * y, //
+ T(1, 0) + T(1, 2) * y, //
+ T(2, 0) + T(2, 2) * y};
+
+ const TinyVector<Dimension> dyT{T(0, 1) + T(0, 2) * x, //
+ T(1, 1) + T(1, 2) * x, //
+ T(2, 1) + T(2, 2) * x};
+
+ return crossProduct(dxT, dyT);
+ }
+
PUGS_INLINE double
areaVariationNorm(const TinyVector<2>& X) const
{
diff --git a/src/geometry/TriangleTransformation.hpp b/src/geometry/TriangleTransformation.hpp
index df8444d894a7ed56c0fe3825ffbf1cc14a742531..c49db368684d8c96aeb5061041828799e8509fac 100644
--- a/src/geometry/TriangleTransformation.hpp
+++ b/src/geometry/TriangleTransformation.hpp
@@ -66,6 +66,7 @@ class TriangleTransformation
private:
TinyMatrix<Dimension, 2> m_jacobian;
TinyVector<Dimension> m_shift;
+ TinyVector<Dimension> m_area_variation;
double m_area_variation_norm;
public:
@@ -76,6 +77,21 @@ class TriangleTransformation
return m_jacobian * x + m_shift;
}
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation() const
+ {
+ return m_area_variation;
+ }
+
+ PUGS_INLINE
+ TinyVector<Dimension>
+ areaVariation(const TinyVector<2>&) const
+ {
+ return m_area_variation;
+ }
+
+ PUGS_INLINE
double
areaVariationNorm() const
{
@@ -92,7 +108,8 @@ class TriangleTransformation
m_shift = a;
- m_area_variation_norm = l2Norm(crossProduct(b - a, c - a));
+ m_area_variation = crossProduct(b - a, c - a);
+ m_area_variation_norm = l2Norm(m_area_variation);
}
~TriangleTransformation() = default;
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index c646dc4fbec060bf5983a24ec33ca8bdb2718725..923458cde39814ef9c1a841cf0cb30cf5eacbbb7 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -824,11 +824,7 @@ PolynomialReconstruction::_build(
return this->_build<CellCenterReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
}
case IntegrationMethodType::boundary: {
- if constexpr (MeshType::Dimension == 2) {
- return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh,
- discrete_function_variant_list);
- }
- [[fallthrough]];
+ return this->_build<BoundaryIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
}
case IntegrationMethodType::element: {
return this->_build<ElementIntegralReconstructionMatrixBuilder<MeshType>>(p_mesh, discrete_function_variant_list);
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
index d030ac6d51d22d70eb209fa437a8684064d57529..7a9f7eded202a3cb98b74769cf4d23930a32c60f 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.cpp
@@ -1,20 +1,183 @@
#include <scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
+#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineCubicTransformation.hpp>
#include <geometry/LineParabolicTransformation.hpp>
+#include <geometry/LineTransformation.hpp>
+#include <geometry/SquareTransformation.hpp>
#include <geometry/SymmetryUtils.hpp>
-#include <mesh/Connectivity.hpp>
+#include <geometry/TriangleTransformation.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <mesh/PolynomialMesh.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
template <MeshConcept MeshTypeT>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::SpecificDimensionalData
+{
+ private:
+ const ItemToItemMatrix<ItemType::cell, ItemType::face> m_cell_to_face_matrix;
+ const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix;
+ const FaceValuePerCell<const bool> m_cell_face_is_reversed;
+
+ public:
+ PUGS_INLINE const auto&
+ cellToFaceMatrix() const noexcept
+ {
+ return m_cell_to_face_matrix;
+ }
+
+ PUGS_INLINE
+ const auto&
+ faceToNodeMatrix() const noexcept
+ {
+ return m_face_to_node_matrix;
+ }
+
+ PUGS_INLINE
+ const auto&
+ cellFaceIsReversed() const noexcept
+ {
+ return m_cell_face_is_reversed;
+ }
+
+ SpecificDimensionalData(const Connectivity<MeshType::Dimension>& connectivity)
+ : m_cell_to_face_matrix{connectivity.cellToFaceMatrix()},
+ m_face_to_node_matrix{connectivity.faceToNodeMatrix()},
+ m_cell_face_is_reversed{connectivity.cellFaceIsReversed()}
+ {}
+
+ ~SpecificDimensionalData() = default;
+};
+
+template <>
+class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::SpecificDimensionalData
+{
+ private:
+ const ItemToItemMatrix<ItemType::cell, ItemType::node> m_cell_to_node_matrix;
+
+ public:
+ PUGS_INLINE
+ const auto&
+ cellToNodeMatrix() const noexcept
+ {
+ return m_cell_to_node_matrix;
+ }
+
+ SpecificDimensionalData(const Connectivity<1>& connectivity) : m_cell_to_node_matrix{connectivity.cellToNodeMatrix()}
+ {}
+
+ ~SpecificDimensionalData() = default;
+};
+
+template <>
+template <typename ConformTransformationT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<3>>::_computeEjkBoundaryMean(
+ const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double inv_Vi,
+ SmallArray<double>& mean_of_ejk)
+{
+ for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
+ const double wq = quadrature.weight(i_q);
+ const TinyVector<2> xi_q = quadrature.point(i_q);
+
+ const double area_variation_e1 = T.areaVariation(xi_q)[0] * inv_Vi;
+
+ const Rd X_Xj = T(xi_q) - Xj;
+
+ const double x_xj = X_Xj[0];
+ const double y_yj = X_Xj[1];
+ const double z_zj = X_Xj[2];
+
+ {
+ size_t k = 0;
+ m_tmp_array[k++] = x_xj * wq * area_variation_e1;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+ const size_t begin_i_y_1 = m_yz_row_index[i_y - 1];
+ const size_t nb_monoms = m_yz_row_size[i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_tmp_array[k] = y_yj * m_tmp_array[begin_i_y_1 + l];
+ }
+ }
+
+ for (size_t i_z = 1; i_z <= m_polynomial_degree; ++i_z) {
+ const size_t nb_y = m_yz_row_size[m_z_triangle_index[i_z]];
+ const size_t index_z = m_z_triangle_index[i_z];
+ const size_t index_z_1 = m_z_triangle_index[i_z - 1];
+ for (size_t i_y = 0; i_y < nb_y; ++i_y) {
+ const size_t begin_i_yz_1 = m_yz_row_index[index_z_1 + i_y];
+ const size_t nb_monoms = m_yz_row_size[index_z + i_y];
+ for (size_t l = 0; l < nb_monoms; ++l, ++k) {
+ m_tmp_array[k] = z_zj * m_tmp_array[begin_i_yz_1 + l];
+ }
+ }
+ }
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] += m_tmp_array[k];
+ }
+ }
+}
+
+template <>
template <typename ConformTransformationT>
void
-PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::_computeEjkBoundaryMean(
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::_computeEjkBoundaryMean(
+ const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
+ const ConformTransformationT& T,
+ const Rd& Xj,
+ const double inv_Vi,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double velocity_perp_e1 = T.velocity()[1] * inv_Vi;
+
+ 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;
+ m_tmp_array[k++] = x_xj * wq * velocity_perp_e1;
+ for (; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
+ const size_t begin_i_y_1 = m_y_row_index[i_y - 1];
+ for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
+ m_tmp_array[k++] = y_yj * m_tmp_array[begin_i_y_1 + l];
+ }
+ }
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] += m_tmp_array[k];
+ }
+ }
+}
+
+template <>
+template <typename ConformTransformationT>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<PolynomialMesh<2>>::_computeEjkBoundaryMean(
const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
const ConformTransformationT& T,
const Rd& Xj,
@@ -34,25 +197,26 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const double y_yj = X_Xj[1];
{
- size_t k = 0;
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1;
+ size_t k = 0;
+ m_tmp_array[k++] = x_xj * wq * velocity_perp_e1;
for (; k <= m_polynomial_degree; ++k) {
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
- x_xj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * m_antiderivative_coef[k];
+ m_tmp_array[k] //
+ = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
}
for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
- const size_t begin_i_y_1 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
+ const size_t begin_i_y_1 = m_y_row_index[i_y - 1];
for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
+ m_tmp_array[k++] = y_yj * m_tmp_array[begin_i_y_1 + l];
}
}
}
for (size_t k = 1; k < m_basis_dimension; ++k) {
- mean_of_ejk[k - 1] += m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
+ mean_of_ejk[k - 1] += m_tmp_array[k];
}
}
+
} else {
for (size_t i_q = 0; i_q < quadrature.numberOfPoints(); ++i_q) {
const double wq = quadrature.weight(i_q);
@@ -66,23 +230,22 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const double y_yj = X_Xj[1];
{
- size_t k = 0;
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = x_xj * wq * velocity_perp_e1;
+ size_t k = 0;
+ m_tmp_array[k++] = x_xj * wq * velocity_perp_e1;
for (; k <= m_polynomial_degree; ++k) {
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k] =
- x_xj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k - 1] * m_antiderivative_coef[k];
+ m_tmp_array[k] = x_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k];
}
for (size_t i_y = 1; i_y <= m_polynomial_degree; ++i_y) {
const size_t begin_i_y_1 = ((i_y - 1) * (2 * m_polynomial_degree - i_y + 4)) / 2;
for (size_t l = 0; l <= m_polynomial_degree - i_y; ++l) {
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k++] = y_yj * m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[begin_i_y_1 + l];
+ m_tmp_array[k++] = y_yj * m_tmp_array[begin_i_y_1 + l];
}
}
}
for (size_t k = 1; k < m_basis_dimension; ++k) {
- mean_of_ejk[k - 1] += m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek[k];
+ mean_of_ejk[k - 1] += m_tmp_array[k];
}
}
}
@@ -96,24 +259,69 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
SmallArray<double>& mean_of_ejk)
{
if constexpr (is_polygonal_mesh_v<MeshType>) {
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
-
const double inv_Vi = 1. / m_Vj[cell_id];
+ const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data_ptr->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
+
mean_of_ejk.fill(0);
- auto cell_face_list = m_cell_to_face_matrix[cell_id];
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = m_face_to_node_matrix[face_id];
- if (is_reversed) {
- const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
- } else {
- const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ if constexpr (MeshType::Dimension == 2) {
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ if (face_is_reversed[i_face]) {
+ const LineTransformation<2> T{m_xr[face_node_list[1]], m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+
+ const auto& triangle_quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(m_polynomial_degree + 1));
+ const auto& square_quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ switch (face_node_list.size()) {
+ case 3: {
+ if (face_is_reversed[i_face]) {
+ const TriangleTransformation<3> T{m_xr[face_node_list[2]], m_xr[face_node_list[1]],
+ m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const TriangleTransformation<3> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]],
+ m_xr[face_node_list[2]]};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ case 4: {
+ if (face_is_reversed[i_face]) {
+ const SquareTransformation<3> T{m_xr[face_node_list[3]], m_xr[face_node_list[2]], m_xr[face_node_list[1]],
+ m_xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const SquareTransformation<3> T{m_xr[face_node_list[0]], m_xr[face_node_list[1]], m_xr[face_node_list[2]],
+ m_xr[face_node_list[3]]};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ throw NotImplementedError("invalid face type");
+ }
+ }
}
}
} else {
@@ -126,17 +334,31 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>
const double inv_Vi = 1. / m_Vj[cell_id];
+ const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data_ptr->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
+
auto xr = m_mesh.xr();
auto xl = m_mesh.xl();
mean_of_ejk.fill(0);
- auto cell_face_list = m_cell_to_face_matrix[cell_id];
+
for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
+ bool is_reversed = face_is_reversed[i_face];
const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = m_face_to_node_matrix[face_id];
-#warning rework
+ auto face_node_list = face_to_node_matrix[face_id];
+
switch (m_mesh.degree()) {
+ case 1: {
+ if (is_reversed) {
+ const LineTransformation<2> T{xr[face_node_list[1]], xr[face_node_list[0]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{xr[face_node_list[0]], xr[face_node_list[1]]};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
case 2: {
if (is_reversed) {
const LineParabolicTransformation<2> T{xr[face_node_list[1]], xl[face_id][0], xr[face_node_list[0]]};
@@ -178,29 +400,79 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const CellId& cell_id,
SmallArray<double>& mean_of_ejk)
{
- if constexpr (is_polygonal_mesh_v<MeshType>) {
- const auto& quadrature =
- QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+ const double inv_Vi = 1. / m_Vj[cell_id];
- const double inv_Vi = 1. / m_Vj[cell_id];
+ const auto& face_is_reversed = m_dimensional_data_ptr->cellFaceIsReversed()[cell_id];
+ const auto& cell_face_list = m_dimensional_data_ptr->cellToFaceMatrix()[cell_id];
+ const auto& face_to_node_matrix = m_dimensional_data_ptr->faceToNodeMatrix();
- mean_of_ejk.fill(0);
- auto cell_face_list = m_cell_to_face_matrix[cell_id];
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
+ mean_of_ejk.fill(0);
- const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = m_face_to_node_matrix[face_id];
+ if constexpr (is_polygonal_mesh_v<MeshType>) {
+ if constexpr (MeshType::Dimension == 2) {
+ const auto& quadrature =
+ QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
- const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
- const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
- if (is_reversed) {
- const LineTransformation<2> T{x1, x0};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
- } else {
- const LineTransformation<2> T{x0, x1};
- _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ if (face_is_reversed[i_face]) {
+ const LineTransformation<2> T{x1, x0};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{x0, x1};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ }
+ } else {
+ static_assert(MeshType::Dimension == 3);
+
+ const auto& triangle_quadrature =
+ QuadratureManager::instance().getTriangleFormula(GaussQuadratureDescriptor(m_polynomial_degree + 1));
+ const auto& square_quadrature =
+ QuadratureManager::instance().getSquareFormula(GaussLegendreQuadratureDescriptor(m_polynomial_degree + 1));
+
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const FaceId face_id = cell_face_list[i_face];
+ auto face_node_list = face_to_node_matrix[face_id];
+ switch (face_node_list.size()) {
+ case 3: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[2]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const TriangleTransformation<3> T{x2, x1, x0};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const TriangleTransformation<3> T{x0, x1, x2};
+ _computeEjkBoundaryMean(triangle_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ case 4: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[3]]);
+ const auto x1 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[2]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x3 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (face_is_reversed[i_face]) {
+ const SquareTransformation<3> T{x3, x2, x1, x0};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const SquareTransformation<3> T{x0, x1, x2, x3};
+ _computeEjkBoundaryMean(square_quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
+ default: {
+ throw NotImplementedError("invalid face type");
+ }
+ }
}
}
} else {
@@ -211,19 +483,28 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const auto& quadrature = QuadratureManager::instance().getLineFormula(
GaussLegendreQuadratureDescriptor((m_polynomial_degree + 1) * m_mesh.degree() + (m_mesh.degree() - 1)));
- const double inv_Vi = 1. / m_Vj[cell_id];
-
auto xl = m_mesh.xl();
- mean_of_ejk.fill(0);
- auto cell_face_list = m_cell_to_face_matrix[cell_id];
for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- bool is_reversed = m_cell_face_is_reversed[cell_id][i_face];
+ bool is_reversed = face_is_reversed[i_face];
const FaceId face_id = cell_face_list[i_face];
- auto face_node_list = m_face_to_node_matrix[face_id];
+ auto face_node_list = face_to_node_matrix[face_id];
switch (m_mesh.degree()) {
+ case 1: {
+ const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
+ const auto x2 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[0]]);
+
+ if (is_reversed) {
+ const LineTransformation<2> T{x2, x0};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ } else {
+ const LineTransformation<2> T{x0, x2};
+ _computeEjkBoundaryMean(quadrature, T, Xj, inv_Vi, mean_of_ejk);
+ }
+ break;
+ }
case 2: {
const auto x0 = symmetrize_coordinates(origin, normal, m_xr[face_node_list[1]]);
const auto x1 = symmetrize_coordinates(origin, normal, xl[face_id][0]);
@@ -263,13 +544,87 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
}
}
+template <>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::_computeEjkMeanByBoundary(
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& cell_node_list = m_dimensional_data_ptr->cellToNodeMatrix()[cell_id];
+
+ const double xr1_xj = (m_xr[cell_node_list[1]] - Xj)[0];
+
+ m_tmp_array[0] = xr1_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr1_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] = m_tmp_array[k];
+ }
+
+ const double xr0_xj = (m_xr[cell_node_list[0]] - Xj)[0];
+
+ m_tmp_array[0] = xr0_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr0_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] -= m_tmp_array[k];
+ }
+}
+
+template <>
+void
+PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<1>>::_computeEjkMeanByBoundaryInSymmetricCell(const Rd& origin,
+ const Rd& normal,
+ const Rd& Xj,
+ const CellId& cell_id,
+ SmallArray<double>& mean_of_ejk)
+{
+ const double inv_Vi = 1. / m_Vj[cell_id];
+
+ const auto& cell_node_list = m_dimensional_data_ptr->cellToNodeMatrix()[cell_id];
+
+ const double xr1_xj = (symmetrize_coordinates(origin, normal, m_xr[cell_node_list[0]]) - Xj)[0];
+
+ m_tmp_array[0] = xr1_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr1_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] = m_tmp_array[k];
+ }
+
+ const double xr0_xj = (symmetrize_coordinates(origin, normal, m_xr[cell_node_list[1]]) - Xj)[0];
+
+ m_tmp_array[0] = xr0_xj * inv_Vi;
+ for (size_t k = 1; k <= m_polynomial_degree; ++k) {
+ m_tmp_array[k] //
+ = xr0_xj * m_tmp_array[k - 1] * m_antiderivative_correction_coef[k]; // ((1. * k) / (k + 1));
+ }
+
+ for (size_t k = 1; k < m_basis_dimension; ++k) {
+ mean_of_ejk[k - 1] -= m_tmp_array[k];
+ }
+}
+
template <MeshConcept MeshTypeT>
void
PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<MeshTypeT>::build(
const CellId cell_j_id,
ShrinkMatrixView<SmallMatrix<double>>& A)
{
- if constexpr (MeshType::Dimension == 2) {
+ if constexpr (MeshType::Dimension <= 3) {
const auto& stencil_cell_list = m_stencil_array[cell_j_id];
const Rd& Xj = m_xj[cell_j_id];
@@ -319,12 +674,12 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
m_basis_dimension{
DiscreteFunctionDPk<MeshType::Dimension, double>::BasisViewType::dimensionFromDegree(polynomial_degree)},
m_polynomial_degree{polynomial_degree},
- m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek{m_basis_dimension},
+ m_tmp_array{m_basis_dimension},
m_mean_j_of_ejk{m_basis_dimension - 1},
m_mean_i_of_ejk{m_basis_dimension - 1},
- m_cell_to_face_matrix{mesh.connectivity().cellToFaceMatrix()},
- m_face_to_node_matrix{mesh.connectivity().faceToNodeMatrix()},
- m_cell_face_is_reversed{mesh.connectivity().cellFaceIsReversed()},
+
+ m_dimensional_data_ptr{std::make_shared<SpecificDimensionalData>(mesh.connectivity())},
+
m_stencil_array{stencil_array},
m_symmetry_origin_list{symmetry_origin_list},
m_symmetry_normal_list{symmetry_normal_list},
@@ -332,16 +687,59 @@ PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
m_xj{MeshDataManager::instance().getMeshData(mesh).xj()},
m_xr{mesh.xr()}
{
+ SmallArray<double> antiderivative_correction_coef(m_polynomial_degree + 1);
+ for (size_t k = 0; k < antiderivative_correction_coef.size(); ++k) {
+ // The antiderivative of x^k is k/(k+1) times the antiderivative of x^(k-1)
+ antiderivative_correction_coef[k] = ((1. * k) / (k + 1));
+ }
+
+ m_antiderivative_correction_coef = antiderivative_correction_coef;
+
if constexpr (MeshType::Dimension == 2) {
- SmallArray<double> antiderivative_coef(m_polynomial_degree + 1);
- for (size_t k = 0; k < antiderivative_coef.size(); ++k) {
- antiderivative_coef[k] = ((1. * k) / (k + 1));
+ SmallArray<size_t> y_row_index(m_polynomial_degree + 1);
+
+ size_t i_y = 0;
+
+ y_row_index[i_y++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n > 1; --n, ++i_y) {
+ y_row_index[i_y] = y_row_index[i_y - 1] + n;
}
- m_antiderivative_coef = antiderivative_coef;
+ m_y_row_index = y_row_index;
+
+ } else if constexpr (MeshType::Dimension == 3) {
+ SmallArray<size_t> yz_row_index((m_polynomial_degree + 2) * (m_polynomial_degree + 1) / 2 + 1);
+ SmallArray<size_t> z_triangle_index(m_polynomial_degree + 1);
+
+ {
+ size_t i_z = 0;
+ size_t i_yz = 0;
+
+ yz_row_index[i_yz++] = 0;
+ for (ssize_t n = m_polynomial_degree + 1; n >= 1; --n) {
+ z_triangle_index[i_z++] = i_yz - 1;
+ for (ssize_t i = n; i >= 1; --i) {
+ yz_row_index[i_yz] = yz_row_index[i_yz - 1] + i;
+ ++i_yz;
+ }
+ }
+ }
+
+ SmallArray<size_t> yz_row_size{yz_row_index.size() - 1};
+ for (size_t i = 0; i < yz_row_size.size(); ++i) {
+ yz_row_size[i] = yz_row_index[i + 1] - yz_row_index[i];
+ }
+
+ m_yz_row_index = yz_row_index;
+ m_z_triangle_index = z_triangle_index;
+ m_yz_row_size = yz_row_size;
}
}
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<1>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<2>>::build(
const CellId,
ShrinkMatrixView<SmallMatrix<double>>&);
@@ -350,6 +748,17 @@ template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuil
const CellId,
ShrinkMatrixView<SmallMatrix<double>>&);
+template void PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<Mesh<3>>::build(
+ const CellId,
+ ShrinkMatrixView<SmallMatrix<double>>&);
+
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<1>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
+
template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
Mesh<2>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
const size_t,
@@ -363,3 +772,9 @@ template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
const SmallArray<const Rd>&,
const SmallArray<const Rd>&,
const CellToCellStencilArray&);
+template PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder<
+ Mesh<3>>::BoundaryIntegralReconstructionMatrixBuilder(const MeshType&,
+ const size_t,
+ const SmallArray<const Rd>&,
+ const SmallArray<const Rd>&,
+ const CellToCellStencilArray&);
diff --git a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
index df24d76a73097d13f6379681780d90aed6f49e1d..60707189b909e9d5f17ae2852be655cbfb877351 100644
--- a/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/BoundaryIntegralReconstructionMatrixBuilder.hpp
@@ -3,14 +3,14 @@
#include <algebra/ShrinkMatrixView.hpp>
#include <algebra/SmallMatrix.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <geometry/LineTransformation.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/StencilArray.hpp>
-#include <mesh/SubItemValuePerItem.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>
+template <size_t Dimension>
+class QuadratureFormula;
+
template <MeshConcept MeshTypeT>
class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
{
@@ -26,13 +26,12 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
const size_t m_basis_dimension;
const size_t m_polynomial_degree;
- const SmallArray<double> m_inv_Vj_alpha_p_1_wq_X_prime_orth_ek;
+ const SmallArray<double> m_tmp_array;
SmallArray<double> m_mean_j_of_ejk;
SmallArray<double> m_mean_i_of_ejk;
- const ItemToItemMatrix<ItemType::cell, ItemType::face> m_cell_to_face_matrix;
- const ItemToItemMatrix<ItemType::face, ItemType::node> m_face_to_node_matrix;
- const FaceValuePerCell<const bool> m_cell_face_is_reversed;
+ class SpecificDimensionalData;
+ std::shared_ptr<SpecificDimensionalData> m_dimensional_data_ptr;
const CellToCellStencilArray& m_stencil_array;
@@ -43,7 +42,15 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
const CellValue<const Rd> m_xj;
const NodeValue<const Rd> m_xr;
- SmallArray<const double> m_antiderivative_coef;
+ // 2D
+ SmallArray<const size_t> m_y_row_index;
+
+ // 3D
+ SmallArray<const size_t> m_yz_row_index;
+ SmallArray<const size_t> m_z_triangle_index;
+ SmallArray<const size_t> m_yz_row_size;
+
+ SmallArray<const double> m_antiderivative_correction_coef;
template <typename ConformTransformationT>
void _computeEjkBoundaryMean(const QuadratureFormula<MeshType::Dimension - 1>& quadrature,
@@ -76,7 +83,6 @@ class PolynomialReconstruction::BoundaryIntegralReconstructionMatrixBuilder
const SmallArray<const Rd>& symmetry_normal_list,
const CellToCellStencilArray& stencil_array);
- BoundaryIntegralReconstructionMatrixBuilder() = default;
~BoundaryIntegralReconstructionMatrixBuilder() = default;
};
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
index 198dba9afe1faf1c324a8c03ae68376dc27ac799..24579e0594276b7a1a992aa3d2e8bd09375e333a 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.cpp
@@ -2,6 +2,8 @@
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/GaussQuadratureDescriptor.hpp>
+#include <analysis/QuadratureFormula.hpp>
+#include <analysis/QuadratureManager.hpp>
#include <geometry/CubeTransformation.hpp>
#include <geometry/LineTransformation.hpp>
#include <geometry/PrismTransformation.hpp>
@@ -76,13 +78,7 @@ PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder<MeshTypeT>:
} else if constexpr (MeshType::Dimension == 3) {
static_assert(MeshType::Dimension == 3);
- const double detT = [&] {
- if constexpr (std::is_same_v<TetrahedronTransformation, std::decay_t<decltype(T)>>) {
- return T.jacobianDeterminant();
- } else {
- return T.jacobianDeterminant(xi_q);
- }
- }();
+ const double detT = T.jacobianDeterminant(xi_q);
const double x_xj = X_Xj[0];
const double y_yj = X_Xj[1];
diff --git a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
index 1f4573438260e321dbd1461d2566f854236b25f4..7a3e8f7fc04d6d4ba21d4a36ed02ba24dd0d67d0 100644
--- a/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
+++ b/src/scheme/reconstruction_utils/ElementIntegralReconstructionMatrixBuilder.hpp
@@ -3,14 +3,15 @@
#include <algebra/ShrinkMatrixView.hpp>
#include <algebra/SmallMatrix.hpp>
-#include <analysis/QuadratureFormula.hpp>
-#include <analysis/QuadratureManager.hpp>
#include <mesh/CellType.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/StencilArray.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <utils/SmallArray.hpp>
+template <size_t Dimension>
+class QuadratureFormula;
+
template <MeshConcept MeshTypeT>
class PolynomialReconstruction::ElementIntegralReconstructionMatrixBuilder
{
diff --git a/src/utils/SignalManager.cpp b/src/utils/SignalManager.cpp
index 4b2046a98796fc72cef406f77e7498930089f533..5054622c0d9959c089ddfa3f3a323820614c5806 100644
--- a/src/utils/SignalManager.cpp
+++ b/src/utils/SignalManager.cpp
@@ -52,17 +52,27 @@ SignalManager::signalName(int signal)
void
SignalManager::pauseForDebug(int signal)
{
- if (std::string(PUGS_BUILD_TYPE) != "Release") {
- if (s_pause_on_error) {
- // Each failing process must write
- std::cerr.clear();
+ if (s_pause_on_error) {
+ // Each failing process must write
+ std::cerr.clear();
+ if (std::string(PUGS_BUILD_TYPE) != "Release") {
+ char hostname[HOST_NAME_MAX + 1];
+ gethostname(hostname, HOST_NAME_MAX + 1);
std::cerr << "\n======================================\n"
- << rang::style::reset << rang::fg::reset << rang::style::bold << "to attach gdb to this process run\n"
+ << rang::style::reset << rang::fg::reset << rang::style::bold << "Process paused on host \""
+ << rang::fg::yellow << hostname << rang::fg::reset << "\"\n"
+ << "to attach gdb to this process run\n"
<< "\tgdb -pid " << rang::fg::red << getpid() << rang::fg::reset << '\n'
<< "else press Control-C to exit\n"
<< rang::style::reset << "======================================\n"
<< std::flush;
pause();
+ } else {
+ std::cerr << '\n'
+ << rang::style::bold
+ << "Pausing is useless for Release version.\n"
+ "To attach debugger use Debug built type."
+ << rang::style::reset << '\n';
}
}
std::exit(signal);
@@ -73,56 +83,52 @@ SignalManager::handler(int signal)
{
static std::mutex mutex;
- if (mutex.try_lock()) {
- std::signal(SIGTERM, SIG_DFL);
- std::signal(SIGINT, SIG_DFL);
- std::signal(SIGABRT, SIG_DFL);
-
- // Each failing process must write
- std::cerr.clear();
+ std::lock_guard<std::mutex> lock(mutex);
- std::cerr << BacktraceManager{} << '\n';
+ std::signal(SIGINT, SIG_BLOCK);
- std::cerr << "\n *** " << rang::style::reset << rang::fg::reset << rang::style::bold << "Signal " << rang::fgB::red
- << signalName(signal) << rang::fg::reset << " caught" << rang::style::reset << " ***\n\n";
+ // Each failing process must write
+ std::cerr.clear();
- std::exception_ptr eptr = std::current_exception();
- try {
- if (eptr) {
- std::rethrow_exception(eptr);
- } else {
- std::ostringstream error_msg;
- error_msg << "received " << signalName(signal);
- std::cerr << ASTExecutionStack::getInstance().errorMessageAt(error_msg.str()) << '\n';
- }
- }
- catch (const IBacktraceError& backtrace_error) {
- auto source_location = backtrace_error.sourceLocation();
- std::cerr << rang::fgB::cyan << source_location.file_name() << ':' << source_location.line() << ':'
- << source_location.column() << ':' << rang::fg::reset << rang::fgB::yellow
- << " threw the following exception" << rang::fg::reset << "\n\n";
- std::cerr << ASTExecutionStack::getInstance().errorMessageAt(backtrace_error.what()) << '\n';
- }
- catch (const ParseError& parse_error) {
- auto p = parse_error.positions().front();
- std::cerr << rang::style::bold << p.source << ':' << p.line << ':' << p.column << ':' << rang::style::reset
- << rang::fgB::red << " error: " << rang::fg::reset << rang::style::bold << parse_error.what()
- << rang::style::reset << '\n';
- }
- catch (const IExitError& exit_error) {
- std::cerr << ASTExecutionStack::getInstance().errorMessageAt(exit_error.what()) << '\n';
- }
- catch (const AssertError& assert_error) {
- std::cerr << assert_error << '\n';
- }
- catch (...) {
- std::cerr << "Unknown exception!\n";
- }
+ std::cerr << BacktraceManager{} << '\n';
- SignalManager::pauseForDebug(signal);
+ std::cerr << "\n *** " << rang::style::reset << rang::fg::reset << rang::style::bold << "Signal " << rang::fgB::red
+ << signalName(signal) << rang::fg::reset << " caught" << rang::style::reset << " ***\n\n";
- mutex.unlock();
+ std::exception_ptr eptr = std::current_exception();
+ try {
+ if (eptr) {
+ std::rethrow_exception(eptr);
+ } else {
+ std::ostringstream error_msg;
+ error_msg << "received " << signalName(signal);
+ std::cerr << ASTExecutionStack::getInstance().errorMessageAt(error_msg.str()) << '\n';
+ }
+ }
+ catch (const IBacktraceError& backtrace_error) {
+ auto source_location = backtrace_error.sourceLocation();
+ std::cerr << rang::fgB::cyan << source_location.file_name() << ':' << source_location.line() << ':'
+ << source_location.column() << ':' << rang::fg::reset << rang::fgB::yellow
+ << " threw the following exception" << rang::fg::reset << "\n\n";
+ std::cerr << ASTExecutionStack::getInstance().errorMessageAt(backtrace_error.what()) << '\n';
}
+ catch (const ParseError& parse_error) {
+ auto p = parse_error.positions().front();
+ std::cerr << rang::style::bold << p.source << ':' << p.line << ':' << p.column << ':' << rang::style::reset
+ << rang::fgB::red << " error: " << rang::fg::reset << rang::style::bold << parse_error.what()
+ << rang::style::reset << '\n';
+ }
+ catch (const IExitError& exit_error) {
+ std::cerr << ASTExecutionStack::getInstance().errorMessageAt(exit_error.what()) << '\n';
+ }
+ catch (const AssertError& assert_error) {
+ std::cerr << assert_error << '\n';
+ }
+ catch (...) {
+ std::cerr << "Unknown exception!\n";
+ }
+
+ SignalManager::pauseForDebug(signal);
}
void
diff --git a/tests/test_LinearSolver.cpp b/tests/test_LinearSolver.cpp
index 76ef6aa2bb75c3ea94fcf392f59ddb74350dc50a..286fe11ed4040486a4387105065adb955077baf8 100644
--- a/tests/test_LinearSolver.cpp
+++ b/tests/test_LinearSolver.cpp
@@ -837,77 +837,52 @@ TEST_CASE("LinearSolver", "[algebra]")
SECTION("Dense Matrices")
{
- SECTION("check linear solvers")
+ SECTION("SmallMatrix")
{
- SECTION("symmetric system")
+ SECTION("check linear solvers")
{
- SmallMatrix<double> A{5};
- A = zero;
-
- A(0, 0) = 2;
- A(0, 1) = -1;
-
- A(1, 0) = -1;
- A(1, 1) = 2;
- A(1, 2) = -1;
-
- A(2, 1) = -1;
- A(2, 2) = 2;
- A(2, 3) = -1;
-
- A(3, 2) = -1;
- A(3, 3) = 2;
- A(3, 4) = -1;
-
- A(4, 3) = -1;
- A(4, 4) = 2;
+ SECTION("symmetric system")
+ {
+ SmallMatrix<double> A{5};
+ A = zero;
- SmallVector<const double> x_exact = [] {
- SmallVector<double> y{5};
- y[0] = 1;
- y[1] = 3;
- y[2] = 2;
- y[3] = 4;
- y[4] = 5;
- return y;
- }();
+ A(0, 0) = 2;
+ A(0, 1) = -1;
- SmallVector<double> b = A * x_exact;
+ A(1, 0) = -1;
+ A(1, 1) = 2;
+ A(1, 2) = -1;
- SECTION("builtin")
- {
- SECTION("CG no preconditioner")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::builtin;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
+ A(2, 1) = -1;
+ A(2, 2) = 2;
+ A(2, 3) = -1;
- SmallVector<double> x{5};
- x = zero;
+ A(3, 2) = -1;
+ A(3, 3) = 2;
+ A(3, 4) = -1;
- LinearSolver solver{options};
+ A(4, 3) = -1;
+ A(4, 4) = 2;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
- }
+ SmallVector<const double> x_exact = [] {
+ SmallVector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
- SECTION("Eigen3")
- {
-#ifdef PUGS_HAS_EIGEN3
+ SmallVector<double> b = A * x_exact;
- SECTION("CG")
+ SECTION("builtin")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
- options.verbose() = true;
-
SECTION("CG no preconditioner")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
SmallVector<double> x{5};
@@ -919,89 +894,176 @@ TEST_CASE("LinearSolver", "[algebra]")
SmallVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
+ }
- SECTION("CG Diagonal")
+ SECTION("Eigen3")
+ {
+#ifdef PUGS_HAS_EIGEN3
+
+ SECTION("CG")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("CG no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ SmallVector<double> x{5};
+ x = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("CG Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("CG ICholesky")
+ {
+ options.precond() = LSPrecond::incomplete_cholesky;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+ "error: incomplete cholesky is not available for dense matrices in Eigen3");
+ }
+
+ SECTION("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
+ }
}
- SECTION("CG ICholesky")
+ SECTION("Cholesky")
{
- options.precond() = LSPrecond::incomplete_cholesky;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::cholesky;
+ options.precond() = LSPrecond::none;
SmallVector<double> x{5};
x = zero;
LinearSolver solver{options};
- REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
- "error: incomplete cholesky is not available for dense matrices in Eigen3");
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("CG AMG")
+#else // PUGS_HAS_EIGEN3
+ SECTION("Eigen3 not linked")
{
- options.precond() = LSPrecond::amg;
-
- SmallVector<double> x{5};
- x = zero;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
}
+#endif // PUGS_HAS_EIGEN3
}
- SECTION("Cholesky")
+ SECTION("PETSc")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::cholesky;
- options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
- SmallVector<double> x{5};
- x = zero;
+ SECTION("CG")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- LinearSolver solver{options};
+ SECTION("CG no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
+ SmallVector<double> x{5};
+ x = zero;
-#else // PUGS_HAS_EIGEN3
- SECTION("Eigen3 not linked")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
+ LinearSolver solver{options};
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
- }
-#endif // PUGS_HAS_EIGEN3
- }
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
- SECTION("PETSc")
- {
-#ifdef PUGS_HAS_PETSC
+ SECTION("CG Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
- SECTION("CG")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
- options.verbose() = true;
+ SmallVector<double> x{5};
+ x = zero;
- SECTION("CG no preconditioner")
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("CG ICholesky")
+ {
+ options.precond() = LSPrecond::incomplete_cholesky;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+ }
+
+ SECTION("Cholesky")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cholesky;
options.precond() = LSPrecond::none;
SmallVector<double> x{5};
@@ -1014,133 +1076,448 @@ TEST_CASE("LinearSolver", "[algebra]")
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("CG Diagonal")
+#else // PUGS_HAS_PETSC
+ SECTION("PETSc not linked")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_PETSC
+ }
+ }
- LinearSolver solver{options};
+ SECTION("none symmetric system")
+ {
+ SECTION("Dense matrix")
+ {
+ SmallMatrix<double> A{5};
+ A = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
+ A(0, 0) = 2;
+ A(0, 1) = -1;
- SECTION("CG ICholesky")
+ A(1, 0) = -0.2;
+ A(1, 1) = 2;
+ A(1, 2) = -1;
+
+ A(2, 1) = -1;
+ A(2, 2) = 4;
+ A(2, 3) = -2;
+
+ A(3, 2) = -1;
+ A(3, 3) = 2;
+ A(3, 4) = -0.1;
+
+ A(4, 3) = 1;
+ A(4, 4) = 3;
+
+ SmallVector<const double> x_exact = [] {
+ SmallVector<double> y{5};
+ y[0] = 1;
+ y[1] = 3;
+ y[2] = 2;
+ y[3] = 4;
+ y[4] = 5;
+ return y;
+ }();
+
+ SmallVector<double> b = A * x_exact;
+
+ SECTION("builtin")
{
- options.precond() = LSPrecond::incomplete_cholesky;
+ SECTION("BICGStab no preconditioner")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ SmallVector<double> x{5};
+ x = zero;
- LinearSolver solver{options};
+ LinearSolver solver{options};
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
}
- SECTION("CG AMG")
+ SECTION("Eigen3")
{
- options.precond() = LSPrecond::amg;
+#ifdef PUGS_HAS_EIGEN3
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- LinearSolver solver{options};
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+ "error: incomplete LU is not available for dense matrices in Eigen3");
+ }
+ }
+
+ SECTION("BICGStab2")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
+
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
+ }
+ }
+
+ SECTION("GMRES")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::gmres;
+ options.precond() = LSPrecond::none;
+
+ SECTION("GMRES no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+ "error: incomplete LU is not available for dense matrices in Eigen3");
+ }
+ }
+
+ SECTION("LU")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+#else // PUGS_HAS_EIGEN3
+ SECTION("Eigen3 not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_EIGEN3
}
- }
- SECTION("Cholesky")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cholesky;
- options.precond() = LSPrecond::none;
+ SECTION("PETSc")
+ {
+#ifdef PUGS_HAS_PETSC
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- LinearSolver solver{options};
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+ }
+
+ SECTION("BICGStab2")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
+
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab2 Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+ }
+
+ SECTION("GMRES")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::gmres;
+ options.precond() = LSPrecond::none;
+
+ SECTION("GMRES no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+ }
+
+ SECTION("LU")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
+
+ SmallVector<double> x{5};
+ x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ SmallVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
#else // PUGS_HAS_PETSC
- SECTION("PETSc not linked")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
+ SECTION("PETSc not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
- }
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
#endif // PUGS_HAS_PETSC
+ }
+ }
}
}
+ }
- SECTION("none symmetric system")
+ SECTION("TinyMatrix")
+ {
+ SECTION("check linear solvers")
{
- SECTION("Dense matrix")
+ SECTION("symmetric system")
{
- SmallMatrix<double> A{5};
- A = zero;
+ TinyMatrix<5> A = zero;
A(0, 0) = 2;
A(0, 1) = -1;
- A(1, 0) = -0.2;
+ A(1, 0) = -1;
A(1, 1) = 2;
A(1, 2) = -1;
A(2, 1) = -1;
- A(2, 2) = 4;
- A(2, 3) = -2;
+ A(2, 2) = 2;
+ A(2, 3) = -1;
A(3, 2) = -1;
A(3, 3) = 2;
- A(3, 4) = -0.1;
+ A(3, 4) = -1;
- A(4, 3) = 1;
- A(4, 4) = 3;
+ A(4, 3) = -1;
+ A(4, 4) = 2;
- SmallVector<const double> x_exact = [] {
- SmallVector<double> y{5};
- y[0] = 1;
- y[1] = 3;
- y[2] = 2;
- y[3] = 4;
- y[4] = 5;
- return y;
- }();
+ TinyVector<5> x_exact{1, 3, 2, 4, 5};
- SmallVector<double> b = A * x_exact;
+ TinyVector b = A * x_exact;
SECTION("builtin")
{
- SECTION("BICGStab no preconditioner")
+ SECTION("CG no preconditioner")
{
LinearSolverOptions options;
options.library() = LSLibrary::builtin;
- options.method() = LSMethod::bicgstab;
+ options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
@@ -1149,322 +1526,535 @@ TEST_CASE("LinearSolver", "[algebra]")
{
#ifdef PUGS_HAS_EIGEN3
- SECTION("BICGStab")
+ SECTION("CG")
{
LinearSolverOptions options;
options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::bicgstab;
+ options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
options.verbose() = true;
- SECTION("BICGStab no preconditioner")
+ SECTION("CG no preconditioner")
{
options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("BICGStab Diagonal")
+ SECTION("CG Diagonal")
{
options.precond() = LSPrecond::diagonal;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("BICGStab ILU")
+ SECTION("CG ICholesky")
{
- options.precond() = LSPrecond::incomplete_LU;
+ options.precond() = LSPrecond::incomplete_cholesky;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
- "error: incomplete LU is not available for dense matrices in Eigen3");
+ "error: incomplete cholesky is not available for dense matrices in Eigen3");
+ }
+
+ SECTION("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: AMG is not an Eigen3 preconditioner!");
}
}
- SECTION("BICGStab2")
+ SECTION("Cholesky")
{
LinearSolverOptions options;
options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::bicgstab2;
+ options.method() = LSMethod::cholesky;
options.precond() = LSPrecond::none;
- SECTION("BICGStab2 no preconditioner")
- {
- options.precond() = LSPrecond::none;
+ TinyVector<5> x = zero;
- SmallVector<double> x{5};
- x = zero;
+ LinearSolver solver{options};
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
- }
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("GMRES")
+#else // PUGS_HAS_EIGEN3
+ SECTION("Eigen3 not linked")
{
LinearSolverOptions options;
options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::gmres;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_EIGEN3
+ }
+
+ SECTION("PETSc")
+ {
+#ifdef PUGS_HAS_PETSC
+
+ SECTION("CG")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
+ options.verbose() = true;
- SECTION("GMRES no preconditioner")
+ SECTION("CG no preconditioner")
{
options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("GMRES Diagonal")
+ SECTION("CG Diagonal")
{
options.precond() = LSPrecond::diagonal;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
- SECTION("GMRES ILU")
+ SECTION("CG ICholesky")
{
- options.precond() = LSPrecond::incomplete_LU;
+ options.precond() = LSPrecond::incomplete_cholesky;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
- REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
- "error: incomplete LU is not available for dense matrices in Eigen3");
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("CG AMG")
+ {
+ options.precond() = LSPrecond::amg;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
}
- SECTION("LU")
+ SECTION("Cholesky")
{
LinearSolverOptions options;
- options.library() = LSLibrary::eigen3;
- options.method() = LSMethod::lu;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cholesky;
options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
-#else // PUGS_HAS_EIGEN3
- SECTION("Eigen3 not linked")
+#else // PUGS_HAS_PETSC
+ SECTION("PETSc not linked")
{
LinearSolverOptions options;
- options.library() = LSLibrary::eigen3;
+ options.library() = LSLibrary::petsc;
options.method() = LSMethod::cg;
options.precond() = LSPrecond::none;
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
}
-#endif // PUGS_HAS_EIGEN3
+#endif // PUGS_HAS_PETSC
}
+ }
- SECTION("PETSc")
+ SECTION("none symmetric system")
+ {
+ SECTION("Dense matrix")
{
-#ifdef PUGS_HAS_PETSC
+ TinyMatrix<5> A = zero;
- SECTION("BICGStab")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab;
- options.precond() = LSPrecond::none;
- options.verbose() = true;
+ A(0, 0) = 2;
+ A(0, 1) = -1;
+
+ A(1, 0) = -0.2;
+ A(1, 1) = 2;
+ A(1, 2) = -1;
+
+ A(2, 1) = -1;
+ A(2, 2) = 4;
+ A(2, 3) = -2;
+
+ A(3, 2) = -1;
+ A(3, 3) = 2;
+ A(3, 4) = -0.1;
+
+ A(4, 3) = 1;
+ A(4, 4) = 3;
+
+ const TinyVector<5> x_exact{1, 3, 2, 4, 5
+
+ };
+
+ const TinyVector b = A * x_exact;
+ SECTION("builtin")
+ {
SECTION("BICGStab no preconditioner")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::builtin;
+ options.method() = LSMethod::bicgstab;
options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
+ }
+
+ SECTION("Eigen3")
+ {
+#ifdef PUGS_HAS_EIGEN3
- SECTION("BICGStab Diagonal")
+ SECTION("BICGStab")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::bicgstab;
+ options.precond() = LSPrecond::none;
+ options.verbose() = true;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ TinyVector<5> x = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
+ LinearSolver solver{options};
- SECTION("BICGStab ILU")
- {
- options.precond() = LSPrecond::incomplete_LU;
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
- LinearSolver solver{options};
+ TinyVector<5> x = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+ "error: incomplete LU is not available for dense matrices in Eigen3");
+ }
}
- }
- SECTION("BICGStab2")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::bicgstab2;
- options.precond() = LSPrecond::none;
+ SECTION("BICGStab2")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
- SECTION("BICGStab2 no preconditioner")
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: BICGStab2 is not an Eigen3 linear solver!");
+ }
+ }
+
+ SECTION("GMRES")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::gmres;
options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("GMRES no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ TinyVector<5> x = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ REQUIRE_THROWS_WITH(solver.solveLocalSystem(A, x, b),
+ "error: incomplete LU is not available for dense matrices in Eigen3");
+ }
}
- SECTION("BICGStab2 Diagonal")
+ SECTION("LU")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
LinearSolver solver{options};
solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
+ TinyVector error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
}
+
+#else // PUGS_HAS_EIGEN3
+ SECTION("Eigen3 not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::eigen3;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
+
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: Eigen3 is not linked to pugs. Cannot use it!");
+ }
+#endif // PUGS_HAS_EIGEN3
}
- SECTION("GMRES")
+ SECTION("PETSc")
{
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::gmres;
- options.precond() = LSPrecond::none;
+#ifdef PUGS_HAS_PETSC
- SECTION("GMRES no preconditioner")
+ SECTION("BICGStab")
{
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab;
options.precond() = LSPrecond::none;
+ options.verbose() = true;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ TinyVector<5> x = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
}
- SECTION("GMRES Diagonal")
+ SECTION("BICGStab2")
{
- options.precond() = LSPrecond::diagonal;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::bicgstab2;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("BICGStab2 no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ TinyVector<5> x = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("BICGStab2 Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
}
- SECTION("GMRES ILU")
+ SECTION("GMRES")
{
- options.precond() = LSPrecond::incomplete_LU;
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::gmres;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ SECTION("GMRES no preconditioner")
+ {
+ options.precond() = LSPrecond::none;
- LinearSolver solver{options};
+ TinyVector<5> x = zero;
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES Diagonal")
+ {
+ options.precond() = LSPrecond::diagonal;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
+
+ SECTION("GMRES ILU")
+ {
+ options.precond() = LSPrecond::incomplete_LU;
+
+ TinyVector<5> x = zero;
+
+ LinearSolver solver{options};
+
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
}
- }
- SECTION("LU")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::lu;
- options.precond() = LSPrecond::none;
+ SECTION("LU")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::lu;
+ options.precond() = LSPrecond::none;
- SmallVector<double> x{5};
- x = zero;
+ TinyVector<5> x = zero;
- LinearSolver solver{options};
+ LinearSolver solver{options};
- solver.solveLocalSystem(A, x, b);
- SmallVector error = x - x_exact;
- REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
- }
+ solver.solveLocalSystem(A, x, b);
+ TinyVector error = x - x_exact;
+ REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x_exact, x_exact)));
+ }
#else // PUGS_HAS_PETSC
- SECTION("PETSc not linked")
- {
- LinearSolverOptions options;
- options.library() = LSLibrary::petsc;
- options.method() = LSMethod::cg;
- options.precond() = LSPrecond::none;
+ SECTION("PETSc not linked")
+ {
+ LinearSolverOptions options;
+ options.library() = LSLibrary::petsc;
+ options.method() = LSMethod::cg;
+ options.precond() = LSPrecond::none;
- REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
- }
+ REQUIRE_THROWS_WITH(LinearSolver{options}, "error: PETSc is not linked to pugs. Cannot use it!");
+ }
#endif // PUGS_HAS_PETSC
+ }
}
}
}
diff --git a/tests/test_PolynomialReconstruction_degree_2.cpp b/tests/test_PolynomialReconstruction_degree_2.cpp
index a1c5c058975c17f9a5ec7e40c41f28b9abb18a66..db13612c1d8c5efad94a6f890251e366826d014b 100644
--- a/tests/test_PolynomialReconstruction_degree_2.cpp
+++ b/tests/test_PolynomialReconstruction_degree_2.cpp
@@ -118,7 +118,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
auto dpk_Ah = reconstructions[0]->get<DiscreteFunctionDPk<1, const R3x3>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Ah, std::function(R3x3_exact));
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -139,7 +139,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const double>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -172,7 +172,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
auto dpk_Vh = reconstructions[0]->get<DiscreteFunctionDPkVector<1, const R3>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
@@ -219,7 +219,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
auto dpk_uh = reconstructions[1]->get<DiscreteFunctionDPk<1, const R3>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_uh, std::function(R3_exact));
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
{
@@ -231,7 +231,7 @@ TEST_CASE("PolynomialReconstruction_degree_2", "[scheme]")
{
auto dpk_Vh = reconstructions[3]->get<DiscreteFunctionDPkVector<1, const double>>();
double max_error = test_only::max_reconstruction_error(mesh, dpk_Vh, vector_exact);
- REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-14));
+ REQUIRE(parallel::allReduceMax(max_error) == Catch::Approx(0).margin(1E-12));
}
}
}
diff --git a/tools/plugin-template/CMakeLists.txt-template b/tools/plugin-template/CMakeLists.txt-template
index 9c6478b0711cbbd5e38a41f3eaa6c4736a300e67..369887ba7ffd5dfeb9b3593157c51d7a1beefd2f 100644
--- a/tools/plugin-template/CMakeLists.txt-template
+++ b/tools/plugin-template/CMakeLists.txt-template
@@ -23,14 +23,21 @@ find_package(Pugs REQUIRED)
list(APPEND CMAKE_MODULE_PATH "${PUGS_PREFIX_PATH}/lib/cmake/Kokkos")
include(KokkosConfig)
-set(HDF5_PREFER_PARALLEL TRUE)
-list(APPEND CMAKE_MODULE_PATH "${PUGS_PREFIX_PATH}/lib/cmake/HighFive")
-include(HighFiveConfig)
-
list(APPEND CMAKE_MODULE_PATH "${PUGS_PREFIX_PATH}/lib/cmake/pugs")
-include(PugsTargets)
include(PugsCompileFlags)
+if (${PUGS_HAS_MPI})
+ set(HDF5_PREFER_PARALLEL TRUE)
+ list(APPEND CMAKE_MODULE_PATH "${PUGS_PREFIX_PATH}/lib/cmake/HighFive")
+
+ include(HighFiveConfig)
+ if (TARGET HighFive)
+ set(HIGHFIVE_TARGET HighFive::HighFive)
+ endif()
+endif()
+
+include(PugsTargets)
+
#------------------------------------------------------
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${PUGS_CMAKE_CXX_FLAGS}")
@@ -77,16 +84,11 @@ include_directories("${CMAKE_CURRENT_SOURCE_DIR}")
include_directories(SYSTEM "${PUGS_PREFIX_PATH}/include")
include_directories(SYSTEM "${PUGS_PREFIX_PATH}/include/kokkos")
include_directories(SYSTEM "${PUGS_PREFIX_PATH}/include/tao/")
+if (${PUGS_HAS_MPI})
include_directories(SYSTEM "${MPI_CXX_INCLUDE_DIRS}")
+endif()
-get_target_property(_prop Pugs::PugsAlgebra INTERFACE_INCLUDE_DIRECTORIES)
-set(PUGS_INC_DIR "${PUGS_INC_DIR};${_prop}")
-get_target_property(_prop Pugs::PugsUtils INTERFACE_INCLUDE_DIRECTORIES)
-set(PUGS_INC_DIR "${PUGS_INC_DIR};${_prop}")
-get_target_property(_prop Pugs::pugs INTERFACE_INCLUDE_DIRECTORIES)
-set(PUGS_INC_DIR "${PUGS_INC_DIR};${_prop}")
-
-include_directories(SYSTEM ${PUGS_INC_DIR})
+include_directories(SYSTEM "${PUGS_PREFIX_PATH}/include")
link_directories(${PUGS_PREFIX_PATH}/lib)
#------------------------------------------------------
@@ -143,6 +145,10 @@ add_library(_PLUGIN_NAME_
# add cpp sources files here
)
+target_link_libraries(_PLUGIN_NAME_
+ ${HIGHFIVE_TARGET}
+)
+
#------------------------------------------------------
add_subdirectory(tests)