GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/linear/pdesolver.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 64 69 92.8%
Functions: 53 58 91.4%
Branches: 55 136 40.4%
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 A high-level solver interface for a linear PDE solver
11 */
12 #ifndef DUMUX_LINEAR_PDE_SOLVER_HH
13 #define DUMUX_LINEAR_PDE_SOLVER_HH
14
15 #include <cmath>
16 #include <memory>
17 #include <iostream>
18 #include <type_traits>
19
20 #include <dune/common/timer.hh>
21 #include <dune/common/exceptions.hh>
22 #include <dune/common/parallel/mpicommunication.hh>
23 #include <dune/common/parallel/mpihelper.hh>
24 #include <dune/common/std/type_traits.hh>
25 #include <dune/istl/bvector.hh>
26 #include <dune/istl/multitypeblockvector.hh>
27
28 #include <dumux/common/parameters.hh>
29 #include <dumux/common/exceptions.hh>
30 #include <dumux/common/typetraits/vector.hh>
31 #include <dumux/common/timeloop.hh>
32 #include <dumux/common/pdesolver.hh>
33 #include <dumux/common/variablesbackend.hh>
34
35 #include <dumux/linear/matrixconverter.hh>
36
37 namespace Dumux::Detail::LinearPDESolver {
38
39 template <class Solver, class Matrix>
40 using SetMatrixDetector = decltype(std::declval<Solver>().setMatrix(std::declval<std::shared_ptr<Matrix>>()));
41
42 template<class Solver, class Matrix>
43 static constexpr bool linearSolverHasSetMatrix()
44 { return Dune::Std::is_detected<SetMatrixDetector, Solver, Matrix>::value; }
45
46 } // end namespace Dumux::Detail::LinearPDESolver
47
48 namespace Dumux {
49
50 /*!
51 * \ingroup Linear
52 * \brief An implementation of a linear PDE solver
53 * \tparam Assembler the assembler
54 * \tparam LinearSolver the linear solver
55 * \tparam Comm the communication object used to communicate with all processes
56 * \note If you want to specialize only some methods but are happy with the
57 * defaults of the reference solver, derive your solver from
58 * this class and simply overload the required methods.
59 */
60 template <class Assembler, class LinearSolver,
61 class Comm = Dune::Communication<Dune::MPIHelper::MPICommunicator>>
62 class LinearPDESolver : public PDESolver<Assembler, LinearSolver>
63 {
64 using ParentType = PDESolver<Assembler, LinearSolver>;
65 using Scalar = typename Assembler::Scalar;
66 using JacobianMatrix = typename Assembler::JacobianMatrix;
67 using SolutionVector = typename Assembler::SolutionVector;
68 using ResidualVector = typename Assembler::ResidualType;
69 using TimeLoop = TimeLoopBase<Scalar>;
70 using Backend = VariablesBackend<typename ParentType::Variables>;
71 using LinearAlgebraNativeBackend = VariablesBackend<ResidualVector>;
72 static constexpr bool assemblerExportsVariables = Detail::PDESolver::assemblerExportsVariables<Assembler>;
73
74 public:
75 using typename ParentType::Variables;
76 using Communication = Comm;
77
78 /*!
79 * \brief The Constructor
80 */
81 42 LinearPDESolver(std::shared_ptr<Assembler> assembler,
82 std::shared_ptr<LinearSolver> linearSolver,
83
2/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 42 const Communication& comm = Dune::MPIHelper::getCommunication(),
84 const std::string& paramGroup = "")
85 : ParentType(assembler, linearSolver)
86 , paramGroup_(paramGroup)
87 , reuseMatrix_(false)
88
3/8
✓ Branch 3 taken 42 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 42 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 42 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 168 , comm_(comm)
89 {
90
1/2
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
 42 initParams_(paramGroup);
91
92 // set the linear system (matrix & residual) in the assembler
93
2/4
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 42 times.
✗ Branch 5 not taken.
 84 this->assembler().setLinearSystem();
94 42 }
95
96 /*!
97 * \brief The Constructor
98 */
99 LinearPDESolver(std::shared_ptr<Assembler> assembler,
100 std::shared_ptr<LinearSolver> linearSolver,
101 const std::string& paramGroup)
102 : LinearPDESolver(assembler, linearSolver, Dune::MPIHelper::getCommunication(), paramGroup)
103 {}
104
105 /*!
106 * \brief Solve a linear PDE system
107 */
108 2983 bool apply(Variables& vars) override
109 {
110 2983 Dune::Timer assembleTimer(false);
111 2983 Dune::Timer solveTimer(false);
112 2983 Dune::Timer updateTimer(false);
113
114
3/4
✓ Branch 0 taken 2928 times.
✓ Branch 1 taken 55 times.
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
 2983 if (verbosity_ >= 1 && enableDynamicOutput_)
115 2928 std::cout << "Assemble: r(x^k) = dS/dt + div F - q; M = grad r"
116 2928 << std::flush;
117
118 ///////////////
119 // assemble
120 ///////////////
121
122 // linearize the problem at the current solution
123
1/2
✓ Branch 0 taken 2983 times.
✗ Branch 1 not taken.
 2983 assembleTimer.start();
124
2/2
✓ Branch 0 taken 950 times.
✓ Branch 1 taken 2033 times.
 2983 if (reuseMatrix_)
125 1900 this->assembler().assembleResidual(vars);
126 else
127 4066 this->assembler().assembleJacobianAndResidual(vars);
128 2983 assembleTimer.stop();
129
130 ///////////////
131 // linear solve
132 ///////////////
133
134 // Clear the current line using an ansi escape
135 // sequence. for an explanation see
136 // http://en.wikipedia.org/wiki/ANSI_escape_code
137 2983 const char clearRemainingLine[] = { 0x1b, '[', 'K', 0 };
138
139
3/4
✓ Branch 0 taken 2928 times.
✓ Branch 1 taken 55 times.
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
 2983 if (verbosity_ >= 1 && enableDynamicOutput_)
140 std::cout << "\rSolve: M deltax^k = r"
141 8784 << clearRemainingLine << std::flush;
142
143 // solve the resulting linear equation system
144
1/2
✓ Branch 0 taken 2983 times.
✗ Branch 1 not taken.
 2983 solveTimer.start();
145
146 // set the delta vector to zero
147
1/4
✓ Branch 2 taken 2983 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 11932 ResidualVector deltaU = LinearAlgebraNativeBackend::zeros(Backend::size(Backend::dofs(vars)));
148
149 // solve by calling the appropriate implementation depending on whether the linear solver
150 // is capable of handling MultiType matrices or not
151
1/2
✓ Branch 1 taken 2983 times.
✗ Branch 2 not taken.
 2983 bool converged = solveLinearSystem_(deltaU);
152 2983 solveTimer.stop();
153
154
1/2
✓ Branch 0 taken 2983 times.
✗ Branch 1 not taken.
 2983 if (!converged)
155 return false;
156
157 ///////////////
158 // update
159 ///////////////
160
3/4
✓ Branch 0 taken 2928 times.
✓ Branch 1 taken 55 times.
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
 2983 if (verbosity_ >= 1 && enableDynamicOutput_)
161 std::cout << "\rUpdate: x^(k+1) = x^k - deltax^k"
162
1/2
✓ Branch 3 taken 2928 times.
✗ Branch 4 not taken.
 8784 << clearRemainingLine << std::flush;
163
164 // update the current solution and secondary variables
165
1/2
✓ Branch 0 taken 2983 times.
✗ Branch 1 not taken.
 2983 updateTimer.start();
166
1/4
✓ Branch 1 taken 2983 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2983 auto uCurrent = Backend::dofs(vars);
167 2983 Backend::axpy(-1.0, deltaU, uCurrent);
168
1/2
✓ Branch 1 taken 2983 times.
✗ Branch 2 not taken.
 2983 Backend::update(vars, uCurrent);
169 if constexpr (!assemblerExportsVariables)
170
2/4
✓ Branch 1 taken 2950 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2950 times.
✗ Branch 5 not taken.
 5966 this->assembler().updateGridVariables(Backend::dofs(vars));
171 2983 updateTimer.stop();
172
173
2/2
✓ Branch 0 taken 2928 times.
✓ Branch 1 taken 55 times.
 2983 if (verbosity_ >= 1)
174 {
175 2928 const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
176
1/2
✓ Branch 0 taken 2928 times.
✗ Branch 1 not taken.
 2928 if (enableDynamicOutput_)
177
1/2
✓ Branch 1 taken 2928 times.
✗ Branch 2 not taken.
 2928 std::cout << '\r';
178 2928 std::cout << "Assemble/solve/update time: "
179
2/4
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 2928 times.
✗ Branch 8 not taken.
 5856 << assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
180
2/4
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 2928 times.
✗ Branch 8 not taken.
 5856 << solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
181
2/4
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
✓ Branch 7 taken 2928 times.
✗ Branch 8 not taken.
 5856 << updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
182
1/2
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
 5856 << "\n";
183 }
184
185
1/2
✓ Branch 0 taken 2983 times.
✗ Branch 1 not taken.
 2983 return true;
186 }
187
188 /*!
189 * \brief Solve a linear PDE system
190 */
191 983 void solve(Variables& vars) override
192 {
193 983 bool converged = apply(vars);
194
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 983 times.
 983 if (!converged)
195 ✗ DUNE_THROW(NumericalProblem, "Linear solver didn't converge.\n");
196 983 }
197
198 /*!
199 * \brief output statistics / report
200 * \todo Implement some solver statistics output
201 */
202 void report(std::ostream& sout = std::cout) const
203 {}
204
205 /*!
206 * \brief Suggest a new time-step size based on the old time-step size.
207 * \note For compatibility with other PDE solvers (e.g. Newton)
208 */
209 Scalar suggestTimeStepSize(Scalar oldTimeStep) const
210 {
211 return oldTimeStep;
212 }
213
214 /*!
215 * \brief Specifies if the solver ought to be chatty.
216 */
217 [[deprecated("Will be removed after release 3.9. Use setVerbosity(int).")]]
218 void setVerbose(bool val)
219 { verbosity_ = val ? 1 : 0; }
220
221 /*!
222 * \brief Returns true if the solver ought to be chatty.
223 */
224 [[deprecated("Will be removed after release 3.9. Use int verbsity().")]]
225 bool verbose() const
226 { return verbosity_ > 0; }
227
228 /*!
229 * \brief Specifies if the solver ought to be chatty.
230 */
231 ✗ void setVerbosity(int val)
232
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 { verbosity_ = val; }
233
234 /*!
235 * \brief Returns true if the solver ought to be chatty.
236 */
237 ✗ int verbosity() const
238 ✗ { return verbosity_ ; }
239
240 /*!
241 * \brief Returns the parameter group
242 */
243 const std::string& paramGroup() const
244 { return paramGroup_; }
245
246 /*!
247 * \brief Set whether the matrix should be reused
248 * \note If this is set to true, the matrix will not be assembled. Make
249 * sure there is an assembled matrix that can be reused before
250 * setting this flag to true.
251 */
252 void reuseMatrix(bool reuse = true)
253 {
254 11 reuseMatrix_ = reuse;
255
256 if constexpr (Detail::LinearPDESolver::linearSolverHasSetMatrix<LinearSolver, JacobianMatrix>())
257 if (reuseMatrix_)
258
4/8
✓ 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.
✓ Branch 10 taken 11 times.
✗ Branch 11 not taken.
 44 this->linearSolver().setMatrix(this->assembler().jacobian());
259 }
260
261 private:
262
263 2983 virtual bool solveLinearSystem_(ResidualVector& deltaU)
264 {
265 if constexpr (Detail::LinearPDESolver::linearSolverHasSetMatrix<LinearSolver, JacobianMatrix>())
266 {
267
2/2
✓ Branch 0 taken 950 times.
✓ Branch 1 taken 33 times.
 983 if (reuseMatrix_)
268 3800 return this->linearSolver().solve(deltaU, this->assembler().residual());
269 }
270 else
271 {
272
1/6
✗ Branch 0 not taken.
✓ Branch 1 taken 2000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2000 if (reuseMatrix_ && comm_.size() > 1)
273 ✗ DUNE_THROW(Dune::NotImplemented,
274 "Reuse matrix for parallel runs with a solver that doesn't support the setMatrix interface"
275 );
276 }
277
278 4066 assert(this->checkSizesOfSubMatrices(this->assembler().jacobian()) && "Matrix blocks have wrong sizes!");
279
280 2033 return this->linearSolver().solve(
281 4066 this->assembler().jacobian(),
282 deltaU,
283 2033 this->assembler().residual()
284 8132 );
285 }
286
287 //! initialize the parameters by reading from the parameter tree
288 42 void initParams_(const std::string& group = "")
289 {
290
4/4
✓ Branch 0 taken 41 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 41 times.
✓ Branch 3 taken 1 times.
 84 verbosity_ = comm_.rank() == 0 ? getParamFromGroup<int>(group, "LinearPDESolver.Verbosity", 2) : 0;
291 42 enableDynamicOutput_ = getParamFromGroup<bool>(group, "LinearPDESolver.EnableDynamicOutput", true);
292 42 }
293
294 //! sets verbosity-level
295 int verbosity_;
296
297 //! further parameters
298 bool enableDynamicOutput_;
299
300 //! the parameter group for getting parameters from the parameter tree
301 std::string paramGroup_;
302
303 //! check if the matrix is supposed to be reused
304 bool reuseMatrix_;
305
306 //! The communication object (for distributed memory parallelism)
307 Communication comm_;
308 };
309
310 } // end namespace Dumux
311
312 #endif
313