Skip to content
Snippets Groups Projects
Commit b004a1c2 authored by PATELA Julie's avatar PATELA Julie
Browse files

Merge branch 'develop' into feature/Navier-Stokes

parents 97395ce9 8ff64fab
Branches
No related tags found
No related merge requests found
...@@ -5,14 +5,19 @@ ...@@ -5,14 +5,19 @@
#include <iomanip> #include <iomanip>
#include <iostream> #include <iostream>
#include <rang.hpp>
template <bool verbose = true>
struct BiCGStab struct BiCGStab
{ {
template <typename VectorType, typename MatrixType> template <typename VectorType, typename MatrixType>
BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6) BiCGStab(const VectorType& b, const MatrixType& A, VectorType& x, const size_t max_iter, const double epsilon = 1e-6)
{ {
if constexpr (verbose) {
std::cout << "- bi-conjugate gradient stabilized\n"; std::cout << "- bi-conjugate gradient stabilized\n";
std::cout << " epsilon = " << epsilon << '\n'; std::cout << " epsilon = " << epsilon << '\n';
std::cout << " maximum number of iterations: " << max_iter << '\n'; std::cout << " maximum number of iterations: " << max_iter << '\n';
}
VectorType r_k_1{b.size()}; VectorType r_k_1{b.size()};
...@@ -33,10 +38,14 @@ struct BiCGStab ...@@ -33,10 +38,14 @@ struct BiCGStab
VectorType r_k{x.size()}; VectorType r_k{x.size()};
if constexpr (verbose) {
std::cout << " initial residu: " << resid0 << '\n'; std::cout << " initial residu: " << resid0 << '\n';
}
for (size_t i = 1; i <= max_iter; ++i) { for (size_t i = 1; i <= max_iter; ++i) {
if constexpr (verbose) {
std::cout << " - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0 std::cout << " - iteration: " << std::setw(6) << i << "\tresidu: " << residu / resid0
<< "\tabsolute: " << residu << '\n'; << "\tabsolute: " << residu << '\n';
}
Ap_k = A * p_k; Ap_k = A * p_k;
...@@ -64,14 +73,14 @@ struct BiCGStab ...@@ -64,14 +73,14 @@ struct BiCGStab
} }
if (residu / resid0 > epsilon) { 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 << " - epsilon: " << epsilon << '\n';
std::cout << " - relative residu : " << residu / resid0 << '\n'; std::cout << " - relative residu : " << residu / resid0 << '\n';
std::cout << " - absolute residu : " << residu << '\n'; std::cout << " - absolute residu : " << residu << '\n';
} }
} }
std::cout << "- bi-conjugate gradient stabilized: done\n";
} }
}; };
......
...@@ -4,6 +4,9 @@ ...@@ -4,6 +4,9 @@
#include <iomanip> #include <iomanip>
#include <iostream> #include <iostream>
#include <rang.hpp>
template <bool verbose = true>
struct PCG struct PCG
{ {
template <typename VectorType, typename MatrixType, typename PreconditionerType> template <typename VectorType, typename MatrixType, typename PreconditionerType>
...@@ -14,14 +17,16 @@ struct PCG ...@@ -14,14 +17,16 @@ struct PCG
const size_t maxiter, const size_t maxiter,
const double epsilon = 1e-6) const double epsilon = 1e-6)
{ {
if constexpr (verbose) {
std::cout << "- conjugate gradient\n"; std::cout << "- conjugate gradient\n";
std::cout << " epsilon = " << epsilon << '\n'; std::cout << " epsilon = " << epsilon << '\n';
std::cout << " maximum number of iterations: " << maxiter << '\n'; std::cout << " maximum number of iterations: " << maxiter << '\n';
}
VectorType h{f.size()}; VectorType h{f.size()};
VectorType b = copy(f); VectorType b = copy(f);
if (true) { if constexpr (verbose) {
h = A * x; h = A * x;
h -= f; h -= f;
std::cout << "- initial *real* residu : " << (h, h) << '\n'; std::cout << "- initial *real* residu : " << (h, h) << '\n';
...@@ -69,10 +74,14 @@ struct PCG ...@@ -69,10 +74,14 @@ struct PCG
if ((i == 1) && (gcg != 0)) { if ((i == 1) && (gcg != 0)) {
relativeEpsilon = epsilon * gcg; relativeEpsilon = epsilon * gcg;
gcg0 = gcg; gcg0 = gcg;
if constexpr (verbose) {
std::cout << " initial residu: " << gcg << '\n'; std::cout << " initial residu: " << gcg << '\n';
} }
}
if constexpr (verbose) {
std::cout << " - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0; std::cout << " - iteration " << std::setw(6) << i << "\tresidu: " << gcg / gcg0;
std::cout << "\tabsolute: " << gcg << '\n'; std::cout << "\tabsolute: " << gcg << '\n';
}
if (gcg < relativeEpsilon) { if (gcg < relativeEpsilon) {
break; break;
...@@ -85,7 +94,7 @@ struct PCG ...@@ -85,7 +94,7 @@ struct PCG
} }
if (gcg > relativeEpsilon) { 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 << " - epsilon: " << epsilon << '\n';
std::cout << " - relative residu : " << gcg / gcg0 << '\n'; std::cout << " - relative residu : " << gcg / gcg0 << '\n';
std::cout << " - absolute residu : " << gcg << '\n'; std::cout << " - absolute residu : " << gcg << '\n';
......
...@@ -155,7 +155,7 @@ class MeshData : public IMeshData ...@@ -155,7 +155,7 @@ class MeshData : public IMeshData
PUGS_INLINE PUGS_INLINE
void void
_computeCellCentroid() _computeCellCentroid()
{ // Computes vertices isobarycenter {
const CellValue<const Rd> cell_iso_barycenter = this->cellIsoBarycenter(); const CellValue<const Rd> cell_iso_barycenter = this->cellIsoBarycenter();
if constexpr (Dimension == 1) { if constexpr (Dimension == 1) {
m_cell_centroid = cell_iso_barycenter; m_cell_centroid = cell_iso_barycenter;
...@@ -190,8 +190,53 @@ class MeshData : public IMeshData ...@@ -190,8 +190,53 @@ class MeshData : public IMeshData
}); });
m_cell_centroid = cell_centroid; m_cell_centroid = cell_centroid;
} else { } else {
std::cout << rang::fgB::yellow << "Centroids are not computed correctly in 3d" << rang::style::reset << '\n'; const auto& face_center = this->xl();
m_cell_centroid = cell_iso_barycenter; 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;
} }
} }
} }
......
...@@ -45,7 +45,7 @@ TEST_CASE("BiCGStab", "[algebra]") ...@@ -45,7 +45,7 @@ TEST_CASE("BiCGStab", "[algebra]")
Vector<double> x{5}; Vector<double> x{5};
x = 0; x = 0;
BiCGStab{b, A, x, 10, 1e-12}; BiCGStab<false>{b, A, x, 10, 1e-12};
Vector error = x - x_exact; Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x))); REQUIRE(std::sqrt((error, error)) < 1E-10 * std::sqrt((x, x)));
} }
......
...@@ -64,7 +64,7 @@ TEST_CASE("PCG", "[algebra]") ...@@ -64,7 +64,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5}; Vector<double> x{5};
x = 0; 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); REQUIRE(std::sqrt((x, x)) == 0);
} }
...@@ -106,7 +106,7 @@ TEST_CASE("PCG", "[algebra]") ...@@ -106,7 +106,7 @@ TEST_CASE("PCG", "[algebra]")
Vector<double> x{5}; Vector<double> x{5};
x = 0; x = 0;
PCG{b, A, A, x, 1, 1e-12}; PCG<false>{b, A, A, x, 1, 1e-12};
Vector error = x - x_exact; Vector error = x - x_exact;
REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x))); REQUIRE(std::sqrt((error, error)) > 1E-10 * std::sqrt((x, x)));
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment