GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/linear/pdesolver.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 72 74 97.3%
Functions: 53 55 96.4%
Branches: 52 133 39.1%
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 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
1/2
✓ Branch 2 taken 40 times.
✗ Branch 3 not taken.
 40 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 42 const Communication& comm = Dune::MPIHelper::getCommunication(),
84 const std::string& paramGroup = "")
85 : ParentType(assembler, linearSolver)
86 84 , paramGroup_(paramGroup)
87 42 , reuseMatrix_(false)
88
3/6
✓ 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.
 126 , 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
1/2
✓ Branch 1 taken 42 times.
✗ Branch 2 not taken.
 42 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 2983 << std::flush;
117
118 ///////////////
119 // assemble
120 ///////////////
121
122 // linearize the problem at the current solution
123 2983 assembleTimer.start();
124
2/2
✓ Branch 0 taken 950 times.
✓ Branch 1 taken 2033 times.
 2983 if (reuseMatrix_)
125 950 this->assembler().assembleResidual(vars);
126 else
127 2033 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 2983 << clearRemainingLine << std::flush;
142
143 // solve the resulting linear equation system
144 2983 solveTimer.start();
145
146 // set the delta vector to zero
147 2983 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
2/4
✓ Branch 1 taken 2928 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2928 times.
✗ Branch 5 not taken.
 2983 << clearRemainingLine << std::flush;
163
164 // update the current solution and secondary variables
165 2983 updateTimer.start();
166
1/2
✓ Branch 1 taken 2983 times.
✗ Branch 2 not taken.
 2983 auto uCurrent = Backend::dofs(vars);
167
2/3
✓ Branch 1 taken 183 times.
✓ Branch 2 taken 2800 times.
✗ Branch 3 not taken.
 2983 Backend::axpy(-1.0, deltaU, uCurrent);
168 2983 Backend::update(vars, uCurrent);
169 if constexpr (!assemblerExportsVariables)
170
1/2
✓ Branch 1 taken 2950 times.
✗ Branch 2 not taken.
 2983 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
4/8
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2928 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 2928 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 2928 times.
✗ Branch 13 not taken.
 2928 << assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
180
4/8
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2928 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 2928 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 2928 times.
✗ Branch 13 not taken.
 2928 << solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
181
4/8
✓ Branch 2 taken 2928 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2928 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 2928 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 2928 times.
✗ Branch 13 not taken.
 2928 << updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
182
1/2
✓ Branch 1 taken 2928 times.
✗ Branch 2 not taken.
 2928 << "\n";
183 }
184
185
1/2
✓ Branch 0 taken 2967 times.
✗ Branch 1 not taken.
 2983 return true;
186 5966 }
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 1 void setVerbosity(int val)
218
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 { verbosity_ = val; }
219
220 /*!
221 * \brief Returns true if the solver ought to be chatty.
222 */
223 2 int verbosity() const
224
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 { return verbosity_ ; }
225
226 /*!
227 * \brief Returns the parameter group
228 */
229 const std::string& paramGroup() const
230 { return paramGroup_; }
231
232 /*!
233 * \brief Set whether the matrix should be reused
234 * \note If this is set to true, the matrix will not be assembled. Make
235 * sure there is an assembled matrix that can be reused before
236 * setting this flag to true.
237 */
238 11 void reuseMatrix(bool reuse = true)
239 {
240 11 reuseMatrix_ = reuse;
241
242 if constexpr (Detail::LinearPDESolver::linearSolverHasSetMatrix<LinearSolver, JacobianMatrix>())
243 if (reuseMatrix_)
244
1/2
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
 11 this->linearSolver().setMatrix(this->assembler().jacobian());
245 }
246
247 private:
248
249 2983 virtual bool solveLinearSystem_(ResidualVector& deltaU)
250 {
251 if constexpr (Detail::LinearPDESolver::linearSolverHasSetMatrix<LinearSolver, JacobianMatrix>())
252 {
253
2/2
✓ Branch 0 taken 950 times.
✓ Branch 1 taken 33 times.
 983 if (reuseMatrix_)
254 950 return this->linearSolver().solve(deltaU, this->assembler().residual());
255 }
256 else
257 {
258
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2000 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2000 if (reuseMatrix_ && comm_.size() > 1)
259 ✗ DUNE_THROW(Dune::NotImplemented,
260 "Reuse matrix for parallel runs with a solver that doesn't support the setMatrix interface"
261 );
262 }
263
264 2033 assert(this->checkSizesOfSubMatrices(this->assembler().jacobian()) && "Matrix blocks have wrong sizes!");
265
266 2033 return this->linearSolver().solve(
267 2033 this->assembler().jacobian(),
268 deltaU,
269 2033 this->assembler().residual()
270 2000 );
271 }
272
273 //! initialize the parameters by reading from the parameter tree
274 42 void initParams_(const std::string& group = "")
275 {
276
2/2
✓ Branch 0 taken 41 times.
✓ Branch 1 taken 1 times.
 42 verbosity_ = comm_.rank() == 0 ? getParamFromGroup<int>(group, "LinearPDESolver.Verbosity", 2) : 0;
277 42 enableDynamicOutput_ = getParamFromGroup<bool>(group, "LinearPDESolver.EnableDynamicOutput", true);
278 42 }
279
280 //! sets verbosity-level
281 int verbosity_;
282
283 //! further parameters
284 bool enableDynamicOutput_;
285
286 //! the parameter group for getting parameters from the parameter tree
287 std::string paramGroup_;
288
289 //! check if the matrix is supposed to be reused
290 bool reuseMatrix_;
291
292 //! The communication object (for distributed memory parallelism)
293 Communication comm_;
294 };
295
296 } // end namespace Dumux
297
298 #endif
299