Skip to content
Snippets Groups Projects
Commit d9245809 authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

Add polynomial reconstruction

Works in parallel for degree 1 polynomials in 1d
parent 83022e21
No related branches found
No related tags found
1 merge request!205High-order polynomial reconstruction
#ifndef GIVENS_HPP
#define GIVENS_HPP
#include <cmath>
#include <iomanip>
#include <iostream>
#include <rang.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
class Givens
{
private:
static void
_givens(const double a, const double b, double& c, double& s)
{
if (b == 0) {
c = 1;
s = 0;
} else {
if (std::abs(b) > std::abs(a)) {
const double tau = -a / b;
s = 1. / sqrt(1 + tau * tau);
c = s * tau;
} else {
Assert(a != 0);
const double tau = -b / a;
c = 1 / sqrt(1 + tau * tau);
s = c * tau;
}
}
}
static void
_rotate(const double Ai_1j, const double Aij, const double c, const double s, double& ui_1, double& ui)
{
ui_1 = c * Ai_1j - s * Aij;
ui = s * Ai_1j + c * Aij;
}
template <typename MatrixType, typename VectorType, typename RHSVectorType>
static void
_lift(const MatrixType& A, const RHSVectorType& b, VectorType& x)
{
x.fill(0);
for (size_t i = x.size(); i >= 1; --i) {
double sum = 0;
for (size_t j = i; j < x.size(); ++j) {
sum += A(i - 1, j) * x[j];
}
x[i - 1] = (b[i - 1] - sum) / A(i - 1, i - 1);
}
}
public:
template <typename MatrixType, typename VectorType, typename RHSVectorType>
static void
solveInPlace(MatrixType& A, VectorType& x, RHSVectorType& b)
{
Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
for (size_t j = 0; j < A.numberOfColumns(); ++j) {
for (size_t i = A.numberOfRows() - 1; i > j; --i) {
double c(0), s(0);
_givens(A(i - 1, j), A(i, j), c, s);
for (size_t k = j; k < A.numberOfColumns(); ++k) {
double xi_1(0), xi(0);
_rotate(A(i - 1, k), A(i, k), c, s, xi_1, xi);
A(i - 1, k) = xi_1;
A(i, k) = xi;
}
double xi_1(0), xi(0);
_rotate(b[i - 1], b[i], c, s, xi_1, xi);
b[i - 1] = xi_1;
b[i] = xi;
}
}
_lift(A, b, x);
}
template <typename MatrixType, typename VectorType, typename RHSVectorType>
static void
solve(const MatrixType& A, VectorType& x, const RHSVectorType& b)
{
Assert(x.size() == A.numberOfColumns(), "The number of columns of A must be the size of x");
Assert(b.size() == A.numberOfRows(), "The number of rows of A must be the size of b");
MatrixType rotateA = copy(A);
RHSVectorType rotateb = copy(b);
solveInPlace(rotateA, x, rotateb);
}
Givens() = delete;
~Givens() = delete;
};
#endif // GIVENS_HPP
#ifndef POLYNOMIAL_RECONSTRUCTION_HPP
#define POLYNOMIAL_RECONSTRUCTION_HPP
#include <mesh/MeshTraits.hpp>
#include <mesh/StencilBuilder.hpp>
#include <scheme/DiscreteFunctionDPk.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
#include <algebra/Givens.hpp>
#include <analysis/GaussLegendreQuadratureDescriptor.hpp>
#include <analysis/QuadratureFormula.hpp>
#include <analysis/QuadratureManager.hpp>
#include <geometry/LineTransformation.hpp>
#include <mesh/MeshData.hpp>
#include <mesh/MeshDataManager.hpp>
class PolynomialReconstruction
{
public:
template <MeshConcept MeshType, typename DataType>
PUGS_INLINE DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>>
build(const MeshType& mesh, const DiscreteFunctionP0<DataType> p0_function)
{
const size_t degree = 1;
const Stencil stencil = StencilBuilder{}.build(mesh.connectivity());
auto xr = mesh.xr();
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
auto cell_is_owned = mesh.connectivity().cellIsOwned();
DiscreteFunctionDPk<MeshType::Dimension, std::remove_const_t<DataType>> dpk{mesh.meshVariant(), degree};
if constexpr ((is_polygonal_mesh_v<MeshType>)and(MeshType::Dimension == 1)) {
auto qf = QuadratureManager::instance().getLineFormula(GaussLegendreQuadratureDescriptor(1));
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_j_id) {
if (cell_is_owned[cell_j_id]) {
auto stencil_cell_list = stencil[cell_j_id];
SmallVector<double> b(stencil_cell_list.size());
const double p_j = p0_function[cell_j_id];
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
const CellId cell_i_id = stencil_cell_list[i];
b[i] = p0_function[cell_i_id] - p_j;
}
SmallMatrix<double> A(stencil_cell_list.size(), degree);
using R1 = TinyVector<1>;
const R1& Xj = xj[cell_j_id];
auto X_Xj = [&](const R1& X) { return X - Xj; };
for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
const CellId cell_i_id = stencil_cell_list[i];
A(i, 0) = X_Xj(xj[cell_i_id])[0];
}
SmallVector<double> x(1);
Givens::solveInPlace(A, x, b);
dpk.coefficients(cell_j_id)[0] = p_j;
dpk.coefficients(cell_j_id)[1] = x[0];
}
});
}
synchronize(dpk.cellArrays());
return dpk;
}
PolynomialReconstruction() {}
};
#endif // POLYNOMIAL_RECONSTRUCTION_HPP
......@@ -96,6 +96,7 @@ add_executable (unit_tests
test_GaussLegendreQuadratureDescriptor.cpp
test_GaussLobattoQuadratureDescriptor.cpp
test_GaussQuadratureDescriptor.cpp
test_Givens.cpp
test_IfProcessor.cpp
test_IncDecExpressionProcessor.cpp
test_IntegrateCellArray.cpp
......@@ -203,6 +204,7 @@ add_executable (mpi_unit_tests
test_OFStream.cpp
test_ParallelChecker_read.cpp
test_Partitioner.cpp
test_PolynomialReconstruction.cpp
test_RandomEngine.cpp
test_StencilBuilder.cpp
test_SubItemArrayPerItemVariant.cpp
......
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
#include <algebra/Givens.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("Givens", "[algebra]")
{
SECTION("square matrix")
{
SmallMatrix<double> A{5, 5};
A.fill(0);
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;
// SmallMatrix<double> A{S.getMatrix()};
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;
SmallVector<double> x{5};
x = zero;
Givens::solve(A, x, b);
SmallVector<double> error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
}
SECTION("rectangular matrix")
{
SmallMatrix<double> A{5, 3};
A.fill(0);
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(3, 2) = -1;
A(4, 0) = 0.5;
A(4, 1) = 1;
A(4, 2) = 1;
// SmallMatrix<double> A{S.getMatrix()};
SmallVector<const double> x_exact = [] {
SmallVector<double> y{3};
y[0] = 1;
y[1] = 3;
y[2] = 2;
return y;
}();
SmallVector<double> b = A * x_exact;
SmallVector<double> x{3};
x = zero;
Givens::solve(A, x, b);
SmallVector<double> error = x - x_exact;
REQUIRE(std::sqrt(dot(error, error)) < 1E-10 * std::sqrt(dot(x, x)));
}
}
#include <catch2/catch_approx.hpp>
#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>
#include <Kokkos_Core.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/Types.hpp>
#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
#include <mesh/Mesh.hpp>
#include <mesh/MeshDataManager.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/PolynomialReconstruction.hpp>
#include <MeshDataBaseForTests.hpp>
// clazy:excludeall=non-pod-global-static
TEST_CASE("PolynomialReconstruction", "[scheme]")
{
SECTION("1D")
{
using R1 = TinyVector<1>;
for (auto named_mesh : MeshDataBaseForTests::get().all1DMeshes()) {
SECTION(named_mesh.name())
{
auto p_mesh = named_mesh.mesh()->get<Mesh<1>>();
auto& mesh = *p_mesh;
auto affine = [](const R1& x) { return 2.3 + 1.7 * x[0]; };
auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
DiscreteFunctionP0<double> fh{p_mesh};
parallel_for(
mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) { fh[cell_id] = affine(xj[cell_id]); });
PolynomialReconstruction reconstruction;
auto dpk = reconstruction.build(mesh, fh);
{
double max_mean_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
max_mean_error = std::max(max_mean_error, std::abs(dpk[cell_id](xj[cell_id]) - affine(xj[cell_id])));
}
REQUIRE(max_mean_error == Catch::Approx(0).margin(1E-14));
}
{
double max_slope_error = 0;
for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
const double reconstructed_slope =
(dpk[cell_id](R1{0.1} + xj[cell_id]) - dpk[cell_id](xj[cell_id] - R1{0.1})) / 0.2;
max_slope_error = std::max(max_slope_error, std::abs(reconstructed_slope - 1.7));
}
REQUIRE(max_slope_error == Catch::Approx(0).margin(1E-14));
}
}
}
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment