| 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 Provides a generic linear solver based on the ISTL that chooses the | ||
| 11 | * solver and preconditioner at runtime | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_LINEAR_ISTL_SOLVERFACTORYBACKEND_HH | ||
| 15 | #define DUMUX_LINEAR_ISTL_SOLVERFACTORYBACKEND_HH | ||
| 16 | |||
| 17 | #include <memory> | ||
| 18 | |||
| 19 | #include <dune/common/parallel/mpihelper.hh> | ||
| 20 | #include <dune/common/parametertree.hh> | ||
| 21 | |||
| 22 | #include "linearsolverparameters.hh" | ||
| 23 | |||
| 24 | // preconditioners | ||
| 25 | #include <dune/istl/preconditioners.hh> | ||
| 26 | #include <dune/istl/paamg/amg.hh> | ||
| 27 | #include "preconditioners.hh" | ||
| 28 | |||
| 29 | // solvers | ||
| 30 | #include <dune/istl/solvers.hh> | ||
| 31 | #include <dune/istl/solverfactory.hh> | ||
| 32 | |||
| 33 | #include <dumux/common/typetraits/matrix.hh> | ||
| 34 | #include <dumux/common/typetraits/vector.hh> | ||
| 35 | #include <dumux/linear/solver.hh> | ||
| 36 | #include <dumux/linear/parallelhelpers.hh> | ||
| 37 | #include <dumux/linear/istlsolverregistry.hh> | ||
| 38 | #include <dumux/linear/solvercategory.hh> | ||
| 39 | |||
| 40 | namespace Dumux { | ||
| 41 | |||
| 42 | // initSolverFactoriesForMultiTypeBlockMatrix differs in different compilation units, | ||
| 43 | // so we have it in an anonymous namespace | ||
| 44 | namespace { | ||
| 45 | |||
| 46 | /*! | ||
| 47 | * \brief Initializes the direct solvers, preconditioners and iterative solvers in | ||
| 48 | * the factories with the corresponding Matrix and Vector types. | ||
| 49 | * \note We currently consider only direct solvers and preconditioners provided by | ||
| 50 | * Dumux which, unlike the ones implemented in Dune, also support MultiTypeBlockMatrices. | ||
| 51 | * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices. | ||
| 52 | * \tparam LinearOperator the linear operator type | ||
| 53 | */ | ||
| 54 | template<class LinearOperator> | ||
| 55 | 1 | int initSolverFactoriesForMultiTypeBlockMatrix() | |
| 56 | { | ||
| 57 | using M = typename LinearOperator::matrix_type; | ||
| 58 | using X = typename LinearOperator::range_type; | ||
| 59 | using Y = typename LinearOperator::domain_type; | ||
| 60 | using TL = Dune::TypeList<M,X,Y>; | ||
| 61 | 1 | auto& dsfac = Dune::DirectSolverFactory<M,X,Y>::instance(); | |
| 62 | 1 | Dune::addRegistryToFactory<TL>(dsfac, Dumux::MultiTypeBlockMatrixDirectSolverTag{}); | |
| 63 | 1 | auto& pfac = Dune::PreconditionerFactory<LinearOperator,X,Y>::instance(); | |
| 64 | 1 | Dune::addRegistryToFactory<TL>(pfac, Dumux::MultiTypeBlockMatrixPreconditionerTag{}); | |
| 65 | using TLS = Dune::TypeList<X,Y>; | ||
| 66 | 1 | auto& isfac = Dune::IterativeSolverFactory<X,Y>::instance(); | |
| 67 | 1 | return Dune::addRegistryToFactory<TLS>(isfac, Dune::IterativeSolverTag{}); | |
| 68 | } | ||
| 69 | } // end namespace | ||
| 70 | |||
| 71 | /*! | ||
| 72 | * \brief Initialize the solver factories for regular matrices or MultiTypeBlockMatrices | ||
| 73 | * \tparam Matrix the matrix | ||
| 74 | * \tparam LinearOperator the linear operator | ||
| 75 | * | ||
| 76 | * \note This function could be removed once Dune::initSolverFactories supports MultiTypeBlockMatrices. | ||
| 77 | */ | ||
| 78 | template<class Matrix, class LinearOperator> | ||
| 79 | void initSolverFactories() | ||
| 80 | { | ||
| 81 | if constexpr (isMultiTypeBlockMatrix<Matrix>::value) | ||
| 82 | 1 | initSolverFactoriesForMultiTypeBlockMatrix<LinearOperator>(); | |
| 83 | else | ||
| 84 | 44 | Dune::initSolverFactories<LinearOperator>(); | |
| 85 | } | ||
| 86 | |||
| 87 | /*! | ||
| 88 | * \ingroup Linear | ||
| 89 | * \brief A linear solver using the dune-istl solver factory | ||
| 90 | * to choose the solver and preconditioner at runtime. | ||
| 91 | * \note the solvers are configured via the input file | ||
| 92 | */ | ||
| 93 | namespace Detail { | ||
| 94 | |||
| 95 | template <class LinearSolverTraits> | ||
| 96 | class OldIstlSolverFactoryBackend : public LinearSolver | ||
| 97 | { | ||
| 98 | public: | ||
| 99 | |||
| 100 | /*! | ||
| 101 | * \brief Update the solver after grid adaption | ||
| 102 | * | ||
| 103 | * \param gridView the grid view on which we are performing the multi-grid | ||
| 104 | * \param dofMapper an index mapper for dof entities | ||
| 105 | */ | ||
| 106 | void updateAfterGridAdaption(const typename LinearSolverTraits::GridView& gridView, | ||
| 107 | const typename LinearSolverTraits::DofMapper& dofMapper) | ||
| 108 | { | ||
| 109 | #if HAVE_MPI | ||
| 110 | if (isParallel_) | ||
| 111 | parallelHelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper); | ||
| 112 | #endif | ||
| 113 | } | ||
| 114 | |||
| 115 | /*! | ||
| 116 | * \brief Solve a linear system. | ||
| 117 | * | ||
| 118 | * \param A the matrix | ||
| 119 | * \param x the seeked solution vector, containing the initial solution upon entry | ||
| 120 | * \param b the right hand side vector | ||
| 121 | */ | ||
| 122 | template<class Matrix, class Vector> | ||
| 123 | bool solve(Matrix& A, Vector& x, Vector& b) | ||
| 124 | { | ||
| 125 | #if HAVE_MPI | ||
| 126 | solveSequentialOrParallel_(A, x, b); | ||
| 127 | #else | ||
| 128 | solveSequential_(A, x, b); | ||
| 129 | #endif | ||
| 130 | firstCall_ = false; | ||
| 131 | return result_.converged; | ||
| 132 | } | ||
| 133 | |||
| 134 | const Dune::InverseOperatorResult& result() const | ||
| 135 | { | ||
| 136 | return result_; | ||
| 137 | } | ||
| 138 | |||
| 139 | const std::string& name() const | ||
| 140 | { | ||
| 141 | return name_; | ||
| 142 | } | ||
| 143 | |||
| 144 | private: | ||
| 145 | |||
| 146 | void initializeParameters_() | ||
| 147 | { | ||
| 148 | params_ = Dumux::LinearSolverParameters<LinearSolverTraits>::createParameterTree(paramGroup_); | ||
| 149 | checkMandatoryParameters_(); | ||
| 150 | name_ = params_.get<std::string>("preconditioner.type") + "-preconditioned " + params_.get<std::string>("type"); | ||
| 151 | if (params_.get<int>("verbose", 0) > 0) | ||
| 152 | std::cout << "Initialized linear solver of type: " << name_ << std::endl; | ||
| 153 | } | ||
| 154 | |||
| 155 | void checkMandatoryParameters_() | ||
| 156 | { | ||
| 157 | if (!params_.hasKey("type")) | ||
| 158 | DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Type\" parameter to select the solver"); | ||
| 159 | |||
| 160 | if (!params_.hasKey("preconditioner.type")) | ||
| 161 | DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Preconditioner.Type\" parameter to select the preconditioner"); | ||
| 162 | } | ||
| 163 | |||
| 164 | #if HAVE_MPI | ||
| 165 | template<class Matrix, class Vector> | ||
| 166 | void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b) | ||
| 167 | { | ||
| 168 | // Dune::MultiTypeBlockMatrix does not provide a ColIterator which is needed by Dune::NonoverlappingSchwarzOperator. | ||
| 169 | // We therefore can only solve these types of systems sequentially. | ||
| 170 | // TODO: This can be adapted once the situation in Dune ISTL changes. | ||
| 171 | if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockMatrix<Matrix>::value) | ||
| 172 | { | ||
| 173 | if (isParallel_) | ||
| 174 | { | ||
| 175 | if (LinearSolverTraits::isNonOverlapping(parallelHelper_->gridView())) | ||
| 176 | { | ||
| 177 | using PTraits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>; | ||
| 178 | solveParallel_<PTraits>(A, x, b); | ||
| 179 | } | ||
| 180 | else | ||
| 181 | { | ||
| 182 | using PTraits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>; | ||
| 183 | solveParallel_<PTraits>(A, x, b); | ||
| 184 | } | ||
| 185 | } | ||
| 186 | else | ||
| 187 | solveSequential_(A, x, b); | ||
| 188 | } | ||
| 189 | else | ||
| 190 | { | ||
| 191 | solveSequential_(A, x, b); | ||
| 192 | } | ||
| 193 | } | ||
| 194 | |||
| 195 | template<class ParallelTraits, class Matrix, class Vector> | ||
| 196 | void solveParallel_(Matrix& A, Vector& x, Vector& b) | ||
| 197 | { | ||
| 198 | using Comm = typename ParallelTraits::Comm; | ||
| 199 | using LinearOperator = typename ParallelTraits::LinearOperator; | ||
| 200 | using ScalarProduct = typename ParallelTraits::ScalarProduct; | ||
| 201 | |||
| 202 | if (firstCall_) | ||
| 203 | initSolverFactories<Matrix, LinearOperator>(); | ||
| 204 | |||
| 205 | std::shared_ptr<Comm> comm; | ||
| 206 | std::shared_ptr<LinearOperator> linearOperator; | ||
| 207 | std::shared_ptr<ScalarProduct> scalarProduct; | ||
| 208 | prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, comm, linearOperator, scalarProduct, *parallelHelper_); | ||
| 209 | |||
| 210 | // construct solver | ||
| 211 | auto solver = getSolverFromFactory_(linearOperator); | ||
| 212 | |||
| 213 | // solve linear system | ||
| 214 | solver->apply(x, b, result_); | ||
| 215 | } | ||
| 216 | #endif // HAVE_MPI | ||
| 217 | |||
| 218 | template<class Matrix, class Vector> | ||
| 219 | void solveSequential_(Matrix& A, Vector& x, Vector& b) | ||
| 220 | { | ||
| 221 | // construct linear operator | ||
| 222 | using Traits = typename LinearSolverTraits::template Sequential<Matrix, Vector>; | ||
| 223 | using LinearOperator = typename Traits::LinearOperator; | ||
| 224 | auto linearOperator = std::make_shared<LinearOperator>(A); | ||
| 225 | |||
| 226 | if (firstCall_) | ||
| 227 | initSolverFactories<Matrix, LinearOperator>(); | ||
| 228 | |||
| 229 | // construct solver | ||
| 230 | auto solver = getSolverFromFactory_(linearOperator); | ||
| 231 | |||
| 232 | // solve linear system | ||
| 233 | solver->apply(x, b, result_); | ||
| 234 | } | ||
| 235 | |||
| 236 | template<class LinearOperator> | ||
| 237 | auto getSolverFromFactory_(std::shared_ptr<LinearOperator>& fop) | ||
| 238 | { | ||
| 239 | try { return Dune::getSolverFromFactory(fop, params_); } | ||
| 240 | catch(Dune::Exception& e) | ||
| 241 | { | ||
| 242 | std::cerr << "Could not create solver with factory" << std::endl; | ||
| 243 | std::cerr << e.what() << std::endl; | ||
| 244 | throw e; | ||
| 245 | } | ||
| 246 | } | ||
| 247 | |||
| 248 | const std::string paramGroup_; | ||
| 249 | #if HAVE_MPI | ||
| 250 | std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> parallelHelper_; | ||
| 251 | #endif | ||
| 252 | bool isParallel_ = false; | ||
| 253 | bool firstCall_; | ||
| 254 | |||
| 255 | Dune::InverseOperatorResult result_; | ||
| 256 | Dune::ParameterTree params_; | ||
| 257 | std::string name_; | ||
| 258 | }; | ||
| 259 | |||
| 260 | template <class LinearSolverTraits, class LinearAlgebraTraits> | ||
| 261 | class IstlSolverFactoryBackend : public LinearSolver | ||
| 262 | { | ||
| 263 | using Matrix = typename LinearAlgebraTraits::Matrix; | ||
| 264 | using Vector = typename LinearAlgebraTraits::Vector; | ||
| 265 | using ScalarProduct = Dune::ScalarProduct<Vector>; | ||
| 266 | #if HAVE_MPI | ||
| 267 | using Comm = Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<96>, int>; | ||
| 268 | #endif | ||
| 269 | public: | ||
| 270 | |||
| 271 | /*! | ||
| 272 | * \brief Construct the backend for the sequential case only | ||
| 273 | * | ||
| 274 | * \param paramGroup the parameter group for parameter lookup | ||
| 275 | */ | ||
| 276 | 3 | IstlSolverFactoryBackend(const std::string& paramGroup = "") | |
| 277 | : paramGroup_(paramGroup) | ||
| 278 |
14/42✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 3 times.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 3 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 3 times.
✗ Branch 19 not taken.
✓ Branch 21 taken 3 times.
✗ Branch 22 not taken.
✓ Branch 24 taken 3 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 3 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 3 times.
✗ Branch 34 not taken.
✓ Branch 36 taken 3 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✓ Branch 39 taken 3 times.
✗ Branch 40 not taken.
✗ Branch 41 not taken.
✗ Branch 44 not taken.
✗ Branch 45 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.
|
6 | , isParallel_(Dune::MPIHelper::getCommunication().size() > 1) |
| 279 | { | ||
| 280 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (isParallel_) |
| 281 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Using sequential constructor for parallel run. Use signature with gridView and dofMapper!"); | |
| 282 | |||
| 283 | 3 | firstCall_ = true; | |
| 284 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | initializeParameters_(); |
| 285 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
|
3 | scalarProduct_ = std::make_shared<ScalarProduct>(); |
| 286 | 3 | solverCategory_ = Dune::SolverCategory::sequential; | |
| 287 | 3 | } | |
| 288 | |||
| 289 | /*! | ||
| 290 | * \brief Construct the backend for parallel or sequential runs | ||
| 291 | * | ||
| 292 | * \param gridView the grid view for parallel communication via the grid | ||
| 293 | * \param dofMapper an index mapper for dof entities | ||
| 294 | * \param paramGroup the parameter group for parameter lookup | ||
| 295 | */ | ||
| 296 | 42 | IstlSolverFactoryBackend(const typename LinearSolverTraits::GridView& gridView, | |
| 297 | const typename LinearSolverTraits::DofMapper& dofMapper, | ||
| 298 | const std::string& paramGroup = "") | ||
| 299 | : paramGroup_(paramGroup) | ||
| 300 | #if HAVE_MPI | ||
| 301 |
15/45✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 42 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 42 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 42 times.
✓ Branch 12 taken 42 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 42 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 26 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 26 times.
✓ Branch 21 taken 16 times.
✓ Branch 22 taken 26 times.
✓ Branch 24 taken 42 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 42 times.
✗ Branch 28 not taken.
✓ Branch 30 taken 42 times.
✗ Branch 31 not taken.
✓ Branch 33 taken 42 times.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 36 taken 26 times.
✗ Branch 37 not taken.
✗ Branch 38 not taken.
✗ Branch 39 not taken.
✗ Branch 40 not taken.
✗ 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.
|
84 | , isParallel_(gridView.comm().size() > 1) |
| 302 | #endif | ||
| 303 | { | ||
| 304 | 42 | firstCall_ = true; | |
| 305 |
1/2✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
|
42 | initializeParameters_(); |
| 306 | #if HAVE_MPI | ||
| 307 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
42 | solverCategory_ = Detail::solverCategory<LinearSolverTraits>(gridView); |
| 308 |
2/2✓ Branch 0 taken 30 times.
✓ Branch 1 taken 12 times.
|
42 | if (solverCategory_ != Dune::SolverCategory::sequential) |
| 309 | { | ||
| 310 |
3/6✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 30 times.
|
30 | parallelHelper_ = std::make_unique<ParallelISTLHelper<LinearSolverTraits>>(gridView, dofMapper); |
| 311 |
7/9✗ Branch 0 not taken.
✓ Branch 1 taken 18 times.
✓ Branch 2 taken 12 times.
✓ Branch 3 taken 18 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 12 times.
✓ Branch 6 taken 18 times.
✓ Branch 7 taken 12 times.
✓ Branch 8 taken 18 times.
|
30 | comm_ = std::make_shared<Comm>(gridView.comm(), solverCategory_); |
| 312 |
4/8✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 30 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 30 times.
|
60 | scalarProduct_ = Dune::createScalarProduct<Vector>(*comm_, solverCategory_); |
| 313 |
3/6✓ Branch 1 taken 30 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 30 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 30 times.
✗ Branch 8 not taken.
|
90 | parallelHelper_->createParallelIndexSet(*comm_); |
| 314 | } | ||
| 315 | else | ||
| 316 |
2/4✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 12 times.
|
12 | scalarProduct_ = std::make_shared<ScalarProduct>(); |
| 317 | #else | ||
| 318 | solverCategory_ = Dune::SolverCategory::sequential; | ||
| 319 | scalarProduct_ = std::make_shared<ScalarProduct>(); | ||
| 320 | #endif | ||
| 321 | 42 | } | |
| 322 | |||
| 323 | /*! | ||
| 324 | * \brief Solve a linear system. | ||
| 325 | * | ||
| 326 | * \param A the matrix | ||
| 327 | * \param x the seeked solution vector, containing the initial solution upon entry | ||
| 328 | * \param b the right hand side vector | ||
| 329 | */ | ||
| 330 | bool solve(Matrix& A, Vector& x, Vector& b) | ||
| 331 | { | ||
| 332 | #if HAVE_MPI | ||
| 333 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
1113 | solveSequentialOrParallel_(A, x, b); |
| 334 | #else | ||
| 335 | solveSequential_(A, x, b); | ||
| 336 | #endif | ||
| 337 | 1062 | firstCall_ = false; | |
| 338 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
1062 | return result_.converged; |
| 339 | } | ||
| 340 | |||
| 341 | const Dune::InverseOperatorResult& result() const | ||
| 342 | { | ||
| 343 | return result_; | ||
| 344 | } | ||
| 345 | |||
| 346 | const std::string& name() const | ||
| 347 | { | ||
| 348 | return name_; | ||
| 349 | } | ||
| 350 | |||
| 351 | 186 | Scalar norm(const Vector& x) const | |
| 352 | { | ||
| 353 | #if HAVE_MPI | ||
| 354 | if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockVector<Vector>::value) | ||
| 355 | { | ||
| 356 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 136 times.
|
186 | if (solverCategory_ == Dune::SolverCategory::nonoverlapping) |
| 357 | { | ||
| 358 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
100 | auto y(x); // make a copy because the vector needs to be made consistent |
| 359 | using GV = typename LinearSolverTraits::GridView; | ||
| 360 | using DM = typename LinearSolverTraits::DofMapper; | ||
| 361 |
4/8✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 50 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 50 times.
✗ Branch 11 not taken.
|
200 | ParallelVectorHelper<GV, DM, LinearSolverTraits::dofCodim> vectorHelper(parallelHelper_->gridView(), parallelHelper_->dofMapper()); |
| 362 |
1/2✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
|
50 | vectorHelper.makeNonOverlappingConsistent(y); |
| 363 |
3/6✓ Branch 1 taken 50 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
|
100 | return scalarProduct_->norm(y); |
| 364 | } | ||
| 365 | } | ||
| 366 | #endif | ||
| 367 |
2/4✓ Branch 3 taken 17 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 17 times.
✗ Branch 7 not taken.
|
408 | return scalarProduct_->norm(x); |
| 368 | } | ||
| 369 | |||
| 370 | private: | ||
| 371 | |||
| 372 | 45 | void initializeParameters_() | |
| 373 | { | ||
| 374 | 45 | params_ = Dumux::LinearSolverParameters<LinearSolverTraits>::createParameterTree(paramGroup_); | |
| 375 | 45 | checkMandatoryParameters_(); | |
| 376 |
18/42✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 45 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 45 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 45 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 45 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 45 times.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✓ Branch 26 taken 45 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 45 times.
✗ Branch 29 not taken.
✓ Branch 30 taken 45 times.
✓ Branch 31 taken 45 times.
✗ Branch 32 not taken.
✓ Branch 33 taken 1 times.
✓ Branch 34 taken 44 times.
✓ Branch 35 taken 1 times.
✓ Branch 36 taken 44 times.
✗ Branch 37 not taken.
✓ Branch 38 taken 45 times.
✓ Branch 40 taken 45 times.
✗ Branch 41 not taken.
✗ 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.
|
181 | name_ = params_.get<std::string>("preconditioner.type") + "-preconditioned " + params_.get<std::string>("type"); |
| 377 |
8/14✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 45 times.
✓ Branch 11 taken 3 times.
✓ Branch 12 taken 42 times.
✓ Branch 13 taken 3 times.
✓ Branch 14 taken 42 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
90 | if (params_.get<int>("verbose", 0) > 0) |
| 378 | 6 | std::cout << "Initialized linear solver of type: " << name_ << std::endl; | |
| 379 | 45 | } | |
| 380 | |||
| 381 | 45 | void checkMandatoryParameters_() | |
| 382 | { | ||
| 383 |
6/14✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 45 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 45 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 45 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
90 | if (!params_.hasKey("type")) |
| 384 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Type\" parameter to select the solver"); | |
| 385 | |||
| 386 |
6/14✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 45 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 45 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 45 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
|
180 | if (!params_.hasKey("preconditioner.type")) |
| 387 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Solver factory needs \"LinearSolver.Preconditioner.Type\" parameter to select the preconditioner"); | |
| 388 | 45 | } | |
| 389 | |||
| 390 | #if HAVE_MPI | ||
| 391 | 1011 | void solveSequentialOrParallel_(Matrix& A, Vector& x, Vector& b) | |
| 392 | { | ||
| 393 | // Dune::MultiTypeBlockMatrix does not provide a ColIterator which is needed by Dune::NonoverlappingSchwarzOperator. | ||
| 394 | // We therefore can only solve these types of systems sequentially. | ||
| 395 | // TODO: This can be adapted once the situation in Dune ISTL changes. | ||
| 396 | if constexpr (LinearSolverTraits::canCommunicate && !isMultiTypeBlockMatrix<Matrix>::value) | ||
| 397 | { | ||
| 398 |
2/2✓ Branch 0 taken 724 times.
✓ Branch 1 taken 287 times.
|
1011 | if (isParallel_) |
| 399 | { | ||
| 400 |
4/6✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 196 times.
✓ Branch 4 taken 98 times.
✓ Branch 5 taken 106 times.
|
1832 | if (LinearSolverTraits::isNonOverlapping(parallelHelper_->gridView())) |
| 401 | { | ||
| 402 | using PTraits = typename LinearSolverTraits::template ParallelNonoverlapping<Matrix, Vector>; | ||
| 403 | 286 | solveParallel_<PTraits>(A, x, b); | |
| 404 | } | ||
| 405 | else | ||
| 406 | { | ||
| 407 | using PTraits = typename LinearSolverTraits::template ParallelOverlapping<Matrix, Vector>; | ||
| 408 | 438 | solveParallel_<PTraits>(A, x, b); | |
| 409 | } | ||
| 410 | } | ||
| 411 | else | ||
| 412 | 287 | solveSequential_(A, x, b); | |
| 413 | } | ||
| 414 | else | ||
| 415 | { | ||
| 416 | 51 | solveSequential_(A, x, b); | |
| 417 | } | ||
| 418 | 1011 | } | |
| 419 | |||
| 420 | template<class ParallelTraits> | ||
| 421 | 1448 | void solveParallel_(Matrix& A, Vector& x, Vector& b) | |
| 422 | { | ||
| 423 | using LinearOperator = typename ParallelTraits::LinearOperator; | ||
| 424 | |||
| 425 |
2/2✓ Branch 0 taken 30 times.
✓ Branch 1 taken 694 times.
|
1448 | if (firstCall_) |
| 426 | 60 | initSolverFactories<Matrix, LinearOperator>(); | |
| 427 | |||
| 428 | 2896 | prepareLinearAlgebraParallel<LinearSolverTraits, ParallelTraits>(A, b, *parallelHelper_); | |
| 429 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2896 | auto linearOperator = std::make_shared<LinearOperator>(A, *comm_); |
| 430 | |||
| 431 | // construct solver | ||
| 432 |
3/8✓ Branch 1 taken 724 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 724 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 724 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4344 | auto solver = getSolverFromFactory_(linearOperator); |
| 433 | |||
| 434 | // solve linear system | ||
| 435 |
2/4✓ Branch 1 taken 724 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 724 times.
✗ Branch 5 not taken.
|
2896 | solver->apply(x, b, result_); |
| 436 | 1448 | } | |
| 437 | #endif // HAVE_MPI | ||
| 438 | |||
| 439 | 338 | void solveSequential_(Matrix& A, Vector& x, Vector& b) | |
| 440 | { | ||
| 441 | // construct linear operator | ||
| 442 | using Traits = typename LinearSolverTraits::template Sequential<Matrix, Vector>; | ||
| 443 | using LinearOperator = typename Traits::LinearOperator; | ||
| 444 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
338 | auto linearOperator = std::make_shared<LinearOperator>(A); |
| 445 | |||
| 446 |
2/2✓ Branch 0 taken 15 times.
✓ Branch 1 taken 323 times.
|
338 | if (firstCall_) |
| 447 |
1/2✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
|
15 | initSolverFactories<Matrix, LinearOperator>(); |
| 448 | |||
| 449 | // construct solver | ||
| 450 |
3/8✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 338 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 338 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
1014 | auto solver = getSolverFromFactory_(linearOperator); |
| 451 | |||
| 452 | // solve linear system | ||
| 453 |
2/4✓ Branch 1 taken 338 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 338 times.
✗ Branch 5 not taken.
|
676 | solver->apply(x, b, result_); |
| 454 | 338 | } | |
| 455 | |||
| 456 | template<class LinearOperator> | ||
| 457 | 1062 | auto getSolverFromFactory_(std::shared_ptr<LinearOperator>& fop) | |
| 458 | { | ||
| 459 |
3/12✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 51 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 51 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
|
3186 | try { return Dune::getSolverFromFactory(fop, params_); } |
| 460 | ✗ | catch(Dune::Exception& e) | |
| 461 | { | ||
| 462 | ✗ | std::cerr << "Could not create solver with factory" << std::endl; | |
| 463 | ✗ | std::cerr << e.what() << std::endl; | |
| 464 | ✗ | throw e; | |
| 465 | } | ||
| 466 | } | ||
| 467 | |||
| 468 | const std::string paramGroup_; | ||
| 469 | #if HAVE_MPI | ||
| 470 | std::unique_ptr<ParallelISTLHelper<LinearSolverTraits>> parallelHelper_; | ||
| 471 | std::shared_ptr<Comm> comm_; | ||
| 472 | #endif | ||
| 473 | bool isParallel_ = false; | ||
| 474 | bool firstCall_; | ||
| 475 | |||
| 476 | std::shared_ptr<ScalarProduct> scalarProduct_; | ||
| 477 | Dune::SolverCategory::Category solverCategory_; | ||
| 478 | Dune::InverseOperatorResult result_; | ||
| 479 | Dune::ParameterTree params_; | ||
| 480 | std::string name_; | ||
| 481 | }; | ||
| 482 | |||
| 483 | template<int i> | ||
| 484 | struct IstlSolverFactoryBackendImplHelper | ||
| 485 | { | ||
| 486 | template<class ...T> | ||
| 487 | using type = IstlSolverFactoryBackend<T...>; | ||
| 488 | }; | ||
| 489 | |||
| 490 | template<> | ||
| 491 | struct IstlSolverFactoryBackendImplHelper<1> | ||
| 492 | { | ||
| 493 | template<class T> | ||
| 494 | using type = OldIstlSolverFactoryBackend<T>; | ||
| 495 | }; | ||
| 496 | |||
| 497 | } // end namespace Detail | ||
| 498 | |||
| 499 | /*! | ||
| 500 | * \ingroup Linear | ||
| 501 | * \brief A linear solver using the dune-istl solver factory | ||
| 502 | * to choose the solver and preconditioner at runtime. | ||
| 503 | * \note the solvers are configured via the input file | ||
| 504 | * \note The template arguments are LinearSolverTraits, LinearAlgebraTraits | ||
| 505 | */ | ||
| 506 | template<class ...T> | ||
| 507 | using IstlSolverFactoryBackend = typename Detail::IstlSolverFactoryBackendImplHelper<sizeof...(T)>::template type<T...>; | ||
| 508 | |||
| 509 | } // end namespace Dumux | ||
| 510 | |||
| 511 | #endif | ||
| 512 |