3.4
DUNE for Multi-{Phase, Component, Scale, Physics, ...} flow and transport in porous media
preconditioners.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*****************************************************************************
4 * See the file COPYING for full copying permissions. *
5 * *
6 * This program is free software: you can redistribute it and/or modify *
7 * it under the terms of the GNU General Public License as published by *
8 * the Free Software Foundation, either version 3 of the License, or *
9 * (at your option) any later version. *
10 * *
11 * This program is distributed in the hope that it will be useful, *
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
14 * GNU General Public License for more details. *
15 * *
16 * You should have received a copy of the GNU General Public License *
17 * along with this program. If not, see <http://www.gnu.org/licenses/>. *
18 *****************************************************************************/
24#ifndef DUMUX_LINEAR_PRECONDITIONERS_HH
25#define DUMUX_LINEAR_PRECONDITIONERS_HH
26
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>
33
34#if HAVE_UMFPACK
35#include <dune/istl/umfpack.hh>
36#endif
37
41
42namespace Dumux {
43
78template<class M, class X, class Y, int l = 1>
79class SeqUzawa : public Dune::Preconditioner<X,Y>
80{
81 static_assert(Dumux::isMultiTypeBlockMatrix<M>::value && M::M() == 2 && M::N() == 2, "SeqUzawa expects a 2x2 MultiTypeBlockMatrix.");
82 static_assert(l== 1, "SeqUzawa expects a block level of 1.");
83
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])>;
86
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>;
91
92public:
94 using matrix_type = M;
96 using domain_type = X;
98 using range_type = Y;
100 using field_type = typename X::field_type;
102 using scalar_field_type = Dune::Simd::Scalar<field_type>;
103
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())
113#else
114 SeqUzawa(const M& mat, const Dune::ParameterTree& params)
115 : matrix_(mat)
116#endif
117 , numIterations_(params.get<std::size_t>("iterations"))
118 , relaxationFactor_(params.get<scalar_field_type>("relaxation"))
119 , verbosity_(params.get<int>("verbosity"))
120 , paramGroup_(params.get<std::string>("ParameterGroup"))
121 , useDirectVelocitySolverForA_(getParamFromGroup<bool>(paramGroup_, "LinearSolver.Preconditioner.DirectSolverForA", false))
122 {
123 const bool determineRelaxationFactor = getParamFromGroup<bool>(paramGroup_, "LinearSolver.Preconditioner.DetermineRelaxationFactor", true);
124
125 // AMG is needed for determination of omega
126 if (determineRelaxationFactor || !useDirectVelocitySolverForA_)
127 initAMG_(params);
128
129 if (useDirectVelocitySolverForA_)
130 initUMFPack_();
131
132 if (determineRelaxationFactor)
133 relaxationFactor_ = estimateOmega_();
134 }
135
139 virtual void pre(X& x, Y& b) {}
140
147 virtual void apply(X& update, const Y& currentDefect)
148 {
149 using namespace Dune::Indices;
150
151 auto& A = matrix_[_0][_0];
152 auto& B = matrix_[_0][_1];
153 auto& C = matrix_[_1][_0];
154 auto& D = matrix_[_1][_1];
155
156 const auto& f = currentDefect[_0];
157 const auto& g = currentDefect[_1];
158 auto& u = update[_0];
159 auto& p = update[_1];
160
161 // incorporate Dirichlet cell values
162 // TODO: pass Dirichlet constraint handler from outside
163 for (std::size_t i = 0; i < D.N(); ++i)
164 {
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()];
169 }
170
171 // the actual Uzawa iteration
172 for (std::size_t k = 0; k < numIterations_; ++k)
173 {
174 // u_k+1 = u_k + Q_A^−1*(f − (A*u_k + B*p_k)),
175 auto uRhs = f;
176 A.mmv(u, uRhs);
177 B.mmv(p, uRhs);
178 auto uIncrement = u;
179 applySolverForA_(uIncrement, uRhs);
180 u += uIncrement;
181
182 // p_k+1 = p_k + omega*(g - C*u_k+1 - D*p_k)
183 auto pIncrement = g;
184 C.mmv(u, pIncrement);
185 D.mmv(p, pIncrement);
186 pIncrement *= relaxationFactor_;
187 p += pIncrement;
188
189 if (verbosity_ > 1)
190 {
191 std::cout << "Uzawa iteration " << k
192 << ", residual: " << uRhs.two_norm() + pIncrement.two_norm()/relaxationFactor_ << std::endl;
193 }
194 }
195 }
196
200 virtual void post(X& x) {}
201
203 virtual Dune::SolverCategory::Category category() const
204 {
205 return Dune::SolverCategory::sequential;
206 }
207
208private:
209
210 void initAMG_(const Dune::ParameterTree& params)
211 {
212 using namespace Dune::Indices;
213 auto linearOperator = std::make_shared<LinearOperator>(matrix_[_0][_0]);
214 amgSolverForA_ = std::make_unique<AMGSolverForA>(linearOperator, params);
215 }
216
217 void initUMFPack_()
218 {
219#if HAVE_UMFPACK
220 using namespace Dune::Indices;
221 umfPackSolverForA_ = std::make_unique<Dune::UMFPack<A>>(matrix_[_0][_0]);
222#else
223 DUNE_THROW(Dune::InvalidStateException, "UMFPack not available. Use LinearSolver.Preconditioner.DirectVelocitySolver = false.");
224#endif
225 }
226
246 scalar_field_type estimateOmega_() const
247 {
248 using namespace Dune::Indices;
249 auto& A = matrix_[_0][_0];
250 auto& B = matrix_[_0][_1];
251 auto& C = matrix_[_1][_0];
252
253 U x(A.M());
254 x = 1.0;
255
256 scalar_field_type omega = 0.0;
257 scalar_field_type lambdaMax = 0.0;
258
259 static const auto iterations = Dumux::getParamFromGroup<std::size_t>(paramGroup_, "LinearSolver.Preconditioner.PowerLawIterations", 5);
260
261 // apply power iteration x_k+1 = M*x_k/|M*x_k| for the matrix M = -C*Ainv*B
262 for (std::size_t i = 0; i < iterations; ++i)
263 {
264 // bx = B*x
265 U bx(x.size());
266 B.mv(x, bx);
267
268 // ainvbx = Ainv*(B*x)
269 auto ainvbx = x;
270 applySolverForA_(ainvbx, bx);
271
272 // v = M*x = -C*(Ainv*B*x)
273 U v(x.size());
274 C.mv(ainvbx, v);
275 v *= -1.0;
276
277 // eigenvalue lambdaMax = xt*M*x/(xt*x) = xt*v/(xt*x);
278 lambdaMax = x.dot(v)/(x.dot(x));
279
280 // relaxation factor omega = 1/lambda;
281 omega = 1.0/lambdaMax;
282
283 // new iterate x = M*x/|M*x| = v/|v|
284 x = v;
285 x /= v.two_norm();
286 }
287
288 if (verbosity_ > 0)
289 {
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;
294 }
295
296 return omega;
297 }
298
299 template<class Sol, class Rhs>
300 void applySolverForA_(Sol& sol, Rhs& rhs) const
301 {
302 if (useDirectVelocitySolverForA_)
303 {
304#if HAVE_UMFPACK
305 Dune::InverseOperatorResult res;
306 umfPackSolverForA_->apply(sol, rhs, res);
307#endif
308 }
309 else
310 {
311 amgSolverForA_->pre(sol, rhs);
312 amgSolverForA_->apply(sol, rhs);
313 amgSolverForA_->post(sol);
314 }
315 }
316
318 const M& matrix_;
320 const std::size_t numIterations_;
322 scalar_field_type relaxationFactor_;
324 const int verbosity_;
325
326 std::unique_ptr<AMGSolverForA> amgSolverForA_;
327#if HAVE_UMFPACK
328 std::unique_ptr<Dune::UMFPack<A>> umfPackSolverForA_;
329#endif
330 const std::string paramGroup_;
331 const bool useDirectVelocitySolverForA_;
332};
333
334DUMUX_REGISTER_PRECONDITIONER("uzawa", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator<Dumux::SeqUzawa, 1>());
335
336} // end namespace Dumux
337
338#endif
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:358
Definition: adapt.hh:29
DUMUX_REGISTER_PRECONDITIONER("uzawa", Dumux::MultiTypeBlockMatrixPreconditionerTag, Dune::defaultPreconditionerBlockLevelCreator< Dumux::SeqUzawa, 1 >())
Definition: common/pdesolver.hh:36
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 &params)
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 &currentDefect)
Apply the preconditioner.
Definition: preconditioners.hh:147