| 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 Preconditioned iterative solver for the incompressible Stokes problem | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_LINEAR_STOKES_SOLVER_HH | ||
| 13 | #define DUMUX_LINEAR_STOKES_SOLVER_HH | ||
| 14 | |||
| 15 | #include <type_traits> | ||
| 16 | #include <memory> | ||
| 17 | #include <tuple> | ||
| 18 | |||
| 19 | #include <dune/common/parametertree.hh> | ||
| 20 | #include <dune/common/hybridutilities.hh> | ||
| 21 | #include <dune/common/exceptions.hh> | ||
| 22 | |||
| 23 | #include <dune/istl/matrixindexset.hh> | ||
| 24 | #include <dune/istl/preconditioner.hh> | ||
| 25 | #include <dune/istl/preconditioners.hh> | ||
| 26 | #include <dune/istl/solvers.hh> | ||
| 27 | #include <dune/istl/paamg/amg.hh> | ||
| 28 | |||
| 29 | #include <dumux/common/math.hh> | ||
| 30 | #include <dumux/common/exceptions.hh> | ||
| 31 | #include <dumux/linear/solver.hh> | ||
| 32 | #include <dumux/linear/preconditioners.hh> | ||
| 33 | #include <dumux/linear/linearsolverparameters.hh> | ||
| 34 | #include <dumux/linear/parallelhelpers.hh> | ||
| 35 | #include <dumux/linear/linearsolvertraits.hh> | ||
| 36 | #include <dumux/linear/parallelmatrixadapter.hh> | ||
| 37 | #include <dumux/discretization/extrusion.hh> | ||
| 38 | #include <dumux/linear/symmetrize_constraints.hh> | ||
| 39 | #include <dumux/assembly/jacobianpattern.hh> | ||
| 40 | |||
| 41 | namespace Dumux::Detail { | ||
| 42 | |||
| 43 | /*! | ||
| 44 | * \ingroup Linear | ||
| 45 | * \brief A Stokes preconditioner (saddle-point problem) for the problem | ||
| 46 | * \f$ | ||
| 47 | \begin{pmatrix} A & B \\ C & 0 \end{pmatrix} | ||
| 48 | \begin{pmatrix} u \\ p \end{pmatrix} = | ||
| 49 | \begin{pmatrix} f \\ g \end{pmatrix}, | ||
| 50 | * \f$ | ||
| 51 | * | ||
| 52 | * where A is the discrete velocity operator and B and C are discrete gradient and divergence operators. | ||
| 53 | * This preconditioner is especially suited for solving the incompressible Stokes equations. | ||
| 54 | * | ||
| 55 | * \tparam M Type of the matrix. | ||
| 56 | * \tparam X Type of the update. | ||
| 57 | * \tparam Y Type of the defect. | ||
| 58 | * \tparam l Preconditioner block level | ||
| 59 | */ | ||
| 60 | template<class M, class X, class Y, int l = 2> | ||
| 61 | class StokesPreconditioner : public Dune::Preconditioner<X,Y> | ||
| 62 | { | ||
| 63 | static_assert(Dumux::isMultiTypeBlockMatrix<M>::value && M::M() == 2 && M::N() == 2, "Expects a 2x2 MultiTypeBlockMatrix."); | ||
| 64 | static_assert(l == 2, "StokesPreconditioner expects a block level of 2 for Matrix type M."); | ||
| 65 | |||
| 66 | using A = std::decay_t<decltype(std::declval<M>()[Dune::Indices::_0][Dune::Indices::_0])>; | ||
| 67 | using U = std::decay_t<decltype(std::declval<X>()[Dune::Indices::_0])>; | ||
| 68 | |||
| 69 | using P = std::decay_t<decltype(std::declval<M>()[Dune::Indices::_1][Dune::Indices::_1])>; | ||
| 70 | using V = std::decay_t<decltype(std::declval<X>()[Dune::Indices::_1])>; | ||
| 71 | |||
| 72 | enum class Mode { symmetric, triangular, diagonal }; | ||
| 73 | |||
| 74 | public: | ||
| 75 | //! \brief The matrix type the preconditioner is for. | ||
| 76 | using matrix_type = M; | ||
| 77 | //! \brief The domain type of the preconditioner. | ||
| 78 | using domain_type = X; | ||
| 79 | //! \brief The range type of the preconditioner. | ||
| 80 | using range_type = Y; | ||
| 81 | //! \brief The field type of the preconditioner. | ||
| 82 | using field_type = typename X::field_type; | ||
| 83 | //! \brief Scalar type underlying the field_type. | ||
| 84 | using scalar_field_type = Dune::Simd::Scalar<field_type>; | ||
| 85 | //! \brief the type of the pressure operator | ||
| 86 | using PressureLinearOperator = Dune::MatrixAdapter<P,V,V>; | ||
| 87 | |||
| 88 | /*! | ||
| 89 | * \brief Constructor | ||
| 90 | * \param fullLinearOperator the Stokes linear operator | ||
| 91 | * \param pressureLinearOperator the linear operator to be used for the pressure space preconditioner | ||
| 92 | * \param params a parameter tree for the preconditioner configuration | ||
| 93 | */ | ||
| 94 | 12 | StokesPreconditioner( | |
| 95 | const std::shared_ptr<const Dune::AssembledLinearOperator<M,X,Y>>& fullLinearOperator, | ||
| 96 | const std::shared_ptr<const PressureLinearOperator>& pressureLinearOperator, | ||
| 97 | const Dune::ParameterTree& params | ||
| 98 | ) | ||
| 99 | 12 | : matrix_(fullLinearOperator->getmat()) | |
| 100 | 12 | , pmatrix_(pressureLinearOperator->getmat()) | |
| 101 |
1/2✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
|
12 | , verbosity_(params.get<int>("verbosity")) |
| 102 |
2/4✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 12 times.
✗ Branch 6 not taken.
|
24 | , paramGroup_(params.get<std::string>("ParameterGroup")) |
| 103 | { | ||
| 104 | 12 | const auto mode = getParamFromGroup<std::string>( | |
| 105 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | paramGroup_, "LinearSolver.Preconditioner.Mode", "Diagonal" |
| 106 | ); | ||
| 107 | |||
| 108 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
|
12 | if (mode == "Symmetric") |
| 109 | ✗ | mode_ = Mode::symmetric; | |
| 110 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 10 times.
|
12 | else if (mode == "Triangular") |
| 111 | 2 | mode_ = Mode::triangular; | |
| 112 | else | ||
| 113 | 10 | mode_ = Mode::diagonal; | |
| 114 | |||
| 115 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | initPreconditioner_(params); |
| 116 | 12 | } | |
| 117 | |||
| 118 | /*! | ||
| 119 | * \brief Prepare the preconditioner. | ||
| 120 | */ | ||
| 121 | 12 | void pre(X& update, Y& currentDefect) override {} | |
| 122 | |||
| 123 | /*! | ||
| 124 | * \brief Apply the preconditioner | ||
| 125 | * | ||
| 126 | * \param update The update to be computed. | ||
| 127 | * \param currentDefect The current defect. | ||
| 128 | * | ||
| 129 | * The currentDefect has be be in a consistent representation, | ||
| 130 | * Definition 2.3 Blatt and Bastian (2009) https://doi.org/10.1504/IJCSE.2008.021112 | ||
| 131 | * The update is initially zero. At exit the update has to be | ||
| 132 | * in a consistent representation. This usually requires communication. | ||
| 133 | */ | ||
| 134 | 1711 | void apply(X& update, const Y& currentDefect) override | |
| 135 | { | ||
| 136 | using namespace Dune::Indices; | ||
| 137 | |||
| 138 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1711 times.
|
1711 | if (mode_ == Mode::symmetric) |
| 139 | { | ||
| 140 | ✗ | update = applyRightBlock_(currentDefect); | |
| 141 | ✗ | update = applyDiagBlock_(update); | |
| 142 | ✗ | update = applyLeftBlock_(update); | |
| 143 | } | ||
| 144 |
2/2✓ Branch 0 taken 229 times.
✓ Branch 1 taken 1482 times.
|
1711 | else if (mode_ == Mode::triangular) |
| 145 | { | ||
| 146 | 458 | update = applyRightBlock_(currentDefect); | |
| 147 | 458 | update = applyDiagBlock_(update); | |
| 148 | } | ||
| 149 | else | ||
| 150 | 2964 | update = applyDiagBlock_(currentDefect); | |
| 151 | 1711 | } | |
| 152 | |||
| 153 | /*! | ||
| 154 | * \brief Clean up. | ||
| 155 | */ | ||
| 156 | 12 | void post(X& update) override {} | |
| 157 | |||
| 158 | //! Category of the preconditioner (see SolverCategory::Category) | ||
| 159 | 12 | Dune::SolverCategory::Category category() const override | |
| 160 | { | ||
| 161 | 12 | return Dune::SolverCategory::sequential; | |
| 162 | } | ||
| 163 | |||
| 164 | private: | ||
| 165 | 229 | X applyRightBlock_(const Y& d) | |
| 166 | { | ||
| 167 | using namespace Dune::Indices; | ||
| 168 | |||
| 169 | // right bit of LDU decomposition | ||
| 170 | // applied to d | ||
| 171 | // | ||
| 172 | // | I 0 | | ||
| 173 | // | -C*inv(A) I | | ||
| 174 | // | ||
| 175 | |||
| 176 | 229 | auto dTmp0 = d[_0]; | |
| 177 |
1/5✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 2 taken 229 times.
|
229 | auto vTmp = d; vTmp = 0.0; |
| 178 | |||
| 179 | // invert velocity block (apply to first block of d) | ||
| 180 |
1/2✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
|
229 | applyPreconditionerForA_(vTmp[_0], dTmp0); |
| 181 | |||
| 182 | // then multiply with C | ||
| 183 |
1/2✓ Branch 2 taken 229 times.
✗ Branch 3 not taken.
|
229 | matrix_[_1][_0].mv(vTmp[_0], vTmp[_1]); |
| 184 | |||
| 185 | // and subtract from d | ||
| 186 |
1/2✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
|
229 | auto v = d; |
| 187 |
1/2✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
|
229 | v[_0] = vTmp[_0]; // already do A^-1 d of the diagonal block because we already computed it here |
| 188 | 229 | v[_1] -= vTmp[_1]; | |
| 189 | 229 | return v; | |
| 190 | 458 | } | |
| 191 | |||
| 192 | 1711 | X applyDiagBlock_(const Y& d) | |
| 193 | { | ||
| 194 | using namespace Dune::Indices; | ||
| 195 | |||
| 196 | // diagonal middle bit of LDU decomposition | ||
| 197 |
1/2✓ Branch 1 taken 1711 times.
✗ Branch 2 not taken.
|
1711 | auto dTmp = d; |
| 198 |
2/2✓ Branch 0 taken 6717600 times.
✓ Branch 1 taken 818 times.
|
6719311 | auto v = d; v = 0.0; |
| 199 | |||
| 200 | // invert velocity block | ||
| 201 |
2/2✓ Branch 0 taken 1482 times.
✓ Branch 1 taken 229 times.
|
1711 | if (mode_ == Mode::diagonal) |
| 202 |
1/2✓ Branch 1 taken 1482 times.
✗ Branch 2 not taken.
|
1482 | applyPreconditionerForA_(v[_0], dTmp[_0]); |
| 203 | |||
| 204 | // reuse the already computed A^-1 d (see applyRightBlock_) | ||
| 205 | else | ||
| 206 |
1/2✓ Branch 1 taken 229 times.
✗ Branch 2 not taken.
|
229 | v[_0] = dTmp[_0]; |
| 207 | |||
| 208 | // invert pressure block | ||
| 209 |
1/2✓ Branch 1 taken 1711 times.
✗ Branch 2 not taken.
|
1711 | applyPreconditionerForP_(v[_1], dTmp[_1]); |
| 210 | |||
| 211 | 1711 | return v; | |
| 212 | 1711 | } | |
| 213 | |||
| 214 | ✗ | X applyLeftBlock_(const Y& d) | |
| 215 | { | ||
| 216 | using namespace Dune::Indices; | ||
| 217 | |||
| 218 | // left bit of LDU decomposition | ||
| 219 | // applied to d | ||
| 220 | // | ||
| 221 | // | I -inv(A)*B | | ||
| 222 | // | 0 I | | ||
| 223 | // | ||
| 224 | |||
| 225 | ✗ | auto dTmp = d; | |
| 226 | ✗ | auto vTmp = d; vTmp = 0.0; | |
| 227 | |||
| 228 | // multiply with B | ||
| 229 | ✗ | matrix_[_0][_1].umv(dTmp[_1], vTmp[_0]); | |
| 230 | |||
| 231 | // invert velocity block (apply to first block of d) | ||
| 232 | ✗ | auto vTmp0 = vTmp[_0]; vTmp0 = 0.0; | |
| 233 | ✗ | applyPreconditionerForA_(vTmp0, vTmp[_0]); | |
| 234 | |||
| 235 | // and subtract from d | ||
| 236 | ✗ | auto v = d; | |
| 237 | ✗ | v[_0] -= vTmp0; | |
| 238 | |||
| 239 | ✗ | return v; | |
| 240 | ✗ | } | |
| 241 | |||
| 242 | 12 | void initPreconditioner_(const Dune::ParameterTree& params) | |
| 243 | { | ||
| 244 | using namespace Dune::Indices; | ||
| 245 | |||
| 246 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
|
12 | if (getParamFromGroup<bool>(paramGroup_, "LinearSolver.DirectSolverForVelocity", false)) |
| 247 | { | ||
| 248 | #if HAVE_UMFPACK | ||
| 249 | ✗ | directSolverForA_ = std::make_shared<Dune::UMFPack<A>>(matrix_[_0][_0], verbosity_); | |
| 250 | using Wrap = Dune::InverseOperator2Preconditioner<Dune::InverseOperator<U, U>>; | ||
| 251 | ✗ | preconditionerForA_ = std::make_shared<Wrap>(*directSolverForA_); | |
| 252 | #else | ||
| 253 | DUNE_THROW(Dune::InvalidStateException, "Selected direct solver but UMFPack is not available."); | ||
| 254 | #endif | ||
| 255 | } | ||
| 256 | else | ||
| 257 | { | ||
| 258 | using VelLinearOperator = Dune::MatrixAdapter<A, U, U>; | ||
| 259 |
1/2✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
|
12 | auto lopV = std::make_shared<VelLinearOperator>(matrix_[_0][_0]); |
| 260 |
2/6✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 12 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
12 | preconditionerForA_ = std::make_shared< |
| 261 | Dune::Amg::AMG<VelLinearOperator, U, Dumux::ParMTSSOR<A,U,U>> | ||
| 262 | >(lopV, params); | ||
| 263 | 12 | } | |
| 264 | |||
| 265 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
|
12 | if (getParamFromGroup<bool>(paramGroup_, "LinearSolver.DirectSolverForPressure", false)) |
| 266 | { | ||
| 267 | #if HAVE_UMFPACK | ||
| 268 | ✗ | directSolverForP_ = std::make_shared<Dune::UMFPack<P>>(pmatrix_, verbosity_); | |
| 269 | using Wrap = Dune::InverseOperator2Preconditioner<Dune::InverseOperator<V, V>>; | ||
| 270 | ✗ | preconditionerForP_ = std::make_shared<Wrap>(*directSolverForP_); | |
| 271 | #else | ||
| 272 | DUNE_THROW(Dune::InvalidStateException, "Selected direct solver but UMFPack is not available."); | ||
| 273 | #endif | ||
| 274 | } | ||
| 275 | else | ||
| 276 | { | ||
| 277 | using PressJacobi = Dumux::ParMTJac<P, V, V>; | ||
| 278 |
2/2✓ Branch 1 taken 1 times.
✓ Branch 2 taken 11 times.
|
12 | std::size_t numIterations = pmatrix_.nonzeroes() == pmatrix_.N() ? 1 : 10; |
| 279 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 12 times.
|
12 | preconditionerForP_ = std::make_shared<PressJacobi>(pmatrix_, numIterations, 1.0); |
| 280 | } | ||
| 281 | 12 | } | |
| 282 | |||
| 283 | template<class Sol, class Rhs> | ||
| 284 | 1711 | void applyPreconditionerForA_(Sol& sol, Rhs& rhs) const | |
| 285 | { | ||
| 286 | 1711 | preconditionerForA_->pre(sol, rhs); | |
| 287 | 1711 | preconditionerForA_->apply(sol, rhs); | |
| 288 | 1711 | preconditionerForA_->post(sol); | |
| 289 | 1711 | } | |
| 290 | |||
| 291 | template<class Sol, class Rhs> | ||
| 292 | 1711 | void applyPreconditionerForP_(Sol& sol, Rhs& rhs) const | |
| 293 | { | ||
| 294 | 1711 | preconditionerForP_->pre(sol, rhs); | |
| 295 | 1711 | preconditionerForP_->apply(sol, rhs); | |
| 296 | 1711 | preconditionerForP_->post(sol); | |
| 297 | 1711 | } | |
| 298 | |||
| 299 | //! \brief The matrix we operate on. | ||
| 300 | const M& matrix_; | ||
| 301 | //! \brief The matrix we operate on. | ||
| 302 | const P& pmatrix_; | ||
| 303 | //! \brief The verbosity level | ||
| 304 | const int verbosity_; | ||
| 305 | |||
| 306 | std::shared_ptr<Dune::Preconditioner<U, U>> preconditionerForA_; | ||
| 307 | std::shared_ptr<Dune::Preconditioner<V, V>> preconditionerForP_; | ||
| 308 | std::shared_ptr<Dune::InverseOperator<U, U>> directSolverForA_; | ||
| 309 | std::shared_ptr<Dune::InverseOperator<V, V>> directSolverForP_; | ||
| 310 | |||
| 311 | const std::string paramGroup_; | ||
| 312 | Mode mode_; | ||
| 313 | }; | ||
| 314 | |||
| 315 | } // end namespace Detail | ||
| 316 | |||
| 317 | namespace Dumux { | ||
| 318 | |||
| 319 | /*! | ||
| 320 | * \ingroup Linear | ||
| 321 | * \brief Preconditioned iterative solver for the incompressible Stokes problem | ||
| 322 | * \note Uses StokesPreconditioner as preconditioner (tailored to the incompressible Stokes problem) | ||
| 323 | * \note No MPI parallelization implemented, some shared-memory parallelism is enabled | ||
| 324 | */ | ||
| 325 | template<class Matrix, class Vector, class VelocityGG, class PressureGG> | ||
| 326 | class StokesSolver | ||
| 327 | : public LinearSolver | ||
| 328 | { | ||
| 329 | using Preconditioner = Detail::StokesPreconditioner<Matrix, Vector, Vector>; | ||
| 330 | public: | ||
| 331 | /*! | ||
| 332 | * \brief Constructor | ||
| 333 | * \param vGridGeometry grid geometry of the velocity discretization | ||
| 334 | * \param pGridGeometry grid geometry of the pressure discretization | ||
| 335 | * \param dirichletDofs a vector (same size and shape as right hand side) where the dirichlet dofs are marked with 1.0 | ||
| 336 | * and all other entries are 0.0. | ||
| 337 | * \param paramGroup group prefix when looking up keys in the parameter tree | ||
| 338 | */ | ||
| 339 | 11 | StokesSolver(std::shared_ptr<const VelocityGG> vGridGeometry, | |
| 340 | std::shared_ptr<const PressureGG> pGridGeometry, | ||
| 341 | const Vector& dirichletDofs, | ||
| 342 | const std::string& paramGroup = "") | ||
| 343 | : LinearSolver(paramGroup) | ||
| 344 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | , vGridGeometry_(std::move(vGridGeometry)) |
| 345 | 11 | , pGridGeometry_(std::move(pGridGeometry)) | |
| 346 |
2/4✓ Branch 2 taken 11 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 11 times.
✗ Branch 6 not taken.
|
11 | , dirichletDofs_(dirichletDofs) |
| 347 | { | ||
| 348 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | params_ = LinearSolverParameters<LinearSolverTraits<VelocityGG>>::createParameterTree(this->paramGroup()); |
| 349 |
2/4✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
|
11 | density_ = getParamFromGroup<double>(this->paramGroup(), "Component.LiquidDensity"); |
| 350 | |||
| 351 |
3/6✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 11 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 11 times.
✗ Branch 8 not taken.
|
22 | if (hasParamInGroup(this->paramGroup(), "Component.LiquidDynamicViscosity")) |
| 352 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | viscosity_ = getParamFromGroup<double>(this->paramGroup(), "Component.LiquidDynamicViscosity"); |
| 353 | ✗ | else if (hasParamInGroup(this->paramGroup(), "Component.LiquidKinematicViscosity")) | |
| 354 | ✗ | viscosity_ = getParamFromGroup<double>(this->paramGroup(), "Component.LiquidKinematicViscosity") * density_; | |
| 355 | else | ||
| 356 | { | ||
| 357 | ✗ | const std::string group = this->paramGroup() == "" ? "Component" : this->paramGroup() + ".Component"; | |
| 358 | ✗ | DUNE_THROW(ParameterException, "Stokes solver requires the parameters" | |
| 359 | << " LiquidDynamicViscosity or LiquidKinematicViscosity" | ||
| 360 | << " in parameter group [" << group << "]"); | ||
| 361 | ✗ | } | |
| 362 | |||
| 363 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | weight_ = getParamFromGroup<double>(this->paramGroup(), "LinearSolver.Preconditioner.MassMatrixWeight", 1.0); |
| 364 |
2/4✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 11 times.
✗ Branch 7 not taken.
|
11 | solverType_ = getParamFromGroup<std::string>(this->paramGroup(), "LinearSolver.Type", "gmres"); |
| 365 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | scalarProduct_ = std::make_shared<Dune::ScalarProduct<Vector>>(); |
| 366 | 11 | } | |
| 367 | |||
| 368 | 12 | bool solve(const Matrix& A, Vector& x, const Vector& b) | |
| 369 | { | ||
| 370 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | auto bTmp = b; |
| 371 | 12 | auto ATmp = A; | |
| 372 | |||
| 373 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | return applyIterativeSolver_(ATmp, x, bTmp); |
| 374 | 12 | } | |
| 375 | |||
| 376 | 8 | Scalar norm(const Vector& b) const | |
| 377 | { | ||
| 378 |
1/2✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
|
8 | return scalarProduct_->norm(b); |
| 379 | } | ||
| 380 | |||
| 381 | std::string name() const | ||
| 382 | { | ||
| 383 | return "Block-preconditioned Stokes solver"; | ||
| 384 | } | ||
| 385 | |||
| 386 | const Dune::InverseOperatorResult& result() const | ||
| 387 | { | ||
| 388 | return result_; | ||
| 389 | } | ||
| 390 | |||
| 391 | private: | ||
| 392 | 12 | bool applyIterativeSolver_(Matrix& A, Vector& x, Vector& b) | |
| 393 | { | ||
| 394 | // make Dirichlet boundary conditions symmetric | ||
| 395 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | if (getParamFromGroup<bool>(this->paramGroup(), "LinearSolver.SymmetrizeDirichlet", true)) |
| 396 | 12 | symmetrizeConstraints(A, b, dirichletDofs_); | |
| 397 | |||
| 398 | // make Matrix symmetric on the block-scale | ||
| 399 | using namespace Dune::Indices; | ||
| 400 | 12 | A[_1] *= -1.0/density_; | |
| 401 | 24 | b[_1] *= -1.0/density_; | |
| 402 | |||
| 403 | 12 | auto op = std::make_shared<Dumux::ParallelMultiTypeMatrixAdapter<Matrix, Vector, Vector>>(A); | |
| 404 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | auto pop = makePressureLinearOperator_<typename Preconditioner::PressureLinearOperator>(); |
| 405 |
3/8✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 12 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
|
12 | auto preconditioner = std::make_shared<Preconditioner>(op, pop, params_.sub("preconditioner")); |
| 406 |
14/32✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 6 times.
✓ Branch 3 taken 6 times.
✓ Branch 5 taken 6 times.
✓ Branch 6 taken 6 times.
✓ Branch 8 taken 6 times.
✓ Branch 9 taken 6 times.
✓ Branch 11 taken 6 times.
✓ Branch 12 taken 6 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 6 times.
✓ Branch 17 taken 6 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 6 times.
✗ Branch 21 not taken.
✓ Branch 25 taken 6 times.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✓ Branch 28 taken 6 times.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✗ Branch 10 not taken.
✗ Branch 13 not taken.
✗ Branch 16 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
|
30 | params_["verbose"] = pGridGeometry_->gridView().comm().rank() == 0 ? params_["verbose"] : "0"; |
| 407 | |||
| 408 | // defaults to restarted GMRes | ||
| 409 | 12 | std::unique_ptr<Dune::InverseOperator<Vector, Vector>> solver; | |
| 410 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
|
12 | if (solverType_ == "minres") |
| 411 | ✗ | solver = std::make_unique<Dune::MINRESSolver<Vector>>(op, scalarProduct_, preconditioner, params_); | |
| 412 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
|
12 | else if (solverType_ == "bicgstab") |
| 413 | ✗ | solver = std::make_unique<Dune::BiCGSTABSolver<Vector>>(op, scalarProduct_, preconditioner, params_); | |
| 414 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | else if (solverType_ == "gmres") |
| 415 |
1/4✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
12 | solver = std::make_unique<Dune::RestartedGMResSolver<Vector>>(op, scalarProduct_, preconditioner, params_); |
| 416 | else | ||
| 417 |
1/26✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✓ Branch 40 taken 12 times.
✗ Branch 41 not taken.
|
12 | DUNE_THROW(Dune::NotImplemented, "Solver choice " << solverType_ << " is not implemented"); |
| 418 | |||
| 419 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | solver->apply(x, b, result_); |
| 420 | |||
| 421 | 12 | return result_.converged; | |
| 422 | 48 | } | |
| 423 | |||
| 424 | template<class LinearOperator> | ||
| 425 | 12 | std::shared_ptr<LinearOperator> makePressureLinearOperator_() | |
| 426 | { | ||
| 427 | using M = typename LinearOperator::matrix_type; | ||
| 428 | 12 | auto massMatrix = createMassMatrix_<M>(); | |
| 429 |
1/2✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
|
12 | return std::make_shared<LinearOperator>(massMatrix); |
| 430 | 12 | } | |
| 431 | |||
| 432 | template<class M, class DM = typename PressureGG::DiscretizationMethod> | ||
| 433 | 12 | std::shared_ptr<M> createMassMatrix_() | |
| 434 | { | ||
| 435 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | auto massMatrix = std::make_shared<M>(); |
| 436 |
1/2✓ Branch 1 taken 12 times.
✗ Branch 2 not taken.
|
12 | massMatrix->setBuildMode(M::random); |
| 437 |
3/5✓ Branch 1 taken 6 times.
✓ Branch 2 taken 6 times.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✗ Branch 3 not taken.
|
12 | const auto numDofs = pGridGeometry_->numDofs(); |
| 438 | |||
| 439 | if constexpr (DM{} == DiscretizationMethods::cctpfa || DM{} == DiscretizationMethods::ccmpfa) | ||
| 440 | { | ||
| 441 | 11 | Dune::MatrixIndexSet pattern; | |
| 442 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | pattern.resize(numDofs, numDofs); |
| 443 |
2/2✓ Branch 0 taken 54606 times.
✓ Branch 1 taken 11 times.
|
54617 | for (unsigned int globalI = 0; globalI < numDofs; ++globalI) |
| 444 | 54606 | pattern.add(globalI, globalI); | |
| 445 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | pattern.exportIdx(*massMatrix); |
| 446 | |||
| 447 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | const auto& gv = pGridGeometry_->gridView(); |
| 448 |
1/2✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
|
11 | auto fvGeometry = localView(*pGridGeometry_); |
| 449 |
4/6✓ Branch 1 taken 54610 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34612 times.
✓ Branch 4 taken 1 times.
✓ Branch 7 taken 5658 times.
✗ Branch 8 not taken.
|
134891 | for (const auto& element : elements(gv)) |
| 450 | { | ||
| 451 |
1/2✓ Branch 1 taken 34606 times.
✗ Branch 2 not taken.
|
54606 | fvGeometry.bindElement(element); |
| 452 |
4/5✓ Branch 1 taken 54606 times.
✓ Branch 2 taken 28948 times.
✓ Branch 3 taken 54606 times.
✓ Branch 4 taken 25658 times.
✗ Branch 0 not taken.
|
109212 | for (const auto& scv : scvs(fvGeometry)) |
| 453 | { | ||
| 454 | using Extrusion = Extrusion_t<PressureGG>; | ||
| 455 |
1/3✓ Branch 1 taken 54606 times.
✗ Branch 2 not taken.
✗ Branch 0 not taken.
|
54606 | const auto dofIndex = scv.dofIndex(); |
| 456 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 25658 times.
|
54606 | if (element.partitionType() == Dune::GhostEntity) // do not modify ghosts |
| 457 | ✗ | (*massMatrix)[dofIndex][dofIndex] = 1.0; | |
| 458 | else | ||
| 459 |
1/2✓ Branch 1 taken 54606 times.
✗ Branch 2 not taken.
|
54606 | (*massMatrix)[dofIndex][dofIndex] += weight_*Extrusion::volume(fvGeometry, scv)/(2.0*viscosity_); |
| 460 | } | ||
| 461 | } | ||
| 462 | |||
| 463 | 11 | return massMatrix; | |
| 464 | 11 | } | |
| 465 | else if constexpr (DM{} == DiscretizationMethods::box || DM{} == DiscretizationMethods::fcdiamond) | ||
| 466 | { | ||
| 467 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
|
1 | getJacobianPattern<true>(*pGridGeometry_).exportIdx(*massMatrix); |
| 468 |
1/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1 | (*massMatrix) = 0.0; |
| 469 | |||
| 470 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto& gv = pGridGeometry_->gridView(); |
| 471 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | auto fvGeometry = localView(*pGridGeometry_); |
| 472 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | std::vector<Dune::FieldVector<double, 1>> values; |
| 473 | |||
| 474 |
4/6✓ Branch 1 taken 28948 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28948 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 28948 times.
✓ Branch 7 taken 1 times.
|
57897 | for (const auto& element : elements(gv)) |
| 475 | { | ||
| 476 | 28948 | const auto geometry = element.geometry(); | |
| 477 |
3/8✓ Branch 1 taken 28948 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 28948 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 28948 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
28948 | const auto& localBasis = pGridGeometry_->feCache().get(element.type()).localBasis(); |
| 478 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
28948 | const auto numLocalDofs = localBasis.size(); |
| 479 |
1/2✓ Branch 1 taken 28948 times.
✗ Branch 2 not taken.
|
28948 | values.resize(numLocalDofs); |
| 480 | |||
| 481 |
1/2✓ Branch 1 taken 28948 times.
✗ Branch 2 not taken.
|
28948 | fvGeometry.bindElement(element); |
| 482 |
3/4✓ Branch 1 taken 86844 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 86844 times.
✓ Branch 4 taken 28948 times.
|
115792 | for (const auto& scvJ : scvs(fvGeometry)) |
| 483 | { | ||
| 484 | // Use mid-point rule (only works for linear ansatz functions) | ||
| 485 |
1/2✓ Branch 1 taken 86844 times.
✗ Branch 2 not taken.
|
86844 | const auto globalJ = scvJ.dofIndex(); |
| 486 |
1/2✓ Branch 1 taken 86844 times.
✗ Branch 2 not taken.
|
86844 | const auto qWeightJ = PressureGG::Extrusion::volume(fvGeometry, scvJ); |
| 487 |
2/4✓ Branch 1 taken 86844 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 86844 times.
✗ Branch 5 not taken.
|
86844 | const auto quadPos = geometry.local(scvJ.center()); |
| 488 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
86844 | localBasis.evaluateFunction(quadPos, values); |
| 489 | |||
| 490 |
2/2✓ Branch 0 taken 260532 times.
✓ Branch 1 taken 86844 times.
|
347376 | for (const auto& scvI : scvs(fvGeometry)) |
| 491 | { | ||
| 492 |
1/2✓ Branch 1 taken 260532 times.
✗ Branch 2 not taken.
|
260532 | const auto valueIJ = values[scvI.localDofIndex()]*qWeightJ/(2.0*viscosity_); |
| 493 |
1/2✓ Branch 1 taken 260532 times.
✗ Branch 2 not taken.
|
260532 | (*massMatrix)[scvI.dofIndex()][globalJ][0][0] += valueIJ; |
| 494 | } | ||
| 495 | } | ||
| 496 | } | ||
| 497 | |||
| 498 | 1 | return massMatrix; | |
| 499 | 1 | } | |
| 500 | |||
| 501 | DUNE_THROW(Dune::NotImplemented, "Mass matrix for discretization method not implemented"); | ||
| 502 | ✗ | } | |
| 503 | |||
| 504 | double density_, viscosity_, weight_; | ||
| 505 | Dune::InverseOperatorResult result_; | ||
| 506 | Dune::ParameterTree params_; | ||
| 507 | std::shared_ptr<const VelocityGG> vGridGeometry_; | ||
| 508 | std::shared_ptr<const PressureGG> pGridGeometry_; | ||
| 509 | const Vector& dirichletDofs_; | ||
| 510 | std::string solverType_; | ||
| 511 | |||
| 512 | std::shared_ptr<Dune::ScalarProduct<Vector>> scalarProduct_; | ||
| 513 | }; | ||
| 514 | |||
| 515 | } // end namespace Dumux | ||
| 516 | |||
| 517 | #endif | ||
| 518 |