From 611784cbf79629c45aa7c5fc02ce2e31b7cb2339 Mon Sep 17 00:00:00 2001
From: Max Kahnt <max.kahnt@fu-berlin.de>
Date: Sat, 23 Sep 2017 09:20:50 +0200
Subject: [PATCH] Add test for energyNorm.
---
dune/solvers/test/CMakeLists.txt | 1 +
dune/solvers/test/energynormtest.cc | 276 ++++++++++++++++++++++++++++
2 files changed, 277 insertions(+)
create mode 100644 dune/solvers/test/energynormtest.cc
diff --git a/dune/solvers/test/CMakeLists.txt b/dune/solvers/test/CMakeLists.txt
index 22af4820..1c7553ed 100644
--- a/dune/solvers/test/CMakeLists.txt
+++ b/dune/solvers/test/CMakeLists.txt
@@ -1,5 +1,6 @@
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)
diff --git a/dune/solvers/test/energynormtest.cc b/dune/solvers/test/energynormtest.cc
new file mode 100644
index 00000000..f56eaf35
--- /dev/null
+++ b/dune/solvers/test/energynormtest.cc
@@ -0,0 +1,276 @@
+#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;
+ }
+
+ std::cout << "...done.\n" << std::endl;
+}
+
+template<class... Args>
+struct NestableMultiTypeBlockVector : 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 NestableMultiTypeBlockMatrix : public Dune::MultiTypeBlockMatrix<Args...> {
+ constexpr static size_t blocklevel = 1; // fake value needed for BCRSMatrix nesting
+};
+
+// inject field traits for custom multi type block vector to be used by the norms
+namespace Dune {
+template<class T, class... TT>
+struct FieldTraits<NestableMultiTypeBlockVector<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<NestableMultiTypeBlockVector<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 = MultiTypeBlockVector<BVector, BVector>;
+ using MMatrix0 = MultiTypeBlockMatrix<BMatrix, BMatrix>;
+ using MMatrix1 = MultiTypeBlockMatrix<BMatrix, BMatrix>;
+ using MMatrix = MultiTypeBlockMatrix<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 = NestableMultiTypeBlockVector<FVector, FVector>;
+ using NMatrix0 = NestableMultiTypeBlockMatrix<FMatrix, FMatrix>;
+ using NMatrix1 = NestableMultiTypeBlockMatrix<FMatrix, FMatrix>;
+ using NMatrix = NestableMultiTypeBlockMatrix<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);
+ NMatrix nMatrix;
+ fillDiagonalStatic(nMatrix, fMatrix);
+ BNMatrix bnMatrix;
+ fillDiagonal(bnMatrix, nMatrix, bnSize);
+
+// TODO DEBUG
+using namespace Dune::Indices;
+ std::cout << "Entry 0" << bnMatrix[0][0][_0][_0][0][0] << std::endl;
+ std::cout << "Entry 0" << bnMatrix[0][0][_1][_1][0][0] << std::endl;
+ // 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;
+}
+
+
--
GitLab