Loading...
Searching...
No Matches
uzawa_solver.hpp
Go to the documentation of this file.
#ifndef UZAWA_SOLVER_HPP_
#define UZAWA_SOLVER_HPP_
//==============================================================================
//==============================================================================
#include <dune/istl/solvers.hh>
#include <opm/elasticity/matrixops.hpp>
#include <opm/elasticity/mortar_utils.hpp>
namespace Opm {
namespace Elasticity {
template<class X, class Y>
class UzawaSolver : public Dune::InverseOperator<X,Y>
{
public:
typedef std::shared_ptr<Dune::InverseOperator<X,Y> > OperatorPtr;
UzawaSolver(OperatorPtr& innersolver_,
OperatorPtr& outersolver_,
const Matrix& B_) :
innersolver(innersolver_), outersolver(outersolver_), B(B_)
{
}
void apply(X& x, Y& b, double /* reduction */,
Dune::InverseOperatorResult& res) override
{
apply(x, b, res);
}
void apply(X& x, Y& b, Dune::InverseOperatorResult& res) override
{
Vector lambda, lambda2, u, u2;
MortarUtils::extractBlock(u, b, B.N());
Dune::InverseOperatorResult res2;
u2.resize(u.size());
u2 = 0;
innersolver->apply(u2, u, res2);
lambda.resize(B.M());
B.mtv(u2, lambda);
lambda2.resize(lambda.size());
lambda2 = 0;
outersolver->apply(lambda2, lambda, res2);
MortarUtils::extractBlock(u, b, B.N());
B.usmv(-1.0, lambda2, u);
u2 = 0;
innersolver->apply(u2, u, res);
MortarUtils::injectBlock(x, u2, B.N());
MortarUtils::injectBlock(x, lambda2, B.M(), B.N());
}
Dune::SolverCategory::Category category() const override
{
return Dune::SolverCategory::sequential;
}
protected:
OperatorPtr innersolver;
OperatorPtr outersolver;
const Matrix& B;
};
}
}
#endif
Opm::Elasticity::MortarUtils::extractBlock
static void extractBlock(Vector &x, const Vector &y, int len, int start=0)
Extract a range of indices from a vector.
Definition mortar_utils.hpp:27
Opm::Elasticity::MortarUtils::injectBlock
static void injectBlock(Vector &x, const Vector &y, int len, int start=0)
Inject a range of indices into a vector.
Definition mortar_utils.hpp:38
Opm::Elasticity::UzawaSolver::outersolver
OperatorPtr outersolver
The outer solver.
Definition uzawa_solver.hpp:85
Opm::Elasticity::UzawaSolver::B
const Matrix & B
The coupling matrix.
Definition uzawa_solver.hpp:86
Opm::Elasticity::UzawaSolver::apply
void apply(X &x, Y &b, double, Dune::InverseOperatorResult &res) override
Apply the scheme to a vector.
Definition uzawa_solver.hpp:47
Opm::Elasticity::UzawaSolver::apply
void apply(X &x, Y &b, Dune::InverseOperatorResult &res) override
Apply the scheme to a vector.
Definition uzawa_solver.hpp:57
Opm::Elasticity::UzawaSolver::UzawaSolver
UzawaSolver(OperatorPtr &innersolver_, OperatorPtr &outersolver_, const Matrix &B_)
Default constructor.
Definition uzawa_solver.hpp:35
Opm::Elasticity::UzawaSolver::innersolver
OperatorPtr innersolver
The inner solver.
Definition uzawa_solver.hpp:84
matrixops.hpp
Helper class with some matrix operations.
Opm::Elasticity::Matrix
Dune::BCRSMatrix< Dune::FieldMatrix< double, 1, 1 > > Matrix
A sparse matrix holding our operator.
Definition matrixops.hpp:27
Opm::Elasticity::Vector
Dune::BlockVector< Dune::FieldVector< double, 1 > > Vector
A vector holding our RHS.
Definition matrixops.hpp:33
mortar_utils.hpp
Mortar helper class.
Opm
Inverting small matrices.
Definition ImplicitAssembly.hpp:43