From c5a804150bf1e72164d10c452b5dfd113eb9a567 Mon Sep 17 00:00:00 2001
From: podlesny <podlesny@zedat.fu-berlin.de>
Date: Thu, 4 Feb 2021 12:33:03 +0100
Subject: [PATCH] robust implementation of eigenvalue computation, prevents
complex eigenvalues due to round-off errors
---
.../spatial-solving/tnnmg/functional.hh | 13 ++--
dune/tectonic/utils/CMakeLists.txt | 2 +
dune/tectonic/utils/maxeigenvalue.hh | 60 +++++++++++++++++++
3 files changed, 71 insertions(+), 4 deletions(-)
create mode 100644 dune/tectonic/utils/maxeigenvalue.hh
diff --git a/dune/tectonic/spatial-solving/tnnmg/functional.hh b/dune/tectonic/spatial-solving/tnnmg/functional.hh
index 3e9eab1a..2e76d0aa 100755
--- a/dune/tectonic/spatial-solving/tnnmg/functional.hh
+++ b/dune/tectonic/spatial-solving/tnnmg/functional.hh
@@ -16,6 +16,7 @@
#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+#include "../../utils/maxeigenvalue.hh"
#include "../../utils/debugutils.hh"
template<class M, class V, class N, class L, class U, class R>
@@ -423,18 +424,19 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
{
for (size_t i=0; i<this->quadraticPart().N(); ++i) {
const auto& Aii = this->quadraticPart()[i][i];
+ maxEigenvalues_[i] = maxEigenvalue(Aii);
// due to numerical imprecision, Aii might not be symmetric leading to complex eigenvalues
// consider Toeplitz decomposition of Aii and take only symmetric part
- auto symBlock = Aii;
+ /*auto symBlock = Aii;
- /*for (size_t j=0; j<symBlock.N(); j++) {
+ for (size_t j=0; j<symBlock.N(); j++) {
for (size_t k=0; k<symBlock.M(); k++) {
if (symBlock[j][k]/symBlock[j][j] < 1e-14)
symBlock[j][k] = 0;
}
- }*/
+ }
try {
typename Vector::block_type eigenvalues;
@@ -442,8 +444,11 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
maxEigenvalues_[i] =
*std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
} catch (Dune::Exception &e) {
+ print(symBlock, "symBlock");
+ typename Vector::block_type eigenvalues;
+ Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
std::cout << "complex eig in dof: " << i << std::endl;
- }
+ } */
}
}
diff --git a/dune/tectonic/utils/CMakeLists.txt b/dune/tectonic/utils/CMakeLists.txt
index c5498570..9d8acab4 100644
--- a/dune/tectonic/utils/CMakeLists.txt
+++ b/dune/tectonic/utils/CMakeLists.txt
@@ -4,6 +4,7 @@ add_custom_target(tectonic_dune_utils SOURCES
diameter.hh
geocoordinate.hh
index-in-sorted-range.hh
+ maxeigenvalue.hh
tobool.hh
)
@@ -14,5 +15,6 @@ install(FILES
diameter.hh
geocoordinate.hh
index-in-sorted-range.hh
+ maxeigenvalue.hh
tobool.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)
diff --git a/dune/tectonic/utils/maxeigenvalue.hh b/dune/tectonic/utils/maxeigenvalue.hh
new file mode 100644
index 00000000..47cbcb5a
--- /dev/null
+++ b/dune/tectonic/utils/maxeigenvalue.hh
@@ -0,0 +1,60 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+#define DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+
+/** \file
+ * \brief max eigenvalue computations for the FieldMatrix class
+ */
+
+#include <iostream>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 1, 1>& matrix) {
+ return matrix[0][0];
+}
+
+/** \brief calculates the eigenvalues of a symmetric field matrix
+ \param[in] matrix matrix eigenvalues are calculated for
+ \param[out] max eigenvalue
+*/
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 2, 2>& matrix) {
+ const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
+ const K p2 = p - matrix[1][1];
+ K q = p2 * p2 + matrix[1][0] * matrix[0][1];
+ if( q < 0 && q > -1e-14 ) q = 0;
+ if (q < 0) {
+ std::cout << matrix << std::endl;
+ // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
+ DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
+ }
+
+ // get square root
+ q = std::sqrt(q);
+
+ return p + q;
+}
+
+/** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
+ \param[in] matrix eigenvalues are calculated for
+ \param[out] max eigenvalue
+
+ \note If the input matrix is not symmetric the behavior of this method is undefined.
+ */
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 3, 3>& matrix) {
+ Dune::FieldVector<K, 3> eigenvalues;
+ Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+ return eigenvalues[0];
+}
+
+ /** @} end documentation */
+
+#endif
--
GitLab