24#ifndef DUMUX_LINEAR_PRECONDITIONERS_HH
25#define DUMUX_LINEAR_PRECONDITIONERS_HH
27#include <dune/common/exceptions.hh>
28#include <dune/common/float_cmp.hh>
29#include <dune/common/indices.hh>
30#include <dune/common/version.hh>
31#include <dune/istl/preconditioners.hh>
32#include <dune/istl/paamg/amg.hh>
35#include <dune/istl/umfpack.hh>
78template<
class M,
class X,
class Y,
int l = 1>
79class SeqUzawa :
public Dune::Preconditioner<X,Y>
82 static_assert(l== 1,
"SeqUzawa expects a block level of 1.");
84 using A = std::decay_t<decltype(std::declval<M>()[Dune::Indices::_0][Dune::Indices::_0])>;
85 using U = std::decay_t<decltype(std::declval<X>()[Dune::Indices::_0])>;
87 using Comm = Dune::Amg::SequentialInformation;
88 using LinearOperator = Dune::MatrixAdapter<A, U, U>;
89 using Smoother = Dune::SeqSSOR<A, U, U>;
90 using AMGSolverForA = Dune::Amg::AMG<LinearOperator, U, Smoother, Comm>;
110#if DUNE_VERSION_GTE(DUNE_ISTL,2,8)
111 SeqUzawa(
const std::shared_ptr<
const Dune::AssembledLinearOperator<M,X,Y>>& op,
const Dune::ParameterTree& params)
112 : matrix_(op->getmat())
117 , numIterations_(params.get<std::size_t>(
"iterations"))
119 , verbosity_(params.get<int>(
"verbosity"))
120 , paramGroup_(params.get<std::string>(
"ParameterGroup"))
121 , useDirectVelocitySolverForA_(
getParamFromGroup<bool>(paramGroup_,
"LinearSolver.Preconditioner.DirectSolverForA", false))
123 const bool determineRelaxationFactor = getParamFromGroup<bool>(paramGroup_,
"LinearSolver.Preconditioner.DetermineRelaxationFactor",
true);
126 if (determineRelaxationFactor || !useDirectVelocitySolverForA_)
129 if (useDirectVelocitySolverForA_)
132 if (determineRelaxationFactor)
133 relaxationFactor_ = estimateOmega_();
139 virtual void pre(X& x, Y& b) {}
147 virtual void apply(X& update,
const Y& currentDefect)
149 using namespace Dune::Indices;
151 auto& A = matrix_[_0][_0];
152 auto& B = matrix_[_0][_1];
153 auto& C = matrix_[_1][_0];
154 auto& D = matrix_[_1][_1];
156 const auto& f = currentDefect[_0];
157 const auto& g = currentDefect[_1];
158 auto& u = update[_0];
159 auto& p = update[_1];
163 for (std::size_t i = 0; i < D.N(); ++i)
165 const auto& block = D[i][i];
166 for (
auto rowIt = block.begin(); rowIt != block.end(); ++rowIt)
167 if (Dune::FloatCmp::eq<scalar_field_type>(rowIt->one_norm(), 1.0))
168 p[i][rowIt.index()] = g[i][rowIt.index()];
172 for (std::size_t k = 0; k < numIterations_; ++k)
179 applySolverForA_(uIncrement, uRhs);
184 C.mmv(u, pIncrement);
185 D.mmv(p, pIncrement);
186 pIncrement *= relaxationFactor_;
191 std::cout <<
"Uzawa iteration " << k
192 <<
", residual: " << uRhs.two_norm() + pIncrement.two_norm()/relaxationFactor_ << std::endl;
203 virtual Dune::SolverCategory::Category
category()
const
205 return Dune::SolverCategory::sequential;
210 void initAMG_(
const Dune::ParameterTree& params)
212 using namespace Dune::Indices;
213 auto linearOperator = std::make_shared<LinearOperator>(matrix_[_0][_0]);
214 amgSolverForA_ = std::make_unique<AMGSolverForA>(linearOperator, params);
220 using namespace Dune::Indices;
221 umfPackSolverForA_ = std::make_unique<Dune::UMFPack<A>>(matrix_[_0][_0]);
223 DUNE_THROW(Dune::InvalidStateException,
"UMFPack not available. Use LinearSolver.Preconditioner.DirectVelocitySolver = false.");
248 using namespace Dune::Indices;
249 auto& A = matrix_[_0][_0];
250 auto& B = matrix_[_0][_1];
251 auto& C = matrix_[_1][_0];
259 static const auto iterations = Dumux::getParamFromGroup<std::size_t>(paramGroup_,
"LinearSolver.Preconditioner.PowerLawIterations", 5);
262 for (std::size_t i = 0; i < iterations; ++i)
270 applySolverForA_(ainvbx, bx);
278 lambdaMax = x.dot(v)/(x.dot(x));
281 omega = 1.0/lambdaMax;
290 std::cout <<
"\n*** Uzawa Preconditioner ***" << std::endl;
291 std::cout <<
"Estimating relaxation factor based on Schur complement" << std::endl;
292 std::cout <<
"Largest estimated eigenvalue lambdaMax = " << lambdaMax << std::endl;
293 std::cout <<
"Relaxation factor omega = 1/lambdaMax = " << omega << std::endl;
299 template<
class Sol,
class Rhs>
300 void applySolverForA_(Sol& sol, Rhs& rhs)
const
302 if (useDirectVelocitySolverForA_)
305 Dune::InverseOperatorResult res;
306 umfPackSolverForA_->apply(sol, rhs, res);
311 amgSolverForA_->pre(sol, rhs);
312 amgSolverForA_->apply(sol, rhs);
313 amgSolverForA_->post(sol);
320 const std::size_t numIterations_;
324 const int verbosity_;
326 std::unique_ptr<AMGSolverForA> amgSolverForA_;
328 std::unique_ptr<Dune::UMFPack<A>> umfPackSolverForA_;
330 const std::string paramGroup_;
331 const bool useDirectVelocitySolverForA_;
Type traits to be used with matrix types.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
The specialized Dumux macro and tag for the ISTL registry to choose the solver and preconditioner at ...
T getParamFromGroup(Args &&... args)
A free function to get a parameter from the parameter tree singleton with a model group.
Definition: parameters.hh:374
DUMUX_REGISTER_PRECONDITIONER("uzawa", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator< Dumux::SeqUzawa, 1 >())
Definition: common/pdesolver.hh:35
Helper type to determine whether a given type is a Dune::MultiTypeBlockMatrix.
Definition: matrix.hh:49
A preconditioner based on the Uzawa algorithm for saddle-point problems of the form .
Definition: preconditioners.hh:80
virtual Dune::SolverCategory::Category category() const
Category of the preconditioner (see SolverCategory::Category)
Definition: preconditioners.hh:203
M matrix_type
The matrix type the preconditioner is for.
Definition: preconditioners.hh:94
SeqUzawa(const M &mat, const Dune::ParameterTree ¶ms)
Constructor.
Definition: preconditioners.hh:114
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: preconditioners.hh:139
virtual void post(X &x)
Clean up.
Definition: preconditioners.hh:200
Dune::Simd::Scalar< field_type > scalar_field_type
Scalar type underlying the field_type.
Definition: preconditioners.hh:102
typename X::field_type field_type
The field type of the preconditioner.
Definition: preconditioners.hh:100
X domain_type
The domain type of the preconditioner.
Definition: preconditioners.hh:96
Y range_type
The range type of the preconditioner.
Definition: preconditioners.hh:98
virtual void apply(X &update, const Y ¤tDefect)
Apply the preconditioner.
Definition: preconditioners.hh:147