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