Skip to content
Snippets Groups Projects
Commit 611784cb authored by Max Kahnt's avatar Max Kahnt
Browse files

Add test for energyNorm.

parent 425efd72
No related branches found
No related tags found
No related merge requests found
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)
......
#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;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment