GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/linear/istlsolvers.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 143 155 92.3%
Functions: 1724 4205 41.0%
Branches: 265 752 35.2%
Line Branch Exec Source
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 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Linear
10 * \brief Linear solvers from dune-istl
11 */
12 #ifndef DUMUX_LINEAR_ISTL_SOLVERS_HH
13 #define DUMUX_LINEAR_ISTL_SOLVERS_HH
14
15 #include <memory>
16 #include <variant>
17
18 #include <dune/common/exceptions.hh>
19 #include <dune/common/shared_ptr.hh>
20 #include <dune/common/parallel/indexset.hh>
21 #include <dune/common/parallel/mpicommunication.hh>
22 #include <dune/grid/common/capabilities.hh>
23 #include <dune/istl/solvers.hh>
24 #include <dune/istl/solverfactory.hh>
25 #include <dune/istl/owneroverlapcopy.hh>
26 #include <dune/istl/scalarproducts.hh>
27 #include <dune/istl/paamg/amg.hh>
28 #include <dune/istl/paamg/pinfo.hh>
29
30 #include <dumux/common/typetraits/matrix.hh>
31 #include <dumux/common/typetraits/vector.hh>
32 #include <dumux/linear/linearalgebratraits.hh>
33 #include <dumux/linear/preconditioners.hh>
34 #include <dumux/linear/linearsolverparameters.hh>
35 #include <dumux/linear/matrixconverter.hh>
36 #include <dumux/linear/parallelhelpers.hh>
37 #include <dumux/linear/solvercategory.hh>
38 #include <dumux/linear/solver.hh>
39
40 #include <dune/istl/foreach.hh>
41
42 namespace Dumux::Detail::IstlSolvers {
43
44 /*!
45 * \ingroup Linear
46 * \brief Returns the block level for the preconditioner for a given matrix
47 *
48 * \tparam M The matrix.
49 */
50 template<class M>
51 constexpr std::size_t preconditionerBlockLevel() noexcept
52 {
53 return isMultiTypeBlockMatrix<M>::value ? 2 : 1;
54 }
55
56 template<template<class,class,class,int> class Preconditioner, int blockLevel = 1>
57 class IstlDefaultBlockLevelPreconditionerFactory
58 {
59 public:
60 template<class TL, class M>
61 ✗ auto operator() (TL typeList, const M& matrix, const Dune::ParameterTree& config)
62 {
63 using Matrix = typename Dune::TypeListElement<0, decltype(typeList)>::type;
64 using Domain = typename Dune::TypeListElement<1, decltype(typeList)>::type;
65 using Range = typename Dune::TypeListElement<2, decltype(typeList)>::type;
66
10/18
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 854 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 854 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 2464 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 2464 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 16029 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 16029 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 36 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 36 times.
✗ Branch 23 not taken.
 38770 std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner
67 = std::make_shared<Preconditioner<Matrix, Domain, Range, blockLevel>>(matrix, config);
68 ✗ return preconditioner;
69 }
70 };
71
72 template<template<class,class,class> class Preconditioner>
73 class IstlDefaultPreconditionerFactory
74 {
75 template<class TL, class M>
76 auto operator() (TL typeList, const M& matrix, const Dune::ParameterTree& config)
77 {
78 using Matrix = typename Dune::TypeListElement<0, decltype(typeList)>::type;
79 using Domain = typename Dune::TypeListElement<1, decltype(typeList)>::type;
80 using Range = typename Dune::TypeListElement<2, decltype(typeList)>::type;
81 std::shared_ptr<Dune::Preconditioner<Domain, Range>> preconditioner
82 = std::make_shared<Preconditioner<Matrix, Domain, Range>>(matrix, config);
83 return preconditioner;
84 }
85 };
86
87 using IstlAmgPreconditionerFactory = Dune::AMGCreator;
88
89 template<class M, bool convert = false>
90 struct MatrixForSolver { using type = M; };
91
92 template<class M>
93 struct MatrixForSolver<M, true>
94 { using type = std::decay_t<decltype(MatrixConverter<M>::multiTypeToBCRSMatrix(std::declval<M>()))>; };
95
96 template<class V, bool convert = false>
97 struct VectorForSolver { using type = V; };
98
99 template<class V>
100 struct VectorForSolver<V, true>
101 { using type = std::decay_t<decltype(VectorConverter<V>::multiTypeToBlockVector(std::declval<V>()))>; };
102
103 template<class LSTraits, class LATraits, bool convert, bool parallel = LSTraits::canCommunicate>
104 struct MatrixOperator;
105
106 template<class LSTraits, class LATraits, bool convert>
107 struct MatrixOperator<LSTraits, LATraits, convert, true>
108 {
109 using M = typename MatrixForSolver<typename LATraits::Matrix, convert>::type;
110 using V = typename VectorForSolver<typename LATraits::Vector, convert>::type;
111 #if HAVE_MPI
112 using type = std::variant<
113 std::shared_ptr<typename LSTraits::template Sequential<M, V>::LinearOperator>,
114 std::shared_ptr<typename LSTraits::template ParallelOverlapping<M, V>::LinearOperator>,
115 std::shared_ptr<typename LSTraits::template ParallelNonoverlapping<M, V>::LinearOperator>
116 >;
117 #else
118 using type = std::variant<
119 std::shared_ptr<typename LSTraits::template Sequential<M, V>::LinearOperator>
120 >;
121 #endif
122 };
123
124 template<class LSTraits, class LATraits, bool convert>
125 struct MatrixOperator<LSTraits, LATraits, convert, false>
126 {
127 using M = typename MatrixForSolver<typename LATraits::Matrix, convert>::type;
128 using V = typename VectorForSolver<typename LATraits::Vector, convert>::type;
129 using type = std::variant<
130 std::shared_ptr<typename LSTraits::template Sequential<M, V>::LinearOperator>
131 >;
132 };
133
134 } // end namespace Dumux::Detail::IstlSolvers
135
136 namespace Dumux::Detail {
137
138 struct IstlSolverResult : public Dune::InverseOperatorResult
139 {
140 IstlSolverResult() = default;
141 IstlSolverResult(const IstlSolverResult&) = default;
142 IstlSolverResult(IstlSolverResult&&) = default;
143
144 IstlSolverResult(const Dune::InverseOperatorResult& o) : InverseOperatorResult(o) {}
145 61358 IstlSolverResult(Dune::InverseOperatorResult&& o) : InverseOperatorResult(std::move(o)) {}
146
147 61517 operator bool() const { return this->converged; }
148 };
149
150 /*!
151 * \ingroup Linear
152 * \brief Standard dune-istl iterative linear solvers
153 */
154 template<class LinearSolverTraits, class LinearAlgebraTraits,
155 class InverseOperator, class PreconditionerFactory,
156 bool convertMultiTypeLATypes = false>
157 class IstlIterativeLinearSolver
158 {
159 using Matrix = typename LinearAlgebraTraits::Matrix;
160 using XVector = typename LinearAlgebraTraits::Vector;
161 using BVector = typename LinearAlgebraTraits::Vector;
162 using Scalar = typename InverseOperator::real_type;
163
164 using ScalarProduct = Dune::ScalarProduct<typename InverseOperator::domain_type>;
165
166 static constexpr bool convertMultiTypeVectorAndMatrix
167 = convertMultiTypeLATypes && isMultiTypeBlockVector<XVector>::value;
168 using MatrixForSolver = typename Detail::IstlSolvers::MatrixForSolver<Matrix, convertMultiTypeVectorAndMatrix>::type;
169 using BVectorForSolver = typename Detail::IstlSolvers::VectorForSolver<BVector, convertMultiTypeVectorAndMatrix>::type;
170 using XVectorForSolver = typename Detail::IstlSolvers::VectorForSolver<XVector, convertMultiTypeVectorAndMatrix>::type;
171 // a variant type that can hold sequential, overlapping, and non-overlapping operators
172 using MatrixOperatorHolder = typename Detail::IstlSolvers::MatrixOperator<
173 LinearSolverTraits, LinearAlgebraTraits, convertMultiTypeVectorAndMatrix
174 >::type;
175
176 #if HAVE_MPI
177 using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>;
178 using ParallelHelper = ParallelISTLHelper<LinearSolverTraits>;
179 #endif
180
181 using ParameterInitializer = std::variant<std::string, Dune::ParameterTree>;
182 public:
183
184 /*!
185 * \brief Constructor for sequential solvers
186 */
187 95 IstlIterativeLinearSolver(const ParameterInitializer& params = "")
188 570 {
189
3/6
✓ Branch 1 taken 95 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 95 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 95 times.
 95 if (Dune::MPIHelper::getCommunication().size() > 1)
190 ✗ DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!");
191
192
1/2
✓ Branch 1 taken 95 times.
✗ Branch 2 not taken.
 95 initializeParameters_(params);
193 95 solverCategory_ = Dune::SolverCategory::sequential;
194
2/4
✓ Branch 1 taken 95 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 95 times.
 95 scalarProduct_ = std::make_shared<ScalarProduct>();
195 95 }
196
197 /*!
198 * \brief Constructor for parallel and sequential solvers
199 */
200 template <class GridView, class DofMapper>
201 262 IstlIterativeLinearSolver(const GridView& gridView,
202 const DofMapper& dofMapper,
203 const ParameterInitializer& params = "")
204 1572 {
205
1/2
✓ Branch 1 taken 258 times.
✗ Branch 2 not taken.
 262 initializeParameters_(params);
206 #if HAVE_MPI
207
2/3
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 93 times.
✗ Branch 2 not taken.
 262 solverCategory_ = Detail::solverCategory<LinearSolverTraits>(gridView);
208 if constexpr (LinearSolverTraits::canCommunicate)
209 {
210
211
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 213 times.
 239 if (solverCategory_ != Dune::SolverCategory::sequential)
212 {
213
3/6
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 22 times.
 22 parallelHelper_ = std::make_shared<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
214
9/10
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 8 times.
✓ Branch 6 taken 12 times.
✓ Branch 7 taken 10 times.
✓ Branch 8 taken 12 times.
✓ Branch 9 taken 2 times.
 24 communication_ = std::make_shared<Comm>(gridView.comm(), solverCategory_);
215
4/8
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 22 times.
 44 scalarProduct_ = Dune::createScalarProduct<XVector>(*communication_, solverCategory_);
216
3/6
✓ Branch 1 taken 22 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 22 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 22 times.
✗ Branch 8 not taken.
 66 parallelHelper_->createParallelIndexSet(*communication_);
217 }
218 else
219
2/4
✓ Branch 1 taken 213 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 213 times.
 217 scalarProduct_ = std::make_shared<ScalarProduct>();
220 }
221 else
222
2/4
✓ Branch 1 taken 23 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 23 times.
 23 scalarProduct_ = std::make_shared<ScalarProduct>();
223 #else
224 solverCategory_ = Dune::SolverCategory::sequential;
225 scalarProduct_ = std::make_shared<ScalarProduct>();
226 #endif
227 262 }
228
229 #if HAVE_MPI
230 /*!
231 * \brief Constructor with custom scalar product and communication
232 */
233 template <class GridView, class DofMapper>
234 IstlIterativeLinearSolver(std::shared_ptr<Comm> communication,
235 std::shared_ptr<ScalarProduct> scalarProduct,
236 const GridView& gridView,
237 const DofMapper& dofMapper,
238 const ParameterInitializer& params = "")
239 {
240 initializeParameters_(params);
241 solverCategory_ = Detail::solverCategory(gridView);
242 scalarProduct_ = scalarProduct;
243 communication_ = communication;
244 if constexpr (LinearSolverTraits::canCommunicate)
245 {
246 if (solverCategory_ != Dune::SolverCategory::sequential)
247 {
248 parallelHelper_ = std::make_shared<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper);
249 parallelHelper_->createParallelIndexSet(communication);
250 }
251 }
252 }
253 #endif
254
255 /*!
256 * \brief Solve the linear system Ax = b
257 */
258 IstlSolverResult solve(Matrix& A, XVector& x, BVector& b)
259 69631 { return solveSequentialOrParallel_(A, x, b); }
260
261 /*!
262 * \brief Set the matrix A of the linear system Ax = b for reuse
263 */
264 34 void setMatrix(std::shared_ptr<Matrix> A)
265 {
266
3/8
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 34 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 102 linearOperator_ = makeParallelOrSequentialLinearOperator_(std::move(A));
267
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
 34 solver_ = constructPreconditionedSolver_(linearOperator_);
268 34 }
269
270 /*!
271 * \brief Set the matrix A of the linear system Ax = b for reuse
272 * \note The client has to take care of the lifetime management of A
273 */
274 3 void setMatrix(Matrix& A)
275
2/6
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 3 { setMatrix(Dune::stackobject_to_shared_ptr(A)); }
276
277 /*!
278 * \brief Solve the linear system Ax = b where A has been set with \ref setMatrix
279 */
280 181 IstlSolverResult solve(XVector& x, BVector& b) const
281 {
282
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 169 times.
 181 if (!solver_)
283 ✗ DUNE_THROW(Dune::InvalidStateException, "Called solve(x, b) but no linear operator has been set");
284
285 181 return solveSequentialOrParallel_(x, b, *solver_);
286 }
287
288 /*!
289 * \brief Compute the 2-norm of vector x
290 */
291 8886 Scalar norm(const XVector& x) const
292 {
293 #if HAVE_MPI
294 if constexpr (LinearSolverTraits::canCommunicate)
295 {
296
2/2
✓ Branch 0 taken 112 times.
✓ Branch 1 taken 8492 times.
 8604 if (solverCategory_ == Dune::SolverCategory::nonoverlapping)
297 {
298
0/4
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 224 auto y(x); // make a copy because the vector needs to be made consistent
299 using GV = typename LinearSolverTraits::GridView;
300 using DM = typename LinearSolverTraits::DofMapper;
301
4/8
✓ Branch 1 taken 112 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 112 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 112 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 112 times.
✗ Branch 11 not taken.
 448 ParallelVectorHelper<GV, DM, LinearSolverTraits::dofCodim> vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper());
302
1/2
✓ Branch 1 taken 112 times.
✗ Branch 2 not taken.
 112 vectorHelper.makeNonOverlappingConsistent(y);
303
3/6
✓ Branch 1 taken 112 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 112 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 112 times.
✗ Branch 7 not taken.
 224 return scalarProduct_->norm(y);
304 }
305 }
306 #endif
307 if constexpr (convertMultiTypeVectorAndMatrix)
308 {
309
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 282 auto y = VectorConverter<XVector>::multiTypeToBlockVector(x);
310
3/6
✓ Branch 1 taken 282 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 282 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 282 times.
✗ Branch 7 not taken.
 846 return scalarProduct_->norm(y);
311 }
312 else
313
4/7
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 368 times.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 368 times.
✗ Branch 7 not taken.
 17722 return scalarProduct_->norm(x);
314 }
315
316 /*!
317 * \brief The name of the linear solver
318 */
319 const std::string& name() const
320 {
321
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 return name_;
322 }
323
324 /*!
325 * \brief Set the residual reduction tolerance
326 */
327 261 void setResidualReduction(double residReduction)
328 {
329
6/16
✓ Branch 2 taken 261 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 261 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 261 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 261 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 261 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 261 times.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 261 params_["reduction"] = std::to_string(residReduction);
330
331 // reconstruct the solver with new parameters
332
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 261 times.
 261 if (solver_)
333 ✗ solver_ = constructPreconditionedSolver_(linearOperator_);
334 261 }
335
336 /*!
337 * \brief Set the maximum number of linear solver iterations
338 */
339 void setMaxIter(std::size_t maxIter)
340 {
341 params_["maxit"] = std::to_string(maxIter);
342
343 // reconstruct the solver with new parameters
344 if (solver_)
345 solver_ = constructPreconditionedSolver_(linearOperator_);
346 }
347
348 /*!
349 * \brief Set the linear solver parameters
350 * \param params Either a std::string giving a parameter group (parameters are read from input file) or a Dune::ParameterTree
351 * \note In case of a Dune::ParameterTree, the parameters are passed trait to the linear solver and preconditioner
352 */
353 19 void setParams(const ParameterInitializer& params)
354 {
355 19 initializeParameters_(params);
356
357 // reconstruct the solver with new parameters
358
1/2
✓ Branch 0 taken 19 times.
✗ Branch 1 not taken.
 19 if (solver_)
359
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 19 times.
 19 solver_ = constructPreconditionedSolver_(linearOperator_);
360 19 }
361
362 private:
363
364 428 void initializeParameters_(const ParameterInitializer& params)
365 {
366
3/4
✓ Branch 0 taken 372 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 353 times.
✓ Branch 3 taken 19 times.
 428 if (std::holds_alternative<std::string>(params))
367 794 params_ = Dumux::LinearSolverParameters<LinearSolverTraits>::createParameterTree(std::get<std::string>(params));
368 else
369
1/2
✓ Branch 0 taken 19 times.
✗ Branch 1 not taken.
 62 params_ = std::get<Dune::ParameterTree>(params);
370 428 }
371
372 741 MatrixOperatorHolder makeSequentialLinearOperator_(std::shared_ptr<Matrix> A)
373 {
374 using SequentialTraits = typename LinearSolverTraits::template Sequential<MatrixForSolver, XVectorForSolver>;
375 if constexpr (convertMultiTypeVectorAndMatrix)
376 {
377 // create the BCRS matrix the IterativeSolver backend can handle
378
1/4
✓ Branch 3 taken 741 times.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 1482 auto M = std::make_shared<MatrixForSolver>(MatrixConverter<Matrix>::multiTypeToBCRSMatrix(*A));
379
4/8
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 741 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 741 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 741 times.
✗ Branch 8 not taken.
 1482 return std::make_shared<typename SequentialTraits::LinearOperator>(M);
380 }
381 else
382 {
383
8/9
✓ Branch 0 taken 668 times.
✓ Branch 1 taken 28139 times.
✓ Branch 2 taken 668 times.
✓ Branch 3 taken 28139 times.
✓ Branch 4 taken 75 times.
✓ Branch 5 taken 1 times.
✓ Branch 6 taken 75 times.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 57766 return std::make_shared<typename SequentialTraits::LinearOperator>(A);
384 }
385 }
386
387 template<class ParallelTraits>
388 ✗ MatrixOperatorHolder makeParallelLinearOperator_(std::shared_ptr<Matrix> A, ParallelTraits = {})
389 {
390 #if HAVE_MPI
391 // make matrix consistent
392
1/4
✓ Branch 1 taken 4582 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4582 prepareMatrixParallel<LinearSolverTraits, ParallelTraits>(*A, *parallelHelper_);
393
4/16
✓ Branch 1 taken 4582 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4582 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 4582 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 4582 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
 9164 return std::make_shared<typename ParallelTraits::LinearOperator>(std::move(A), *communication_);
394 #else
395 DUNE_THROW(Dune::InvalidStateException, "Calling makeParallelLinearOperator for sequential run");
396 #endif
397 }
398
399 MatrixOperatorHolder makeParallelOrSequentialLinearOperator_(std::shared_ptr<Matrix> A)
400 {
401 return executeSequentialOrParallel_(
402
2/6
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 32 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 64 [&]{ return makeSequentialLinearOperator_(std::move(A)); },
403
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 4 [&](auto traits){ return makeParallelLinearOperator_(std::move(A), traits); }
404
2/4
✓ Branch 1 taken 34 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
 70 );
405 }
406
407 30295 MatrixOperatorHolder makeSequentialLinearOperator_(Matrix& A)
408
2/6
✓ Branch 2 taken 29592 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 29592 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 60590 { return makeSequentialLinearOperator_(Dune::stackobject_to_shared_ptr<Matrix>(A)); }
409
410 MatrixOperatorHolder makeParallelOrSequentialLinearOperator_(Matrix& A)
411 { return makeParallelOrSequentialLinearOperator_(Dune::stackobject_to_shared_ptr<Matrix>(A)); }
412
413 template<class ParallelTraits>
414 9164 MatrixOperatorHolder makeParallelLinearOperator_(Matrix& A, ParallelTraits = {})
415
2/6
✓ Branch 2 taken 4582 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 4582 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 9164 { return makeParallelLinearOperator_<ParallelTraits>(Dune::stackobject_to_shared_ptr<Matrix>(A)); }
416
417 30295 IstlSolverResult solveSequential_(Matrix& A, XVector& x, BVector& b)
418 {
419 // construct solver from linear operator
420
1/2
✓ Branch 1 taken 25 times.
✗ Branch 2 not taken.
 30320 auto linearOperatorHolder = makeSequentialLinearOperator_(A);
421
2/6
✓ Branch 0 taken 29567 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 766 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 31061 auto solver = constructPreconditionedSolver_(linearOperatorHolder);
422
423
1/2
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
 31011 return solveSequential_(x, b, *solver);
424 }
425
426 741 IstlSolverResult solveSequential_(XVector& x, BVector& b, InverseOperator& solver) const
427 {
428
3/4
✓ Branch 1 taken 685 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28141 times.
✓ Branch 5 taken 25 times.
 29592 Dune::InverseOperatorResult result;
429 if constexpr (convertMultiTypeVectorAndMatrix)
430 {
431 // create the vector the IterativeSolver backend can handle
432
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 741 BVectorForSolver bTmp = VectorConverter<BVector>::multiTypeToBlockVector(b);
433
434 // create a block vector to which the linear solver writes the solution
435
3/8
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 741 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 741 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2223 XVectorForSolver y(bTmp.size());
436
437 // solve linear system
438
1/2
✓ Branch 1 taken 741 times.
✗ Branch 2 not taken.
 741 solver.apply(y, bTmp, result);
439
440 // copy back the result y into x
441
1/2
✓ Branch 0 taken 741 times.
✗ Branch 1 not taken.
 741 if (result.converged)
442 741 VectorConverter<XVector>::retrieveValues(x, y);
443 }
444 else
445 {
446 // solve linear system
447
3/4
✓ Branch 1 taken 685 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28141 times.
✓ Branch 5 taken 25 times.
 28920 solver.apply(x, b, result);
448 }
449
450
2/4
✓ Branch 0 taken 28775 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 51 times.
✗ Branch 3 not taken.
 29567 return result;
451 }
452
453 IstlSolverResult solveSequentialOrParallel_(Matrix& A, XVector& x, BVector& b)
454 {
455 return executeSequentialOrParallel_(
456 29592 [&]{ return solveSequential_(A, x, b); },
457 5236 [&](auto traits){ return solveParallel_(A, x, b, traits); }
458
4/6
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 times.
✓ Branch 5 taken 500 times.
✗ Branch 6 not taken.
 37729 );
459 }
460
461 IstlSolverResult solveSequentialOrParallel_(XVector& x, BVector& b, InverseOperator& solver) const
462 {
463 return executeSequentialOrParallel_(
464 243 [&]{ return solveSequential_(x, b, solver); },
465 100 [&](auto traits){ return solveParallel_(x, b, solver, traits); }
466 250 );
467 }
468
469 template<class ParallelTraits>
470 10472 IstlSolverResult solveParallel_(Matrix& A, XVector& x, BVector& b, ParallelTraits = {})
471 {
472 // construct solver from linear operator
473
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 10472 auto linearOperatorHolder = makeParallelLinearOperator_<ParallelTraits>(A);
474
2/6
✓ Branch 0 taken 5236 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 654 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 11780 auto solver = constructPreconditionedSolver_(linearOperatorHolder);
475
1/2
✓ Branch 1 taken 654 times.
✗ Branch 2 not taken.
 11780 return solveParallel_<ParallelTraits>(x, b, *solver);
476 }
477
478 template<class ParallelTraits>
479 ✗ IstlSolverResult solveParallel_(XVector& x, BVector& b, InverseOperator& solver, ParallelTraits = {}) const
480 {
481 #if HAVE_MPI
482 // make right hand side consistent
483
1/5
✓ Branch 1 taken 4582 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4582 prepareVectorParallel<LinearSolverTraits, ParallelTraits>(b, *parallelHelper_);
484
485 // solve linear system
486
1/5
✓ Branch 1 taken 4582 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4582 Dune::InverseOperatorResult result;
487
1/5
✓ Branch 1 taken 4582 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4582 solver.apply(x, b, result);
488
1/4
✓ Branch 0 taken 4582 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 4582 return result;
489 #else
490 DUNE_THROW(Dune::InvalidStateException, "Calling makeParallelLinearOperator for sequential run");
491 #endif
492 }
493
494
495 std::shared_ptr<InverseOperator> constructPreconditionedSolver_(MatrixOperatorHolder& ops)
496 {
497 133694 return std::visit([&](auto&& op)
498 {
499 using LinearOperator = typename std::decay_t<decltype(op)>::element_type;
500
20/60
✓ Branch 1 taken 3571 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3571 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3571 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3571 times.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✓ Branch 17 taken 4582 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 4582 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 4582 times.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 4582 times.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 33 taken 26539 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 26539 times.
✗ Branch 37 not taken.
✓ Branch 39 taken 26539 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✓ Branch 42 taken 26539 times.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✓ Branch 49 taken 36 times.
✗ Branch 50 not taken.
✓ Branch 52 taken 36 times.
✗ Branch 53 not taken.
✓ Branch 55 taken 36 times.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✓ Branch 58 taken 36 times.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 71 not taken.
✗ Branch 72 not taken.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✗ Branch 76 not taken.
✗ Branch 77 not taken.
✓ Branch 81 taken 153 times.
✗ Branch 82 not taken.
✓ Branch 84 taken 153 times.
✗ Branch 85 not taken.
✓ Branch 87 taken 153 times.
✗ Branch 88 not taken.
✗ Branch 89 not taken.
✓ Branch 90 taken 153 times.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
 69762 const auto& params = params_.sub("preconditioner");
501 using Prec = Dune::Preconditioner<typename LinearOperator::domain_type, typename LinearOperator::range_type>;
502 using TL = Dune::TypeList<typename LinearOperator::matrix_type, typename LinearOperator::domain_type, typename LinearOperator::range_type>;
503
0/12
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 18 not taken.
✗ Branch 19 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 69762 std::shared_ptr<Prec> prec = PreconditionerFactory{}(TL{}, op, params);
504
505 #if HAVE_MPI
506
33/96
✓ Branch 1 taken 3571 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3571 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3571 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 3571 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
✓ Branch 13 taken 3569 times.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✓ Branch 23 taken 4582 times.
✗ Branch 24 not taken.
✓ Branch 26 taken 4582 times.
✗ Branch 27 not taken.
✓ Branch 29 taken 4582 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 4582 times.
✗ Branch 33 not taken.
✓ Branch 34 taken 2464 times.
✓ Branch 35 taken 2118 times.
✓ Branch 37 taken 2464 times.
✗ Branch 38 not taken.
✓ Branch 40 taken 2464 times.
✗ Branch 41 not taken.
✓ Branch 42 taken 2464 times.
✗ Branch 43 not taken.
✓ Branch 45 taken 26539 times.
✗ Branch 46 not taken.
✓ Branch 48 taken 26539 times.
✗ Branch 49 not taken.
✓ Branch 51 taken 26539 times.
✗ Branch 52 not taken.
✓ Branch 54 taken 26539 times.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✓ Branch 57 taken 26539 times.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✓ Branch 67 taken 36 times.
✗ Branch 68 not taken.
✓ Branch 70 taken 36 times.
✗ Branch 71 not taken.
✓ Branch 73 taken 36 times.
✗ Branch 74 not taken.
✓ Branch 76 taken 36 times.
✗ Branch 77 not taken.
✗ Branch 78 not taken.
✓ Branch 79 taken 36 times.
✗ Branch 81 not taken.
✗ Branch 82 not taken.
✗ Branch 84 not taken.
✗ Branch 85 not taken.
✗ Branch 86 not taken.
✗ Branch 87 not taken.
✗ Branch 89 not taken.
✗ Branch 90 not taken.
✗ Branch 92 not taken.
✗ Branch 93 not taken.
✗ Branch 95 not taken.
✗ Branch 96 not taken.
✗ Branch 98 not taken.
✗ Branch 99 not taken.
✗ Branch 100 not taken.
✗ Branch 101 not taken.
✗ Branch 103 not taken.
✗ Branch 104 not taken.
✗ Branch 106 not taken.
✗ Branch 107 not taken.
✗ Branch 108 not taken.
✗ Branch 109 not taken.
✓ Branch 111 taken 153 times.
✗ Branch 112 not taken.
✓ Branch 114 taken 153 times.
✗ Branch 115 not taken.
✓ Branch 117 taken 153 times.
✗ Branch 118 not taken.
✓ Branch 120 taken 153 times.
✗ Branch 121 not taken.
✗ Branch 122 not taken.
✓ Branch 123 taken 153 times.
✗ Branch 125 not taken.
✗ Branch 126 not taken.
✗ Branch 128 not taken.
✗ Branch 129 not taken.
✗ Branch 130 not taken.
✗ Branch 131 not taken.
 69762 if (prec->category() != op->category() && prec->category() == Dune::SolverCategory::sequential)
507
6/32
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✓ Branch 8 taken 2464 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 2464 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 2464 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
 4932 prec = Dune::wrapPreconditioner4Parallel(prec, op);
508 #endif
509
10/24
✓ Branch 1 taken 3571 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3571 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 4582 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 4582 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 26539 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 26539 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 36 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 36 times.
✗ Branch 19 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✓ Branch 26 taken 153 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 153 times.
✗ Branch 29 not taken.
 69762 return std::make_shared<InverseOperator>(op, scalarProduct_, prec, params_);
510
13/26
✓ Branch 1 taken 6829 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 7483 times.
✓ Branch 5 taken 7 times.
✓ Branch 7 taken 1315 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 654 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 26468 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 26615 times.
✓ Branch 17 taken 25 times.
✓ Branch 19 taken 184 times.
✓ Branch 20 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 24 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 31 taken 51 times.
✗ Branch 32 not taken.
✓ Branch 34 taken 51 times.
✗ Branch 35 not taken.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
 34881 }, ops);
511 }
512
513 template<class Seq, class Par>
514 32080 decltype(auto) executeSequentialOrParallel_(Seq&& sequentialAction, Par&& parallelAction) const
515 {
516 #if HAVE_MPI
517 // For Dune::MultiTypeBlockMatrix there is currently no generic way
518 // of handling parallelism, we therefore can only solve these types of systems sequentially
519 if constexpr (isMultiTypeBlockMatrix<Matrix>::value || !LinearSolverTraits::canCommunicate)
520
1/5
✓ Branch 2 taken 34 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 2951 return sequentialAction();
521 else
522 {
523
3/4
✓ Branch 0 taken 25955 times.
✓ Branch 1 taken 654 times.
✓ Branch 2 taken 4582 times.
✗ Branch 3 not taken.
 32080 switch (solverCategory_)
524 {
525 26742 case Dune::SolverCategory::sequential:
526 26742 return sequentialAction();
527 2 case Dune::SolverCategory::nonoverlapping:
528 using NOTraits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, XVector>;
529 756 return parallelAction(NOTraits{});
530 ✗ case Dune::SolverCategory::overlapping:
531 using OTraits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, XVector>;
532 4582 return parallelAction(OTraits{});
533 ✗ default: DUNE_THROW(Dune::InvalidStateException, "Unknown solver category");
534 }
535 }
536 #else
537 return sequentialAction();
538 #endif
539 }
540
541 #if HAVE_MPI
542 std::shared_ptr<const ParallelHelper> parallelHelper_;
543 std::shared_ptr<Comm> communication_;
544 #endif
545
546 Dune::SolverCategory::Category solverCategory_;
547 std::shared_ptr<ScalarProduct> scalarProduct_;
548
549 // for stored solvers (reuse matrix)
550 MatrixOperatorHolder linearOperator_;
551 // for stored solvers (reuse matrix)
552 std::shared_ptr<InverseOperator> solver_;
553
554 Dune::ParameterTree params_;
555 std::string name_;
556 };
557
558 } // end namespace Dumux::Detail
559
560 namespace Dumux {
561
562 /*!
563 * \ingroup Linear
564 * \brief An ILU preconditioned BiCGSTAB solver using dune-istl
565 *
566 * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
567 * faster and smoother convergence than the original BiCG. It can be applied to
568 * nonsymmetric matrices.\n
569 * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
570 * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
571 * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
572 *
573 * Preconditioner: ILU(n) incomplete LU factorization. The order n indicates
574 * fill-in. It can be damped by the relaxation parameter
575 * LinearSolver.PreconditionerRelaxation.\n
576 * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
577 */
578 template<class LSTraits, class LATraits>
579 using ILUBiCGSTABIstlSolver =
580 Detail::IstlIterativeLinearSolver<LSTraits, LATraits,
581 Dune::BiCGSTABSolver<typename LATraits::SingleTypeVector>,
582 Detail::IstlSolvers::IstlDefaultBlockLevelPreconditionerFactory<Dune::SeqILU>,
583 // the Dune::ILU preconditioners don't accept multi-type matrices
584 /*convert multi-type istl types?*/ true
585 >;
586
587 /*!
588 * \ingroup Linear
589 * \brief An ILU preconditioned GMres solver using dune-istl
590 *
591 * Solver: The GMRes (generalized minimal residual) method is an iterative
592 * method for the numerical solution of a nonsymmetric system of linear
593 * equations.\n
594 * See: Saad, Y., Schultz, M. H. (1986). "GMRES: A generalized minimal residual
595 * algorithm for solving nonsymmetric linear systems." SIAM J. Sci. and Stat.
596 * Comput. 7: 856–869.
597 *
598 * Preconditioner: ILU(n) incomplete LU factorization. The order n indicates
599 * fill-in. It can be damped by the relaxation parameter
600 * LinearSolver.PreconditionerRelaxation.\n
601 * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
602 */
603 template<class LSTraits, class LATraits>
604 using ILURestartedGMResIstlSolver =
605 Detail::IstlIterativeLinearSolver<LSTraits, LATraits,
606 Dune::RestartedGMResSolver<typename LATraits::SingleTypeVector>,
607 Detail::IstlSolvers::IstlDefaultBlockLevelPreconditionerFactory<Dune::SeqILU>,
608 // the Dune::ILU preconditioners don't accept multi-type matrices
609 /*convert multi-type istl types?*/ true
610 >;
611
612 /*!
613 * \ingroup Linear
614 * \brief An SSOR-preconditioned BiCGSTAB solver using dune-istl
615 *
616 * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
617 * faster and smoother convergence than the original BiCG. While, it can be
618 * applied to nonsymmetric matrices, the preconditioner SSOR assumes symmetry.\n
619 * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
620 * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
621 * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
622 *
623 * Preconditioner: SSOR symmetric successive overrelaxation method. The
624 * relaxation is controlled by the parameter LinearSolver.PreconditionerRelaxation.
625 * In each preconditioning step, it is applied as often as given by the parameter
626 * LinearSolver.PreconditionerIterations.\n
627 * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
628 */
629 template<class LSTraits, class LATraits>
630 using SSORBiCGSTABIstlSolver =
631 Detail::IstlIterativeLinearSolver<LSTraits, LATraits,
632 Dune::BiCGSTABSolver<typename LATraits::Vector>,
633 Detail::IstlSolvers::IstlDefaultBlockLevelPreconditionerFactory<Dune::SeqSSOR>
634 >;
635
636 /*!
637 * \ingroup Linear
638 * \brief An SSOR-preconditioned CG solver using dune-istl
639 *
640 * Solver: CG (conjugate gradient) is an iterative method for solving linear
641 * systems with a symmetric, positive definite matrix.\n
642 * See: Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
643 * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
644 * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
645 *
646 * Preconditioner: SSOR symmetric successive overrelaxation method. The
647 * relaxation is controlled by the parameter LinearSolver.PreconditionerRelaxation.
648 * In each preconditioning step, it is applied as often as given by the parameter
649 * LinearSolver.PreconditionerIterations.\n
650 * See: Golub, G. H., and Van Loan, C. F. (2012). Matrix computations. JHU Press.
651 */
652 template<class LSTraits, class LATraits>
653 using SSORCGIstlSolver =
654 Detail::IstlIterativeLinearSolver<LSTraits, LATraits,
655 Dune::CGSolver<typename LATraits::Vector>,
656 Detail::IstlSolvers::IstlDefaultBlockLevelPreconditionerFactory<Dune::SeqSSOR>
657 >;
658
659 /*!
660 * \ingroup Linear
661 * \brief An AMG preconditioned BiCGSTAB solver using dune-istl
662 *
663 * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
664 * faster and smoother convergence than the original BiCG. While, it can be
665 * applied to nonsymmetric matrices, the preconditioner SSOR assumes symmetry.\n
666 * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
667 * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
668 * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
669 *
670 * Preconditioner: AMG (algebraic multigrid)
671 */
672 template<class LSTraits, class LATraits>
673 using AMGBiCGSTABIstlSolver =
674 Detail::IstlIterativeLinearSolver<LSTraits, LATraits,
675 Dune::BiCGSTABSolver<typename LATraits::SingleTypeVector>,
676 Detail::IstlSolvers::IstlAmgPreconditionerFactory,
677 // the AMG preconditioner doesn't accept multi-type matrices
678 /*convert multi-type istl types?*/ true
679 >;
680
681 /*!
682 * \ingroup Linear
683 * \brief An AMG preconditioned CG solver using dune-istl
684 *
685 * Solver: CG (conjugate gradient) is an iterative method for solving linear
686 * systems with a symmetric, positive definite matrix.\n
687 * See: Helfenstein, R., Koko, J. (2010). "Parallel preconditioned conjugate
688 * gradient algorithm on GPU", Journal of Computational and Applied Mathematics,
689 * Volume 236, Issue 15, Pages 3584–3590, http://dx.doi.org/10.1016/j.cam.2011.04.025.
690 *
691 * Preconditioner: AMG (algebraic multigrid)
692 */
693 template<class LSTraits, class LATraits>
694 using AMGCGIstlSolver =
695 Detail::IstlIterativeLinearSolver<LSTraits, LATraits,
696 Dune::CGSolver<typename LATraits::SingleTypeVector>,
697 Detail::IstlSolvers::IstlAmgPreconditionerFactory,
698 // the AMG preconditioner doesn't accept multi-type matrices
699 /*convert multi-type istl types?*/ true
700 >;
701
702 /*!
703 * \ingroup Linear
704 * \brief An Uzawa preconditioned BiCGSTAB solver using dune-istl
705 *
706 * Solver: The BiCGSTAB (stabilized biconjugate gradients method) solver has
707 * faster and smoother convergence than the original BiCG. While, it can be
708 * applied to nonsymmetric matrices, the preconditioner SSOR assumes symmetry.\n
709 * See: Van der Vorst, H. A. (1992). "Bi-CGSTAB: A Fast and Smoothly Converging
710 * Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems".
711 * SIAM J. Sci. and Stat. Comput. 13 (2): 631–644. doi:10.1137/0913035.
712 *
713 * Preconditioner: Uzawa method for saddle point problems
714 * \note Expects a 2x2 MultiTypeBlockMatrix
715 */
716 template<class LSTraits, class LATraits>
717 using UzawaBiCGSTABIstlSolver =
718 Detail::IstlIterativeLinearSolver<LSTraits, LATraits,
719 Dune::BiCGSTABSolver<typename LATraits::Vector>,
720 Detail::IstlSolvers::IstlDefaultBlockLevelPreconditionerFactory<Dumux::SeqUzawa>
721 >;
722
723 } // end namespace Dumux
724
725 namespace Dumux::Detail {
726
727 /*!
728 * \ingroup Linear
729 * \brief Direct dune-istl linear solvers
730 */
731 template<class LSTraits, class LATraits, template<class M> class Solver,
732 bool convertMultiTypeVectorAndMatrix = isMultiTypeBlockVector<typename LATraits::Vector>::value>
733 class DirectIstlSolver : public LinearSolver
734 {
735 using Matrix = typename LATraits::Matrix;
736 using XVector = typename LATraits::Vector;
737 using BVector = typename LATraits::Vector;
738
739 using MatrixForSolver = typename Detail::IstlSolvers::MatrixForSolver<Matrix, convertMultiTypeVectorAndMatrix>::type;
740 using BVectorForSolver = typename Detail::IstlSolvers::VectorForSolver<BVector, convertMultiTypeVectorAndMatrix>::type;
741 using XVectorForSolver = typename Detail::IstlSolvers::VectorForSolver<XVector, convertMultiTypeVectorAndMatrix>::type;
742 using InverseOperator = Dune::InverseOperator<XVectorForSolver, BVectorForSolver>;
743 public:
744 using LinearSolver::LinearSolver;
745
746 /*!
747 * \brief Solve the linear system Ax = b
748 */
749 IstlSolverResult solve(const Matrix& A, XVector& x, const BVector& b)
750 {
751 26340 return solve_(A, x, b);
752 }
753
754 /*!
755 * \brief Solve the linear system Ax = b using the matrix set with \ref setMatrix
756 */
757 800 IstlSolverResult solve(XVector& x, const BVector& b)
758 {
759
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 800 times.
 800 if (!solver_)
760 ✗ DUNE_THROW(Dune::InvalidStateException, "Called solve(x, b) but no linear operator has been set");
761
762 1600 return solve_(x, b, *solver_);
763 }
764
765 /*!
766 * \brief Set the matrix A of the linear system Ax = b for reuse
767 */
768 8 void setMatrix(std::shared_ptr<Matrix> A)
769 {
770 if constexpr (convertMultiTypeVectorAndMatrix)
771 matrix_ = std::make_shared<MatrixForSolver>(MatrixConverter<Matrix>::multiTypeToBCRSMatrix(A));
772 else
773 8 matrix_ = A;
774
775
2/4
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 8 times.
 16 solver_ = std::make_shared<Solver<MatrixForSolver>>(*matrix_);
776 8 }
777
778 /*!
779 * \brief Set the matrix A of the linear system Ax = b for reuse
780 * \note The client has to take care of the lifetime management of A
781 */
782 8 void setMatrix(Matrix& A)
783
2/6
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 8 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 16 { setMatrix(Dune::stackobject_to_shared_ptr(A)); }
784
785 /*!
786 * \brief name of the linear solver
787 */
788 std::string name() const
789 {
790 return "Direct solver";
791 }
792
793 private:
794 26411 IstlSolverResult solve_(const Matrix& A, XVector& x, const BVector& b)
795 {
796 // support dune-istl multi-type block vector/matrix by copying
797 if constexpr (convertMultiTypeVectorAndMatrix)
798 {
799 10258 const auto AA = MatrixConverter<Matrix>::multiTypeToBCRSMatrix(A);
800
2/4
✓ Branch 1 taken 5129 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5129 times.
✗ Branch 5 not taken.
 15387 Solver<MatrixForSolver> solver(AA, this->verbosity() > 0);
801
1/2
✓ Branch 1 taken 5129 times.
✗ Branch 2 not taken.
 10258 return solve_(x, b, solver);
802 }
803 else
804 {
805 63846 Solver<MatrixForSolver> solver(A, this->verbosity() > 0);
806
1/2
✓ Branch 1 taken 21211 times.
✗ Branch 2 not taken.
 42564 return solve_(x, b, solver);
807 }
808 }
809
810 27211 IstlSolverResult solve_(XVector& x, const BVector& b, InverseOperator& solver) const
811 {
812 27211 Dune::InverseOperatorResult result;
813
814 if constexpr (convertMultiTypeVectorAndMatrix)
815 {
816
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 5129 auto bb = VectorConverter<BVector>::multiTypeToBlockVector(b);
817
2/6
✓ Branch 1 taken 5129 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5129 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 10258 XVectorForSolver xx(bb.size());
818
1/2
✓ Branch 1 taken 5129 times.
✗ Branch 2 not taken.
 5129 solver.apply(xx, bb, result);
819 10258 checkResult_(xx, result);
820 5129 if (result.converged)
821 5125 VectorConverter<XVector>::retrieveValues(x, xx);
822
1/2
✓ Branch 0 taken 5129 times.
✗ Branch 1 not taken.
 10258 return result;
823 }
824 else
825 {
826
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 22082 BVectorForSolver bTmp(b);
827
1/2
✓ Branch 1 taken 22011 times.
✗ Branch 2 not taken.
 22082 solver.apply(x, bTmp, result);
828 44164 checkResult_(x, result);
829 22082 return result;
830 }
831 }
832
833
834 ✗ void checkResult_(XVectorForSolver& x, Dune::InverseOperatorResult& result) const
835 {
836
4/5
✓ Branch 1 taken 51 times.
✓ Branch 2 taken 27065 times.
✓ Branch 3 taken 4 times.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 27140 flatVectorForEach(x, [&](auto&& entry, std::size_t){
837 using std::isnan, std::isinf;
838
23/84
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 2328066 times.
✓ Branch 13 taken 175 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2328066 times.
✗ Branch 16 not taken.
✓ Branch 17 taken 2328066 times.
✓ Branch 18 taken 2328109 times.
✓ Branch 19 taken 132 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 2328109 times.
✗ Branch 22 not taken.
✓ Branch 23 taken 2328109 times.
✓ Branch 24 taken 62245 times.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✓ Branch 27 taken 62245 times.
✗ Branch 28 not taken.
✓ Branch 29 taken 62245 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 23127806 times.
✓ Branch 37 taken 5785 times.
✓ Branch 38 taken 240 times.
✓ Branch 39 taken 23127566 times.
✓ Branch 40 taken 240 times.
✓ Branch 41 taken 23127566 times.
✗ Branch 42 not taken.
✗ Branch 43 not taken.
✗ Branch 44 not taken.
✗ Branch 45 not taken.
✗ Branch 46 not taken.
✗ Branch 47 not taken.
✗ Branch 48 not taken.
✗ Branch 49 not taken.
✗ Branch 50 not taken.
✗ Branch 51 not taken.
✗ Branch 52 not taken.
✗ Branch 53 not taken.
✗ Branch 54 not taken.
✗ Branch 55 not taken.
✗ Branch 56 not taken.
✗ Branch 57 not taken.
✗ Branch 58 not taken.
✗ Branch 59 not taken.
✗ Branch 60 not taken.
✗ Branch 61 not taken.
✗ Branch 62 not taken.
✗ Branch 63 not taken.
✗ Branch 64 not taken.
✗ Branch 65 not taken.
✗ Branch 66 not taken.
✗ Branch 67 not taken.
✗ Branch 68 not taken.
✗ Branch 69 not taken.
✗ Branch 70 not taken.
✗ Branch 71 not taken.
✓ Branch 72 taken 376839 times.
✗ Branch 73 not taken.
✗ Branch 74 not taken.
✓ Branch 75 taken 376839 times.
✗ Branch 76 not taken.
✓ Branch 77 taken 376839 times.
✓ Branch 78 taken 147780 times.
✗ Branch 79 not taken.
✗ Branch 80 not taken.
✓ Branch 81 taken 147780 times.
✗ Branch 82 not taken.
✓ Branch 83 taken 147780 times.
 25986626 if (isnan(entry) || isinf(entry))
839 6332 result.converged = false;
840 });
841 ✗ }
842
843 //! matrix when using the setMatrix interface for matrix reuse
844 std::shared_ptr<MatrixForSolver> matrix_;
845 //! solver when using the setMatrix interface for matrix reuse
846 std::shared_ptr<InverseOperator> solver_;
847 };
848
849 } // end namespace Dumux::Detail
850
851 #if HAVE_SUPERLU
852 #include <dune/istl/superlu.hh>
853
854 namespace Dumux {
855
856 /*!
857 * \ingroup Linear
858 * \brief Direct linear solver using the SuperLU library.
859 *
860 * See: Li, X. S. (2005). "An overview of SuperLU: Algorithms, implementation,
861 * and user interface." ACM Transactions on Mathematical Software (TOMS) 31(3): 302-325.
862 * http://crd-legacy.lbl.gov/~xiaoye/SuperLU/
863 */
864 template<class LSTraits, class LATraits>
865 using SuperLUIstlSolver = Detail::DirectIstlSolver<LSTraits, LATraits, Dune::SuperLU>;
866
867 } // end namespace Dumux
868
869 #endif // HAVE_SUPERLU
870
871 #if HAVE_UMFPACK
872 #include <dune/istl/umfpack.hh>
873 #include <dune/common/version.hh>
874
875 namespace Dumux {
876
877 /*!
878 * \ingroup Linear
879 * \brief Direct linear solver using the UMFPack library.
880 *
881 * See: Davis, Timothy A. (2004). "Algorithm 832". ACM Transactions on
882 * Mathematical Software 30 (2): 196–199. doi:10.1145/992200.992206.
883 * http://faculty.cse.tamu.edu/davis/suitesparse.html
884 */
885 template<class LSTraits, class LATraits>
886 using UMFPackIstlSolver = Detail::DirectIstlSolver<
887 LSTraits, LATraits, Dune::UMFPack
888 #if DUNE_VERSION_GTE(DUNE_ISTL,2,10)
889 , false // no need to convert multi-type matrix anymore
890 #endif
891 >;
892
893 } // end namespace Dumux
894
895 #endif // HAVE_UMFPACK
896
897 #endif
898