Skip to content
Snippets Groups Projects
Commit 265acf0a authored by Elias Pipping's avatar Elias Pipping
Browse files

Tests: Add and use SymmetricSampleProblem class 

Such sample problems are set up in more than one place.
This change avoids code duplication
parent aec19234
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/solvers/test/cgsteptest.cc
+ 43
− 58
View file @ 265acf0a
......@@ -61,52 +61,24 @@ struct CGTestSuite
bool passed = true;
typedef typename Dune::FieldMatrix<double, 1, 1> MatrixBlock;
typedef typename Dune::FieldVector<double, 1> VectorBlock;
typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
typedef typename Dune::BlockVector<VectorBlock> Vector;
typedef typename Dune::BitSetVector<1> BitVector;
typedef typename GridType::LevelGridView GridView;
const GridView gridView = grid.levelGridView(grid.maxLevel());
Matrix A;
constructPQ1Pattern(gridView, A);
A=0.0;
assemblePQ1Stiffness(gridView, A);
EnergyNorm<Matrix, Vector> energyNorm(A);
Vector u(A.N()), u_sol(A.N());
u = 0;
for (size_t i=0; i<u.size(); ++i)
u_sol[i] = (1.0*rand())/RAND_MAX;
Vector rhs(A.N());
A.mv(u_sol, rhs);
BitVector ignore(A.N());
ignore.unsetAll();
markBoundaryDOFs(gridView, ignore);
// set Dirichlet values
for (size_t i=0; i<u.size(); ++i)
if (ignore[i][0])
u[i] = u_sol[i];
using Problem =
SymmetricSampleProblem<1, typename GridType::LevelGridView>;
Problem p(grid.levelGridView(grid.maxLevel()));
const auto verbosity = Solver::REDUCED;
const bool relativeErrors = false;
Analyser<Vector> analyser(u_sol, energyNorm, solveTol);
Analyser<typename Problem::Vector> analyser(p.u_ex, p.energyNorm,
solveTol);
{
std::string const testCase = "CGSolver ";
Vector u_copy = u;
Vector rhs_copy = rhs;
typename Problem::Vector u_copy = p.u;
typename Problem::Vector rhs_copy = p.rhs;
CGSolver<Matrix,Vector> solver(&A, &u_copy, &rhs_copy, nullptr,
maxIterations, stepTol, &energyNorm,
verbosity, relativeErrors);
CGSolver<typename Problem::Matrix, typename Problem::Vector> solver(
&p.A, &u_copy, &rhs_copy, nullptr, maxIterations, stepTol,
&p.energyNorm, verbosity, relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
......@@ -116,16 +88,19 @@ struct CGTestSuite
{
std::string const testCase = "CGSolver, preconditioned ";
Vector u_copy = u;
Vector rhs_copy = rhs;
typename Problem::Vector u_copy = p.u;
typename Problem::Vector rhs_copy = p.rhs;
using BGS = BlockGSStep<Matrix,Vector,BitVector>;
using BGS =
BlockGSStep<typename Problem::Matrix, typename Problem::Vector,
typename Problem::BitVector>;
BGS blockgs(BGS::Direction::SYMMETRIC);
blockgs.ignoreNodes_ = new BitVector(u.size(), false);
blockgs.ignoreNodes_ =
new typename Problem::BitVector(p.u.size(), false);
CGSolver<Matrix,Vector> solver(&A, &u_copy, &rhs_copy, &blockgs,
maxIterations, stepTol, &energyNorm,
verbosity, relativeErrors);
CGSolver<typename Problem::Matrix, typename Problem::Vector> solver(
&p.A, &u_copy, &rhs_copy, &blockgs, maxIterations, stepTol,
&p.energyNorm, verbosity, relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
......@@ -135,12 +110,15 @@ struct CGTestSuite
{
std::string const testCase = "Dune::Solvers::CGStep ";
Vector u_copy = u;
Vector rhs_copy = rhs;
typename Problem::Vector u_copy = p.u;
typename Problem::Vector rhs_copy = p.rhs;
Dune::Solvers::CGStep<Matrix,Vector> cgStep(A, u_copy, rhs_copy);
::LoopSolver<Vector> solver(&cgStep, maxIterations, stepTol, &energyNorm,
verbosity, relativeErrors);
Dune::Solvers::CGStep<typename Problem::Matrix,
typename Problem::Vector> cgStep(p.A, u_copy,
rhs_copy);
::LoopSolver<typename Problem::Vector> solver(
&cgStep, maxIterations, stepTol, &p.energyNorm, verbosity,
relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
......@@ -150,16 +128,23 @@ struct CGTestSuite
{
std::string const testCase = "Dune::Solvers::CGStep, preconditioned";
Vector u_copy = u;
Vector rhs_copy = rhs;
typename Problem::Vector u_copy = p.u;
typename Problem::Vector rhs_copy = p.rhs;
using BGS = BlockGSStep<Matrix,Vector,BitVector>;
using BGS =
BlockGSStep<typename Problem::Matrix, typename Problem::Vector,
typename Problem::BitVector>;
BGS blockgs(BGS::Direction::SYMMETRIC);
blockgs.ignoreNodes_ = new BitVector(u.size(), false);
Dune::Solvers::CGStep<Matrix,Vector> cgStep(A, u_copy, rhs_copy, blockgs);
::LoopSolver<Vector> solver(&cgStep, maxIterations, stepTol, &energyNorm,
verbosity, relativeErrors);
blockgs.ignoreNodes_ =
new typename Problem::BitVector(p.u.size(), false);
Dune::Solvers::CGStep<typename Problem::Matrix,
typename Problem::Vector> cgStep(p.A, u_copy,
rhs_copy,
blockgs);
::LoopSolver<typename Problem::Vector> solver(
&cgStep, maxIterations, stepTol, &p.energyNorm, verbosity,
relativeErrors);
solver.check();
solver.preprocess();
solver.solve();
......
dune/solvers/test/common.hh
+ 55
− 0
View file @ 265acf0a
......@@ -348,6 +348,61 @@ void markBoundaryDOFs(const GridView& gridView, BitVector& isBoundary)
}
}
template <size_t blocksize, class GridView>
class SymmetricSampleProblem {
public:
typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
typedef typename Dune::BlockVector<VectorBlock> Vector;
typedef typename Dune::BitSetVector<blocksize> BitVector;
typedef EnergyNorm<Matrix, Vector> Norm;
SymmetricSampleProblem(GridView const& gridView) {
constructPQ1Pattern(gridView, A);
A = 0.0;
assemblePQ1Stiffness(gridView, A);
energyNorm.setMatrix(&A);
ignore.resize(A.N());
ignore.unsetAll();
markBoundaryDOFs(gridView, ignore);
u.resize(A.N());
u_ex.resize(A.N());
rhs.resize(A.N());
randomize();
}
void randomize() {
// Set up a random 'solution'
u = 0;
for (size_t i = 0; i < u.size(); ++i)
for (size_t j = 0; j < blocksize; ++j) {
// in case of no dirichlet values you can make the
// matrix pd hereby
// A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
u_ex[i][j] = (1.0 * rand()) / RAND_MAX;
}
// Construct right hand side corresponding to the 'solution'
A.mv(u_ex, rhs);
// set dirichlet values
for (size_t i = 0; i < u.size(); ++i)
for (size_t j = 0; j < blocksize; ++j)
if (ignore[i][j])
u[i][j] = u_ex[i][j];
}
Matrix A;
Vector u;
Vector u_ex;
Vector rhs;
BitVector ignore;
Norm energyNorm;
};
template<class GridType, class TestSuite>
bool checkWithGrid(TestSuite& suite, int uniformRefinements)
......
dune/solvers/test/multigridtest.cc
+ 20
− 72
View file @ 265acf0a
......@@ -43,45 +43,13 @@ struct MultigridTestSuite
bool passed = true;
typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
typedef typename Dune::BlockVector<VectorBlock> Vector;
typedef typename Dune::BitSetVector<blocksize> BitVector;
using Problem =
SymmetricSampleProblem<1, typename GridType::LevelGridView>;
Problem p(grid.levelGridView(grid.maxLevel()));
typedef typename GridType::LevelGridView GridView;
const GridView gridView = grid.levelGridView(grid.maxLevel());
Matrix A;
constructPQ1Pattern(gridView, A);
A=0.0;
assemblePQ1Stiffness(gridView, A);
Vector u(A.N()), u_ex(A.N());
u = 0;
for (size_t i=0; i<u.size(); ++i)
for (size_t j=0; j<blocksize; ++j)
{
u_ex[i][j] = (1.0*rand())/RAND_MAX;
// A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]); // in case of no dirichlet values you can make the matrix pd hereby
}
Vector rhs(A.N());
A.mv(u_ex, rhs);
BitVector ignore(A.N());
ignore.unsetAll();
markBoundaryDOFs(gridView, ignore);
// set dirichlet values
for (size_t i=0; i<u.size(); ++i)
for (size_t j=0; j<blocksize; ++j)
if (ignore[i][j])
u[i][j] = u_ex[i][j];
typedef EnergyNorm<Matrix, Vector> Norm;
typedef typename Problem::Vector Vector;
typedef typename Problem::Matrix Matrix;
typedef typename Problem::BitVector BitVector;
typedef ::LoopSolver<Vector> Solver;
typedef MultigridStep<Matrix, Vector, BitVector> MGStep;
typedef MultigridTransfer<Vector, BitVector, Matrix> Transfer;
......@@ -98,71 +66,51 @@ struct MultigridTestSuite
transfer[i] = t;
}
// we want to control errors with respect to the energy norm induced by the matrix
Norm energynorm(A);
// set up smoothers and basesolver
Smoother smoother, basesolver_step;
Norm basenorm(basesolver_step);
typename Problem::Norm basenorm(basesolver_step);
Solver basesolver(&basesolver_step, 1, 1e-10, &basenorm, Solver::QUIET);
// create iteration step and set all kinds of data
MGStep mgStep;
mgStep.setTransferOperators(transfer);
mgStep.setProblem(A,u,rhs);
mgStep.setProblem(p.A,p.u,p.rhs);
mgStep.setSmoother(&smoother);
mgStep.setMGType(1,3,3);
mgStep.ignoreNodes_ = &ignore;
mgStep.ignoreNodes_ = &p.ignore;
mgStep.basesolver_ = &basesolver;
// create loop solver
Solver solver(&mgStep, maxIterations, tol, &energynorm, Solver::FULL);
Solver solver(&mgStep, maxIterations, tol, &p.energyNorm, Solver::FULL);
// solve problem
solver.preprocess();
solver.solve();
u = mgStep.getSol();
p.u = mgStep.getSol();
if (energynorm.diff(u,u_ex)>tol*10)
if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
{
std::cout << "The MultigridStep did not produce a satisfactory result. ||u-u_ex||=" << energynorm.diff(u,u_ex) << std::endl;
std::cout << "||u_ex||=" << energynorm(u_ex) << std::endl;
std::cout << "The MultigridStep did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
passed = false;
}
// now solve with the same operator but different rhs (and dirichlet values)
// create new exact solution
for (size_t i=0; i<u.size(); ++i)
for (size_t j=0; j<blocksize; ++j)
u_ex[i][j] = (1.0*rand())/RAND_MAX;
// ... and its corresponding rhs
A.mv(u_ex, rhs);
// set initial iterate (remember that the iterationStep member points to u)
u = 0.0;
// set dirichlet values
for (size_t i=0; i<u.size(); ++i)
for (size_t j=0; j<blocksize; ++j)
if (ignore[i][j])
u[i][j] = u_ex[i][j];
p.randomize();
// reset rhs in mgStep without calling preprocess subsequently
mgStep.setRhs(rhs);
mgStep.setRhs(p.rhs);
solver.solve();
u = mgStep.getSol();
p.u = mgStep.getSol();
if (energynorm.diff(u,u_ex)>tol*10)
if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
{
std::cout << "The MultigridStep did not produce a satisfactory result after calling setRhs(). ||u-u_ex||=" << energynorm.diff(u,u_ex) << std::endl;
std::cout << "||u_ex||=" << energynorm(u_ex) << std::endl;
std::cout << "The MultigridStep did not produce a satisfactory result after calling setRhs(). ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
passed = false;
}
......
dune/solvers/test/umfpacksolvertest.cc
+ 11
− 49
View file @ 265acf0a
......@@ -34,52 +34,14 @@ struct UMFPackSolverTestSuite
bool passed = true;
typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
typedef typename Dune::BlockVector<VectorBlock> Vector;
typedef typename Dune::BitSetVector<blocksize> BitVector;
typedef EnergyNorm<Matrix, Vector> Norm;
const auto gridView = grid.leafGridView();
Matrix A;
constructPQ1Pattern(gridView, A);
A=0.0;
assemblePQ1Stiffness(gridView, A);
// Set up a random 'solution'
Vector u(A.N()), u_ex(A.N());
u = 0;
for (size_t i=0; i<u.size(); ++i)
for (size_t j=0; j<blocksize; ++j)
{
u_ex[i][j] = (1.0*rand())/RAND_MAX;
// A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]); // in case of no dirichlet values you can make the matrix pd hereby
}
// Construct right hand side corresponding to the 'solution'
Vector rhs(A.N());
A.mv(u_ex, rhs);
BitVector ignore(A.N());
ignore.unsetAll();
markBoundaryDOFs(gridView, ignore);
// set dirichlet values
for (size_t i=0; i<u.size(); ++i)
for (size_t j=0; j<blocksize; ++j)
if (ignore[i][j])
u[i][j] = u_ex[i][j];
// we want to control errors with respect to the energy norm induced by the matrix
Norm energynorm(A);
typedef Solvers::UMFPackSolver<Matrix,Vector> Solver;
// create loop solver
Solver solver(A,u,rhs);
solver.ignoreNodes_ = &ignore;
using Problem =
SymmetricSampleProblem<blocksize, typename GridType::LeafGridView>;
Problem p(grid.leafGridView());
typedef Solvers::UMFPackSolver<typename Problem::Matrix,
typename Problem::Vector> Solver;
Solver solver(p.A,p.u,p.rhs);
solver.ignoreNodes_ = &p.ignore;
// solve problem
solver.preprocess();
......@@ -87,10 +49,10 @@ struct UMFPackSolverTestSuite
//std::cout << u_ex << std::endl;
if (energynorm.diff(u,u_ex)>tol*10)
if (p.energyNorm.diff(p.u,p.u_ex)>tol*10)
{
std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << energynorm.diff(u,u_ex) << std::endl;
std::cout << "||u_ex||=" << energynorm(u_ex) << std::endl;
std::cout << "The UMFPackSolver did not produce a satisfactory result. ||u-u_ex||=" << p.energyNorm.diff(p.u,p.u_ex) << std::endl;
std::cout << "||u_ex||=" << p.energyNorm(p.u_ex) << std::endl;
passed = false;
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment