diff --git a/src/algebra/BiCGStab.hpp b/src/algebra/BiCGStab.hpp
index c931385b214c91ca2807db4f07fca027ebd79b25..bc822668a7f63c4d5d7cc4c45f4d1b1d70084fc0 100644
--- a/src/algebra/BiCGStab.hpp
+++ b/src/algebra/BiCGStab.hpp
@@ -5,14 +5,19 @@
#include <iomanip>
#include <iostream>
+#include <rang.hpp>
+
+template <bool verbose = true>
struct BiCGStab
{
template <typename VectorType, typename MatrixType>
BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6)
{
- std::cout << "- bi-conjugate gradient stabilized\n";
- std::cout << " epsilon = " << epsilon << '\n';
- std::cout << " maximum number of iterations: " << max_iter << '\n';
+ if constexpr (verbose) {
+ std::cout << "- bi-conjugate gradient stabilized\n";
+ std::cout << " epsilon = " << epsilon << '\n';
+ std::cout << " maximum number of iterations: " << max_iter << '\n';
+ }
VectorType r_k_1{b.size()};
@@ -33,10 +38,14 @@ struct BiCGStab
VectorType r_k{x.size()};
- std::cout << " initial residu: " << resid0 << '\n';
+ if constexpr (verbose) {
+ std::cout << " initial residu: " << resid0 << '\n';
+ }
for (size_t i = 1; i <= max_iter; ++i) {
- std::cout << " - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
- << "\tabsolute: " << residu << '\n';
+ if constexpr (verbose) {
+ std::cout << " - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
+ << "\tabsolute: " << residu << '\n';
+ }
Ap_k = A * p_k;
@@ -64,14 +73,14 @@ struct BiCGStab
}
if (residu / resid0 > epsilon) {
- std::cout << " bi_conjugate gradient stabilized: *NOT CONVERGED*\n";
+ std::cout << " bi-conjugate gradient stabilized: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset
+ << '\n';
+ ;
std::cout << " - epsilon: " << epsilon << '\n';
std::cout << " - relative residu : " << residu / resid0 << '\n';
std::cout << " - absolute residu : " << residu << '\n';
}
}
-
- std::cout << "- bi-conjugate gradient stabilized: done\n";
}
};
diff --git a/src/algebra/PCG.hpp b/src/algebra/PCG.hpp
index 74c17ee40971b96c512c196d630f1edfadfa6ddc..36a4635b313d2bba13d97e1b7a4ea1882e2a66a2 100644
--- a/src/algebra/PCG.hpp
+++ b/src/algebra/PCG.hpp
@@ -4,6 +4,9 @@
#include <iomanip>
#include <iostream>
+#include <rang.hpp>
+
+template <bool verbose = true>
struct PCG
{
template <typename VectorType, typename MatrixType, typename PreconditionerType>
@@ -14,14 +17,16 @@ struct PCG
const size_t maxiter,
const double epsilon = 1e-6)
{
- std::cout << "- conjugate gradient\n";
- std::cout << " epsilon = " << epsilon << '\n';
- std::cout << " maximum number of iterations: " << maxiter << '\n';
+ if constexpr (verbose) {
+ std::cout << "- conjugate gradient\n";
+ std::cout << " epsilon = " << epsilon << '\n';
+ std::cout << " maximum number of iterations: " << maxiter << '\n';
+ }
VectorType h{f.size()};
VectorType b = copy(f);
- if (true) {
+ if constexpr (verbose) {
h = A * x;
h -= f;
std::cout << "- initial *real* residu : " << (h, h) << '\n';
@@ -69,10 +74,14 @@ struct PCG
if ((i == 1) && (gcg != 0)) {
relativeEpsilon = epsilon * gcg;
gcg0 = gcg;
- std::cout << " initial residu: " << gcg << '\n';
+ if constexpr (verbose) {
+ std::cout << " initial residu: " << gcg << '\n';
+ }
+ }
+ if constexpr (verbose) {
+ std::cout << " - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
+ std::cout << "\tabsolute: " << gcg << '\n';
}
- std::cout << " - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
- std::cout << "\tabsolute: " << gcg << '\n';
if (gcg < relativeEpsilon) {
break;
@@ -85,7 +94,7 @@ struct PCG
}
if (gcg > relativeEpsilon) {
- std::cout << " conjugate gradient: *NOT CONVERGED*\n";
+ std::cout << " conjugate gradient: " << rang::fgB::red << "*NOT CONVERGED*" << rang::style::reset << '\n';
std::cout << " - epsilon: " << epsilon << '\n';
std::cout << " - relative residu : " << gcg / gcg0 << '\n';
std::cout << " - absolute residu : " << gcg << '\n';
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index df4246a0ee36f2277e50d34398bf6de7a8d5ad3f..c82e9ce06ec8b7d43ce1c6cc2fec282f89dc1fee 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -155,7 +155,7 @@ class MeshData : public IMeshData
PUGS_INLINE
void
_computeCellCentroid()
- { // Computes vertices isobarycenter
+ {
const CellValue<const Rd> cell_iso_barycenter = this->cellIsoBarycenter();
if constexpr (Dimension == 1) {
m_cell_centroid = cell_iso_barycenter;
@@ -190,8 +190,53 @@ class MeshData : public IMeshData
});
m_cell_centroid = cell_centroid;
} else {
- std::cout << rang::fgB::yellow << "Centroids are not computed correctly in 3d" << rang::style::reset << '\n';
- m_cell_centroid = cell_iso_barycenter;
+ const auto& face_center = this->xl();
+ const CellValue<const double> Vj = this->Vj();
+ const NodeValue<const Rd>& xr = m_mesh.xr();
+ const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+ const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+ const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+
+ CellValue<Rd> cell_centroid{m_mesh.connectivity()};
+ parallel_for(
+ m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+ const Rd xj = m_cell_iso_barycenter[j];
+
+ const auto& cell_faces = cell_to_face_matrix[j];
+
+ Rd sum = zero;
+ for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+ const FaceId face_id = cell_faces[i_face];
+
+ const Rd xl = face_center[face_id];
+
+ const Rd xjxl = xl - xj;
+
+ const auto& face_nodes = face_to_node_matrix[face_id];
+
+ const Rd xl_plus_xj = xl + xj;
+
+ double sign = (cell_face_is_reversed(j, i_face)) ? -1 : 1;
+
+ for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
+ const Rd& xr0 = xr[face_nodes[i_face_node]];
+ const Rd& xr1 = xr[face_nodes[(i_face_node + 1) % face_nodes.size()]];
+
+ const Rd xjxr0 = xr0 - xj;
+ const Rd xjxr1 = xr1 - xj;
+
+ const double Vjlr = (crossProduct(xjxr0, xjxr1), xjxl);
+
+ sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1);
+ }
+ }
+
+ sum *= 1 / (24 * Vj[j]);
+
+ cell_centroid[j] = sum;
+ });
+
+ m_cell_centroid = cell_centroid;
}
}
}
diff --git a/tests/test_BiCGStab.cpp b/tests/test_BiCGStab.cpp
index b7cb2d3bd8afe280a7d8c416a662016851e0d820..99b94515ba7795903434db761d509b9052fcaf71 100644
--- a/tests/test_BiCGStab.cpp
+++ b/tests/test_BiCGStab.cpp
@@ -45,7 +45,7 @@ TEST_CASE("BiCGStab", "[algebra]")
Vector<double> x{5};
x = 0;
- BiCGStab{b, A, x, 10, 1e-12};
+ BiCGStab<false>{b, A, x, 10, 1e-12};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
}
diff --git a/tests/test_PCG.cpp b/tests/test_PCG.cpp
index 92c0b4caf10824c25a1bef352dcd86ca5151b1e5..acc48b91e7f362e2f24f830465a9e7fc7f53f6d3 100644
--- a/tests/test_PCG.cpp
+++ b/tests/test_PCG.cpp
@@ -64,7 +64,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5};
x = 0;
- PCG{b, A, A, x, 10, 1e-12};
+ PCG<false>{b, A, A, x, 10, 1e-12};
REQUIRE(std::sqrt((x, x)) == 0);
}
@@ -106,7 +106,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5};
x = 0;
- PCG{b, A, A, x, 1, 1e-12};
+ PCG<false>{b, A, A, x, 1, 1e-12};
Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
}