Commit 33c2455c authored by Ansgar Burchardt's avatar Ansgar Burchardt
Browse files

move implementation details of linear correction out of TNNMGStep

This simplifies the `iterate()` method such that it fits on a single
screen.
parent 5a811eb5
Pipeline #9798 passed with stage
in 4 minutes and 15 seconds
......@@ -2,6 +2,7 @@ install(FILES
acceleratednonlineargsstep.hh
genericnonlineargs.hh
genericnonlinearjacobi.hh
linearcorrection.hh
nonlineargsstep.hh
preconfiguredtnnmgstep.hh
tnnmgacceleration.hh
......
#ifndef DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH
#define DUNE_TNNMG_ITERATIONSTEPS_LINEARCORRECTION_HH 1
#include <functional>
#include <memory>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/solvers/linearsolver.hh>
namespace Dune {
namespace TNNMG {
namespace Impl {
template<typename Vector>
Solvers::DefaultBitVector_t<Vector>
emptyIgnore(const Vector& v)
{
// TNNMGStep assumes that the linearization and the solver for the
// linearized problem will not use the ignoreNodes field
Solvers::DefaultBitVector_t<Vector> ignore;
Solvers::resizeInitialize(ignore, v, false);
return ignore;
}
} /* namespace Impl */
/**
* \brief linear correction step for use by \ref TNNMGStep
*
* The function object should solve the linear equation \f$ A x = b \f$.
* Be aware that full rows or columns of `A` might contain only zeroes.
*/
template<typename Matrix, typename Vector>
using LinearCorrection = std::function<void(const Matrix& A, Vector& x, const Vector& b)>;
template<typename Matrix, typename Vector>
LinearCorrection<Matrix, Vector>
makeLinearCorrection(std::shared_ptr< Solvers::LinearSolver<Matrix, Vector> > linearSolver)
{
return [=](const Matrix& A, Vector& x, const Vector& b) {
linearSolver->setProblem(A, x, b);
linearSolver->preprocess();
linearSolver->solve();
};
}
template<typename Matrix, typename Vector>
LinearCorrection<Matrix, Vector>
makeLinearCorrection(std::shared_ptr< Solvers::IterativeSolver<Vector> > iterativeSolver)
{
return [=](const Matrix& A, Vector& x, const Vector& b) {
using LinearIterationStep = Solvers::LinearIterationStep<Matrix, Vector>;
auto emptyIgnore = Impl::emptyIgnore(x);
auto linearIterationStep = dynamic_cast<LinearIterationStep*>(&iterativeSolver->getIterationStep());
if (not linearIterationStep)
DUNE_THROW(Exception, "iterative solver must use a linear iteration step");
linearIterationStep->setProblem(A, x, b);
iterativeSolver->preprocess();
iterativeSolver->solve();
};
}
template<typename Matrix, typename Vector>
LinearCorrection<Matrix, Vector>
makeLinearCorrection(std::shared_ptr< Solvers::LinearIterationStep<Matrix, Vector> > linearIterationStep, int nIterationSteps = 1)
{
return [=](const Matrix& A, Vector& x, const Vector& b) {
auto emptyIgnore = Impl::emptyIgnore(x);
linearIterationStep->setIgnore(emptyIgnore);
linearIterationStep->setProblem(A, x, b);
linearIterationStep->preprocess();
for (int i = 0; i < nIterationSteps; ++i)
linearIterationStep->iterate();
};
}
} /* namespace TNNMG */
} /* namespace Dune */
#endif
......@@ -14,6 +14,8 @@
#include <dune/solvers/solvers/iterativesolver.hh>
#include <dune/solvers/solvers/linearsolver.hh>
#include <dune/tnnmg/iterationsteps/linearcorrection.hh>
namespace Dune {
namespace TNNMG {
......@@ -58,7 +60,7 @@ public:
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
preSmoothingSteps_(1),
iterativeSolver_(iterativeSolver),
linearCorrection_(makeLinearCorrection(iterativeSolver)),
projection_(projection),
lineSolver_(lineSolver)
{}
......@@ -79,7 +81,7 @@ public:
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
preSmoothingSteps_(1),
linearSolver_(linearSolver),
linearCorrection_(makeLinearCorrection(linearSolver)),
projection_(projection),
lineSolver_(lineSolver)
{}
......@@ -101,8 +103,7 @@ public:
f_(&f),
nonlinearSmoother_(nonlinearSmoother),
preSmoothingSteps_(1),
linearIterationStep_(linearIterationStep),
noOfLinearIterationSteps_(noOfLinearIterationSteps),
linearCorrection_(makeLinearCorrection(linearIterationStep, noOfLinearIterationSteps)),
projection_(projection),
lineSolver_(lineSolver)
{}
......@@ -151,57 +152,7 @@ public:
Solvers::resizeInitializeZero(correction_, x);
Solvers::resizeInitializeZero(constrainedCorrection_, r);
// TNNMGStep assumes that the linearization and the solver for the
// linearized problem will not use the ignoreNodes field
auto emptyIgnore = ConstrainedBitVector();
Solvers::resizeInitialize(emptyIgnore, constrainedCorrection_, false);
// Do the linear correction with a LinearIterationStep object
if (linearIterationStep_)
{
linearIterationStep_->setIgnore(emptyIgnore);
linearIterationStep_->setProblem(A, constrainedCorrection_, r);
linearIterationStep_->preprocess();
for (unsigned int i=0; i<noOfLinearIterationSteps_; i++)
linearIterationStep_->iterate();
}
else if (iterativeSolver_) // Do the linear correction with an iterative Solver object
{
// Hand the linear problem to the iterative solver.
// The IterationStep member needs to be a LinearIterationStep,
// so we can give it the matrix.
using LinearIterationStepType = Solvers::LinearIterationStep<std::decay_t<decltype(A)>,
typename Linearization::ConstrainedVector,
decltype(emptyIgnore) >;
LinearIterationStepType* linearIterationStep;
auto iterativeSolver = std::dynamic_pointer_cast<Solvers::IterativeSolver<typename Linearization::ConstrainedVector> >(iterativeSolver_);
if (iterativeSolver)
{
iterativeSolver->getIterationStep().setIgnore(emptyIgnore);
linearIterationStep = dynamic_cast<LinearIterationStepType*>(&(iterativeSolver->getIterationStep()));
} else
DUNE_THROW(Exception, "Linear solver has to be an IterativeSolver!");
if (linearIterationStep)
linearIterationStep->setProblem(A, constrainedCorrection_, r);
else
DUNE_THROW(Exception, "Linear solver does not accept matrices!");
iterativeSolver_->preprocess();
iterativeSolver_->solve();
}
else // Do the linear correction with a linear Solver object
{
linearSolver_->setProblem(A,constrainedCorrection_, r);
linearSolver_->preprocess();
linearSolver_->solve();
}
linearCorrection_(A, constrainedCorrection_, r);
linearization_->extendCorrection(constrainedCorrection_, correction_);
// Project onto admissible set
......@@ -243,17 +194,8 @@ private:
std::shared_ptr<Linearization> linearization_;
// The following members cannot all be used at once:
// either we have a Dune::Solvers::IterativeSolver that implements the linear correction...
std::shared_ptr<IterativeSolver> iterativeSolver_;
// or we have a Dune::Solvers::LinearSolver that implements the linear correction...
std::shared_ptr<LinearSolver> linearSolver_;
// ... or we have a mere LinearIterationStep, together with a number of times
// it is supposed to be called. You cannot have more than one at once.
std::shared_ptr<LinearIterationStep<typename Linearization::ConstrainedMatrix,typename Linearization::ConstrainedVector> > linearIterationStep_;
unsigned int noOfLinearIterationSteps_;
//! \brief linear correction
LinearCorrection<ConstrainedMatrix, ConstrainedVector> linearCorrection_;
typename Linearization::ConstrainedVector constrainedCorrection_;
Vector correction_;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment