GCC Code Coverage Report

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