Merge branch 'feature/energynorm-matrix-provider' into 'master'

Feature/energynorm matrix provider

See merge request !18

dune/solvers/norms/energynorm.hh

@@ -4,6 +4,7 @@
 #define DUNE_SOLVERS_NORMS_ENERGY_NORM_HH
 
 #include <cmath>
+#include <functional>
 
 #include <dune/matrix-vector/axy.hh>
 
@@ -39,29 +40,34 @@ namespace Solvers {
         /** \brief The type used for the result */
         using typename Base::field_type;
 
-        EnergyNorm(const field_type tol=1e-10 ) : iterationStep_(NULL), matrix_(NULL), tol_(tol) {}
+        EnergyNorm(const field_type tol=1e-10 ) : tol_(tol) {}
 
-        EnergyNorm(LinearIterationStep<MatrixType, VectorType>& it, const field_type tol=1e-10)
-            : iterationStep_(&it), matrix_(NULL), tol_(tol)
-        {}
+        template<class BV>
+        EnergyNorm(LinearIterationStep<MatrixType, VectorType, BV>& it, const field_type tol=1e-10)
+            : EnergyNorm(tol)
+        {
+            setIterationStep(&it);
+        }
 
         EnergyNorm(const MatrixType& matrix, const field_type tol=1e-10)
-            : iterationStep_(NULL), matrix_(&matrix), tol_(tol)
-        {}
+            : EnergyNorm(tol)
+        {
+            setMatrix(&matrix);
+        }
 
+        //! \brief sets the energy norm matrix
         void setMatrix(const MatrixType* matrix) {
-            matrix_ = matrix;
+            matrixProvider_ = [=]() -> const MatrixType& { return *matrix; };
         }
 
-        void setIterationStep(LinearIterationStep<MatrixType, VectorType>* iterationStep) {
-            iterationStep_ = iterationStep;
+        //! \brief sets to use the current problem matrix of the linear iteration step
+        template<class BV>
+        void setIterationStep(LinearIterationStep<MatrixType, VectorType, BV>* step) {
+            matrixProvider_ = [=]() -> const MatrixType& { return *step->getMatrix(); };
         }
 
         //! Compute the norm of the difference of two vectors
         field_type diff(const VectorType& f1, const VectorType& f2) const override {
-            if (iterationStep_ == NULL && matrix_ == NULL)
-                DUNE_THROW(Dune::Exception, "You have supplied neither a matrix nor an IterationStep to the EnergyNorm!");
-
             VectorType tmp_f = f1;
             tmp_f -= f2;
             return (*this)(tmp_f);
@@ -76,51 +82,37 @@ namespace Solvers {
         // \brief Compute the square of the norm of the given vector
         virtual field_type normSquared(const VectorType& f) const override
         {
-            if (iterationStep_ == NULL && matrix_ == NULL)
+            if (not matrixProvider_)
                 DUNE_THROW(Dune::Exception, "You have supplied neither a matrix nor an IterationStep to the EnergyNorm!");
-
-            const MatrixType& A = (iterationStep_)
-                ? *(iterationStep_->getMatrix())
-                : *matrix_;
-
-            const field_type ret = Dune::MatrixVector::Axy(A, f, f);
-
-            if (ret < 0)
-            {
-                if (ret < -tol_)
-                    DUNE_THROW(Dune::RangeError, "Supplied linear operator is not positive (semi-)definite: (u,Au) = " << ret);
-                return 0.0;
-            }
-
-            return ret;
+            const field_type ret = Dune::MatrixVector::Axy(matrixProvider_(), f, f);
+            return checkedValue(ret, tol_);
         }
 
         // \brief Compute the squared norm for a given vector and matrix
         DUNE_DEPRECATED static field_type normSquared(const VectorType& u,
                                                       const MatrixType& A,
-                                                      const double tol=1e-10)
+                                                      const field_type tol=1e-10)
         {
             const field_type ret = Dune::MatrixVector::Axy(A, u, u);
-
-            if (ret < 0)
-            {
-                if (ret < -tol)
-                    DUNE_THROW(Dune::RangeError, "Supplied linear operator is not positive (semi-)definite: (u,Au) = " << ret);
-                return 0.0;
-            }
-
-            return ret;
-
+            return checkedValue(ret, tol);
         }
 
     protected:
 
-        LinearIterationStep<MatrixType, VectorType>* iterationStep_;
-
-        const MatrixType* matrix_;
+        //! yields access to energy matrix
+        std::function<const MatrixType&()> matrixProvider_;
 
         const field_type tol_;
 
+        //! \brief throw an exception if value is below tolerance, project to R^+ otherwise.
+        static field_type checkedValue(field_type value, field_type tolerance)
+        {
+            if (value < 0)
+                if (value < -tolerance)
+                    DUNE_THROW(Dune::RangeError, "Supplied linear operator is not positive (semi-)definite: (u,Au) = " << value);
+                else return 0.0;
+            return value;
+        }
     };
 
 } /* namespace Solvers */

dune/solvers/test/CMakeLists.txt

 dune_add_test(SOURCES cgsteptest.cc)
 dune_add_test(SOURCES gssteptest.cc)
+dune_add_test(SOURCES energynormtest.cc)
 dune_add_test(SOURCES mmgtest.cc)
 dune_add_test(SOURCES multigridtest.cc)
 dune_add_test(SOURCES projectedgradienttest.cc)

dune/solvers/test/energynormtest.cc

+#include <config.h>
+
+#include <cmath>
+#include <iostream>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/typetraits.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+
+#include <dune/solvers/common/defaultbitvector.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+template <class M, class V, class BV>
+struct FakeLinearIterationStep : public LinearIterationStep<M, V, BV> {
+  using Base = LinearIterationStep<M, V, BV>;
+  using Base::Base;
+  void iterate() override {
+    DUNE_THROW(Dune::Exception, "This FakeLinearIterationStep must not be iterated.");
+  }
+};
+
+//! \brief tests with given matrix and vector
+template<class M, class V, class Field>
+bool testWith(const M& m, const V& v, Field normValue, const char* s) {
+  using namespace Dune;
+  using Norm = EnergyNorm<M, V>;
+
+  std::cout << "Test with " << s << "..." << std::endl;
+
+  // test unset matrix exception
+  std::cout << "Testing unset matrix exception.. ";
+  try {
+    Norm norm;
+    norm(v);
+    std::cout << "FAILURE." << std::endl;
+    std::cout << "...FAILED." << std::endl;
+    return false;
+  } catch (...) {
+    // TODO assert that the matrix exception is helpful
+    std::cout << " done." << std::endl;
+  }
+
+
+  // evaluation test lambda
+  auto evaluate = [&v, &normValue](auto&& norm) {
+    if(std::fabs(norm(v)-normValue) > 1e-10) {
+      std::cout << norm(v) <<  " vs " << normValue << std::endl;
+      DUNE_THROW(Exception, "Energy norm evaluation failed.");
+    }
+  };
+
+  try {
+  std::cout << "Testing matrix usage.. ";
+  // test on-construction matrix usage
+  {
+    Norm norm(m);
+    evaluate(norm);
+  }
+  // test setter matrix usage
+  {
+    Norm norm;
+    norm.setMatrix(&m);
+    evaluate(norm);
+  }
+  std::cout << "done." << std::endl;
+
+  std::cout << "Testing step usage.. ";
+  // test on-construction default step usage
+  using DefaultBV = Solvers::DefaultBitVector_t<V>;
+  auto defaultStep = FakeLinearIterationStep<M, V, DefaultBV>();
+  defaultStep.setMatrix(m);
+  {
+    Norm norm(defaultStep);
+    evaluate(norm);
+  }
+  // test setter default step usage
+  {
+    Norm norm;
+    norm.setIterationStep(&defaultStep);
+    evaluate(norm);
+  }
+  // test on-construction custom step usage
+  using CustomBV = double;
+  auto customStep = FakeLinearIterationStep<M, V, CustomBV>();
+  customStep.setMatrix(m);
+  {
+    Norm norm(customStep);
+    evaluate(norm);
+  }
+  // test setter custom step usage
+  {
+    Norm norm;
+    norm.setIterationStep(&customStep);
+    evaluate(norm);
+  }
+  std::cout << "done." << std::endl;
+
+  std::cout << "Testing reset usage.. ";
+  // test resetting matrix
+  {
+    M wrongM = m;
+    wrongM *= 2;
+    Norm norm(wrongM);
+    norm.setMatrix(&m);
+    evaluate(norm);
+  }
+  // test resetting step (with different (ignore node) type even)
+  {
+    auto wrongStep = customStep;
+    M wrongM = m;
+    wrongM *= 2;
+    wrongStep.setMatrix(wrongM);
+    Norm norm(wrongStep);
+    norm.setIterationStep(&defaultStep);
+    evaluate(norm);
+  }
+  std::cout << "done." << std::endl;
+
+  } catch(...) {
+    std::cout << "FAILURE." << std::endl;
+    std::cout << "...FAILED." << std::endl;
+    throw;
+    return false;
+  }
+
+  std::cout << "...done.\n" << std::endl;
+  return true;
+}
+
+template<class... Args>
+struct CustomMultiTypeBlockVector : public Dune::MultiTypeBlockVector<Args...> {
+  constexpr static size_t blocklevel = 1; // fake value need for BlockVector nesting
+  template<class K, typename = std::enable_if_t<Dune::IsNumber<K>::value>>
+  void operator=(K v) { Dune::Hybrid::forEach(*this, [v](auto& m) { m = v; }); } // allow init by scalar
+};
+
+template<class... Args>
+struct CustomMultiTypeBlockMatrix : public Dune::MultiTypeBlockMatrix<Args...> {
+  constexpr static size_t blocklevel = 1; // fake value needed for BCRSMatrix nesting
+  template<class K, typename = std::enable_if_t<Dune::IsNumber<K>::value>>
+  void operator*=(K v) { Dune::Hybrid::forEach(*this, [v](auto& b) { b *= v; }); } // allow multiplication by scalar
+};
+
+// inject field traits for custom multi type block vector to be used by the norms
+namespace Dune {
+template<class T, class... TT>
+struct FieldTraits<CustomMultiTypeBlockVector<T, TT...>> {
+  using field_type = typename FieldTraits<T>::field_type;
+  using real_type = typename FieldTraits<T>::real_type;
+};
+}
+
+// inject fake(!) default bit vector for custom multi type block vector
+namespace Dune { namespace Solvers { namespace Imp {
+template<class T, class... TT>
+struct DefaultBitVector<CustomMultiTypeBlockVector<T, TT...>> {
+  using type = T;
+};
+}}}
+
+
+//! \brief create and fill diagonal of sparse matrix with blocks
+template<class Block>
+void fillDiagonal(Dune::BCRSMatrix<Block>& m, const Block& b, size_t size) {
+  Dune::MatrixIndexSet indices(size, size);
+  for(size_t i=0; i<size; ++i)
+    indices.add(i, i);
+  indices.exportIdx(m);
+  for(size_t i=0; i<size; ++i)
+    m[i][i] = b;
+}
+
+//! \brief fill diagonal of static multi type matrix with block
+template<class M, class Block>
+void fillDiagonalStatic(M& m, const Block& b) {
+  using namespace Dune::Hybrid;
+  forEach(m, [&](auto& row) {
+    forEach(row, [&](auto& col) {
+      col = b;
+    });
+  });
+  using namespace Dune::Indices;
+  m[_0][_1] = 0.0;
+  m[_1][_0] = 0.0;
+}
+
+
+int main() try
+{
+  bool passed = true;
+
+  using namespace Dune;
+  using Field = double;
+
+  /// test with simple vectors
+  // types
+  constexpr size_t fSize = 2;
+  using FVector = FieldVector<Field, fSize>;
+  using FMatrix = FieldMatrix<Field, fSize, fSize>;
+  // instance setup
+  FMatrix fMatrix = {{1.0, 0.0}, {0.0, 1.0}};
+  FVector fVector = {3.0, 4.0};
+  // test
+  passed = passed && testWith(fMatrix, fVector, 5.0, "simple vectors");
+
+  /// test with blocked vectors
+  // types
+  constexpr size_t bSize = 4;
+  using BVector = BlockVector<FVector>;
+  using BMatrix = BCRSMatrix<FMatrix>;
+  // instance setup
+  BMatrix bMatrix;
+  fillDiagonal(bMatrix, fMatrix, bSize);
+  BVector bVector(bSize);
+  for(auto& b : bVector)
+    b = fVector;
+  // test
+  passed = passed && testWith(bMatrix, bVector, std::sqrt(bSize)*5.0, "blocked vectors");
+
+  /// test with multi blocked vectors
+  // types
+  using MVector = CustomMultiTypeBlockVector<BVector, BVector>;
+  using MMatrix0 = CustomMultiTypeBlockMatrix<BMatrix, BMatrix>;
+  using MMatrix1 = CustomMultiTypeBlockMatrix<BMatrix, BMatrix>;
+  using MMatrix = CustomMultiTypeBlockMatrix<MMatrix0, MMatrix1>;
+  // instance setup
+  using namespace Dune::Hybrid;
+  MVector mVector;
+  forEach(mVector, [&](auto& m) { m = bVector; });
+  MMatrix mMatrix;
+  fillDiagonalStatic(mMatrix, bMatrix);
+  // test
+  passed = passed && testWith(mMatrix, mVector, sqrt(2)*sqrt(bSize)*5.0, "multi blocked vectors");
+
+  /// test with blocked multitype vectors
+  // types
+  using NVector = CustomMultiTypeBlockVector<FVector, FVector>;
+  using NMatrix0 = CustomMultiTypeBlockMatrix<FMatrix, FMatrix>;
+  using NMatrix1 = CustomMultiTypeBlockMatrix<FMatrix, FMatrix>;
+  using NMatrix = CustomMultiTypeBlockMatrix<NMatrix0, NMatrix1>;
+  constexpr size_t bnSize = 4;
+  using BNVector = BlockVector<NVector>;
+  using BNMatrix = BCRSMatrix<NMatrix>;
+  // instance setup
+  NVector nVector;
+  forEach(nVector, [&](auto& n) { n = fVector; });
+  BNVector bnVector(bnSize);
+  for(auto& bn : bnVector)
+    bn = nVector;
+  NMatrix nMatrix;
+  fillDiagonalStatic(nMatrix, fMatrix);
+  BNMatrix bnMatrix;
+  fillDiagonal(bnMatrix, nMatrix, bnSize);
+  // test
+  passed = passed && testWith(bnMatrix, bnVector, sqrt(bnSize)*sqrt(2)*5.0, "blocked multitype vectors");
+
+  if(not passed)
+    return 1;
+
+  return 0;
+} catch(Dune::Exception e) {
+  std::cout << e << std::endl;
+  return 1;
+} catch(...) {
+  std::cout << "Unknown exception thrown." << std::endl;
+  return 1;
+}
+
+