From 32220435e1b20f3ffba1044d72b199eae63d1aa8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Mon, 8 Jul 2024 14:23:52 +0200
Subject: [PATCH] Add row weighting for least-square resolution
The idea is to add weight for the values that are provided by cells
that are close to the reconstruction
---
src/scheme/PolynomialReconstruction.cpp | 11 +++++++++
src/scheme/PolynomialReconstruction.hpp | 31 +++++++++++++++++++++++++
tests/test_PolynomialReconstruction.cpp | 2 +-
3 files changed, 43 insertions(+), 1 deletion(-)
diff --git a/src/scheme/PolynomialReconstruction.cpp b/src/scheme/PolynomialReconstruction.cpp
index 0e10c1bf8..788c6df64 100644
--- a/src/scheme/PolynomialReconstruction.cpp
+++ b/src/scheme/PolynomialReconstruction.cpp
@@ -278,6 +278,17 @@ PolynomialReconstruction::_build(
}
}
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
+ for (size_t l = 0; l < basis_dimension - 1; ++l) {
+ A(i, l) *= weight;
+ }
+ for (size_t l = 0; l < number_of_columns; ++l) {
+ B(i, l) *= weight;
+ }
+ }
+
const SmallMatrix<double>& X = X_pool[t];
Givens::solveCollectionInPlace(A, X, B);
diff --git a/src/scheme/PolynomialReconstruction.hpp b/src/scheme/PolynomialReconstruction.hpp
index 454e77c9d..78f7bc5ba 100644
--- a/src/scheme/PolynomialReconstruction.hpp
+++ b/src/scheme/PolynomialReconstruction.hpp
@@ -100,6 +100,15 @@ class PolynomialReconstruction
}
}
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
+ for (size_t l = 0; l < MeshType::Dimension; ++l) {
+ A(i, l) *= weight;
+ }
+ b[i] *= weight;
+ }
+
SmallVector<double> x(MeshType::Dimension);
Givens::solveInPlace(A, x, b);
@@ -150,6 +159,17 @@ class PolynomialReconstruction
}
}
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
+ for (size_t l = 0; l < MeshType::Dimension - 1; ++l) {
+ A(i, l) *= weight;
+ }
+ for (size_t l = 0; l < DataType::Dimension; ++l) {
+ B(i, l) *= weight;
+ }
+ }
+
TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
Givens::solveCollectionInPlace(A, X, B);
@@ -205,6 +225,17 @@ class PolynomialReconstruction
}
}
+ for (size_t i = 0; i < stencil_cell_list.size(); ++i) {
+ const CellId cell_i_id = stencil_cell_list[i];
+ const double weight = 1. / l2Norm(xj[cell_i_id] - Xj);
+ for (size_t l = 0; l < MeshType::Dimension - 1; ++l) {
+ A(i, l) *= weight;
+ }
+ for (size_t l = 0; l < DataType::Dimension; ++l) {
+ B(i, l) *= weight;
+ }
+ }
+
TinyMatrix<MeshType::Dimension, DataType::Dimension> X;
Givens::solveCollectionInPlace(A, X, B);
diff --git a/tests/test_PolynomialReconstruction.cpp b/tests/test_PolynomialReconstruction.cpp
index 7ce59a910..e019cc016 100644
--- a/tests/test_PolynomialReconstruction.cpp
+++ b/tests/test_PolynomialReconstruction.cpp
@@ -506,7 +506,7 @@ TEST_CASE("PolynomialReconstruction", "[scheme]")
max_y_slope_error = std::max(max_y_slope_error, l2Norm(reconstructed_slope - R3{-2.2, 1.3, -1.1}));
}
- REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-14));
+ REQUIRE(max_y_slope_error == Catch::Approx(0).margin(1E-13));
}
}
}
--
GitLab