GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/assembly/fvlocalresidual.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 38 43 88.4%
Functions: 1262 1922 65.7%
Branches: 70 95 73.7%
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 Assembly
10 * \brief The element-wise residual for finite volume schemes
11 */
12 #ifndef DUMUX_FV_LOCAL_RESIDUAL_HH
13 #define DUMUX_FV_LOCAL_RESIDUAL_HH
14
15 #include <dune/common/exceptions.hh>
16 #include <dune/istl/bvector.hh>
17
18 #include <dumux/common/properties.hh>
19 #include <dumux/common/timeloop.hh>
20 #include <dumux/common/reservedblockvector.hh>
21 #include <dumux/common/numeqvector.hh>
22 #include <dumux/discretization/method.hh>
23 #include <dumux/discretization/extrusion.hh>
24
25 namespace Dumux {
26
27 /*!
28 * \ingroup Assembly
29 * \brief The element-wise residual for finite volume schemes
30 * \note This class defines the interface used by the assembler using
31 * static polymorphism. Implementations are specialized for a certain discretization scheme
32 */
33 template<class TypeTag>
34 class FVLocalResidual
35 {
36 using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
37 using Problem = GetPropType<TypeTag, Properties::Problem>;
38 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
39 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
40 using Element = typename GridView::template Codim<0>::Entity;
41 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
42 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
43 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
44 using SubControlVolume = typename GridGeometry::SubControlVolume;
45 using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
46 using Extrusion = Extrusion_t<GridGeometry>;
47 using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
48 using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
49 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
50 using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
51 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
52 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
53 using TimeLoop = TimeLoopBase<Scalar>;
54
55 public:
56 //! the container storing all element residuals
57 using ElementResidualVector = ReservedBlockVector<NumEqVector, FVElementGeometry::maxNumElementScvs>;
58
59 //! the constructor
60 FVLocalResidual(const Problem* problem,
61 const TimeLoop* timeLoop = nullptr)
62 : problem_(problem)
63 , timeLoop_(timeLoop)
64 {}
65
66 /*!
67 * \name User interface
68 * \note The following methods are usually expensive to evaluate
69 * They are useful for outputting / postprocessing residual information.
70 */
71 // \{
72
73 /*!
74 * \brief Compute the storage term for the current solution.
75 *
76 * This can be used to figure out how much of each conservation
77 * quantity is inside the element.
78 *
79 * \param problem The problem to solve
80 * \param element The DUNE Codim<0> entity for which the storage
81 * term ought to be calculated
82 * \param gridGeometry The finite-volume grid geometry
83 * \param gridVariables The grid variables (volume and flux variables)
84 * \param sol The solution vector
85 */
86 ElementResidualVector evalStorage(const Problem& problem,
87 const Element &element,
88 const GridGeometry& gridGeometry,
89 const GridVariables& gridVariables,
90 const SolutionVector& sol) const
91 {
92 // make sure FVElementGeometry and volume variables are bound to the element
93 const auto fvGeometry = localView(gridGeometry).bind(element);
94 const auto elemVolVars = localView(gridVariables.curGridVolVars()).bind(element, fvGeometry, sol);
95
96 ElementResidualVector storage(fvGeometry.numScv());
97
98 // calculate the amount of conservation each quantity inside
99 // all sub control volumes
100 for (auto&& scv : scvs(fvGeometry))
101 {
102 auto localScvIdx = scv.localDofIndex();
103 const auto& volVars = elemVolVars[scv];
104 storage[localScvIdx] = asImp().computeStorage(problem, scv, volVars);
105 storage[localScvIdx] *= Extrusion::volume(fvGeometry, scv) * volVars.extrusionFactor();
106 }
107
108 return storage;
109 }
110
111 // \}
112
113 /*!
114 * \name Main interface
115 * \note Methods used by the assembler to compute derivatives and residual
116 */
117 // \{
118
119 /*!
120 * \brief Compute the storage local residual, i.e. the deviation of the
121 * storage term from zero for instationary problems.
122 *
123 * \param element The DUNE Codim<0> entity for which the residual
124 * ought to be calculated
125 * \param fvGeometry The finite-volume geometry of the element
126 * \param prevElemVolVars The volume averaged variables for all
127 * sub-control volumes of the element at the previous time level
128 * \param curElemVolVars The volume averaged variables for all
129 * sub-control volumes of the element at the current time level
130 */
131 246802391 ElementResidualVector evalStorage(const Element& element,
132 const FVElementGeometry& fvGeometry,
133 const ElementVolumeVariables& prevElemVolVars,
134 const ElementVolumeVariables& curElemVolVars) const
135 {
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 198365708 times.
 246802391 assert(timeLoop_ && "no time loop set for storage term evaluation");
137
138 // initialize the residual vector for all scvs in this element
139 491120174 ElementResidualVector residual(fvGeometry.numScv());
140
141 // evaluate the volume terms (storage + source terms)
142 // forward to the local residual specialized for the discretization methods
143
5/5
✓ Branch 0 taken 455709621 times.
✓ Branch 1 taken 202010120 times.
✓ Branch 2 taken 423661326 times.
✓ Branch 3 taken 170048045 times.
✓ Branch 4 taken 1214804 times.
 1028136996 for (auto&& scv : scvs(fvGeometry))
144 536961822 asImp().evalStorage(residual, this->problem(), element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
145
146 246802391 return residual;
147 }
148
149 307670544 ElementResidualVector evalFluxAndSource(const Element& element,
150 const FVElementGeometry& fvGeometry,
151 const ElementVolumeVariables& elemVolVars,
152 const ElementFluxVariablesCache& elemFluxVarsCache,
153 const ElementBoundaryTypes &bcTypes) const
154 {
155 // initialize the residual vector for all scvs in this element
156 610122180 ElementResidualVector residual(fvGeometry.numScv());
157
158 // evaluate the volume terms (storage + source terms)
159 // forward to the local residual specialized for the discretization methods
160
5/5
✓ Branch 0 taken 534625889 times.
✓ Branch 1 taken 243651684 times.
✓ Branch 2 taken 498420064 times.
✓ Branch 3 taken 209822936 times.
✓ Branch 4 taken 2581782 times.
 1287541182 for (auto&& scv : scvs(fvGeometry))
161 677356932 asImp().evalSource(residual, this->problem(), element, fvGeometry, elemVolVars, scv);
162
163 // forward to the local residual specialized for the discretization methods
164
5/5
✓ Branch 0 taken 1006820250 times.
✓ Branch 1 taken 264527150 times.
✓ Branch 2 taken 1004160037 times.
✓ Branch 3 taken 259630363 times.
✓ Branch 4 taken 2581782 times.
 2663955359 for (auto&& scvf : scvfs(fvGeometry))
165 1655786195 asImp().evalFlux(residual, this->problem(), element, fvGeometry, elemVolVars, bcTypes, elemFluxVarsCache, scvf);
166
167 307670542 return residual;
168 }
169
170 // \}
171
172
173 /*!
174 * \name Model specific interface
175 * \note The following method are the model specific implementations of the local residual
176 */
177 // \{
178
179 /*!
180 * \brief Calculate the source term of the equation
181 *
182 * \param problem The problem to solve
183 * \param scv The sub-control volume over which we integrate the storage term
184 * \param volVars The volume variables associated with the scv
185 * \note has to be implemented by the model specific residual class
186 *
187 */
188 NumEqVector computeStorage(const Problem& problem,
189 const SubControlVolume& scv,
190 const VolumeVariables& volVars) const
191 {
192 DUNE_THROW(Dune::NotImplemented, "This model does not implement a storage method!");
193 }
194
195 /*!
196 * \brief Calculate the source term of the equation
197 *
198 * \param problem The problem to solve
199 * \param element The DUNE Codim<0> entity for which the residual
200 * ought to be calculated
201 * \param fvGeometry The finite-volume geometry of the element
202 * \param elemVolVars The volume variables associated with the element stencil
203 * \param scv The sub-control volume over which we integrate the source term
204 * \note This is the default implementation for all models as sources are computed
205 * in the user interface of the problem
206 *
207 */
208 507047182 NumEqVector computeSource(const Problem& problem,
209 const Element& element,
210 const FVElementGeometry& fvGeometry,
211 const ElementVolumeVariables& elemVolVars,
212 const SubControlVolume &scv) const
213 {
214
4/4
✓ Branch 0 taken 548464 times.
✓ Branch 1 taken 246360 times.
✓ Branch 2 taken 91723 times.
✓ Branch 3 taken 148 times.
 573422903 NumEqVector source(0.0);
215
216 // add contributions from volume flux sources
217
16/20
✓ Branch 0 taken 556142 times.
✓ Branch 1 taken 25367505 times.
✓ Branch 2 taken 11681582 times.
✓ Branch 3 taken 25384905 times.
✓ Branch 4 taken 13435425 times.
✓ Branch 5 taken 30014328 times.
✓ Branch 6 taken 13435425 times.
✓ Branch 7 taken 30014328 times.
✓ Branch 8 taken 1512624 times.
✓ Branch 9 taken 148 times.
✓ Branch 10 taken 1512624 times.
✓ Branch 11 taken 148 times.
✓ Branch 12 taken 4126024 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 4126024 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 3323 times.
✗ Branch 17 not taken.
✓ Branch 18 taken 3323 times.
✗ Branch 19 not taken.
 911633535 source += problem.source(element, fvGeometry, elemVolVars, scv);
218
219 // add contribution from possible point sources
220
24/30
✓ Branch 0 taken 9756987 times.
✓ Branch 1 taken 466828466 times.
✓ Branch 2 taken 9756987 times.
✓ Branch 3 taken 466828466 times.
✓ Branch 4 taken 9756987 times.
✓ Branch 5 taken 466828466 times.
✓ Branch 6 taken 13435425 times.
✓ Branch 7 taken 30014328 times.
✓ Branch 8 taken 13435425 times.
✓ Branch 9 taken 30014328 times.
✓ Branch 10 taken 13435425 times.
✓ Branch 11 taken 30014328 times.
✓ Branch 12 taken 1512624 times.
✓ Branch 13 taken 148 times.
✓ Branch 14 taken 1512624 times.
✓ Branch 15 taken 148 times.
✓ Branch 16 taken 1512624 times.
✓ Branch 17 taken 148 times.
✓ Branch 18 taken 4126024 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 4126024 times.
✗ Branch 21 not taken.
✓ Branch 22 taken 4126024 times.
✗ Branch 23 not taken.
✓ Branch 24 taken 3323 times.
✗ Branch 25 not taken.
✓ Branch 26 taken 3323 times.
✗ Branch 27 not taken.
✓ Branch 28 taken 3323 times.
✗ Branch 29 not taken.
 1745396817 if (!problem.pointSourceMap().empty())
221
0/7
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 39472464 source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv);
222
223 507047180 return source;
224 }
225
226 /*!
227 * \brief Calculate the flux term of the equation
228 *
229 * \param problem The problem to solve
230 * \param element The DUNE Codim<0> entity for which the residual
231 * ought to be calculated
232 * \param fvGeometry The finite-volume geometry of the element
233 * \param elemVolVars The volume variables associated with the element stencil
234 * \param scvf The sub-control volume over which we integrate the flux
235 * \param elemFluxVarsCache the flux variable caches for the element's flux stencils
236 *
237 * \note has to be implemented by the model specific residual class
238 *
239 */
240 NumEqVector computeFlux(const Problem& problem,
241 const Element& element,
242 const FVElementGeometry& fvGeometry,
243 const ElementVolumeVariables& elemVolVars,
244 const SubControlVolumeFace& scvf,
245 const ElementFluxVariablesCache& elemFluxVarsCache) const
246 {
247 DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!");
248 }
249
250 // \}
251
252 /*!
253 * \name Discretization specific interface
254 * \note The following method are the discretization specific wrapper methods
255 */
256 // \{
257
258 /*!
259 * \brief Compute the storage local residual, i.e. the deviation of the
260 * storage term from zero for instationary problems.
261 *
262 * \param residual The residual vector to fill
263 * \param problem The problem to solve
264 * \param element The DUNE Codim<0> entity for which the residual
265 * ought to be calculated
266 * \param fvGeometry The finite-volume geometry of the element
267 * \param prevElemVolVars The volume averaged variables for all
268 * sub-control volumes of the element at the previous time level
269 * \param curElemVolVars The volume averaged variables for all
270 * sub-control volumes of the element at the current time level
271 * \param scv The sub control volume the storage term is integrated over
272 */
273 161603656 void evalStorage(ElementResidualVector& residual,
274 const Problem& problem,
275 const Element& element,
276 const FVElementGeometry& fvGeometry,
277 const ElementVolumeVariables& prevElemVolVars,
278 const ElementVolumeVariables& curElemVolVars,
279 const SubControlVolume& scv) const
280 {
281
1/2
✓ Branch 0 taken 14460000 times.
✗ Branch 1 not taken.
 161603656 const auto& curVolVars = curElemVolVars[scv];
282
1/2
✓ Branch 0 taken 14460000 times.
✗ Branch 1 not taken.
 161603656 const auto& prevVolVars = prevElemVolVars[scv];
283
284 // mass balance within the element. this is the
285 // \f$\frac{m}{\partial t}\f$ term if using implicit or explicit
286 // euler as time discretization.
287 //
288 // TODO: We might need a more explicit way for
289 // doing the time discretization...
290
291 //! Compute storage with the model specific storage residual
292
1/2
✓ Branch 0 taken 21870800 times.
✗ Branch 1 not taken.
 161603656 NumEqVector prevStorage = asImp().computeStorage(problem, scv, prevVolVars);
293
1/2
✓ Branch 0 taken 21870800 times.
✗ Branch 1 not taken.
 161603656 NumEqVector storage = asImp().computeStorage(problem, scv, curVolVars);
294
295 161603656 prevStorage *= prevVolVars.extrusionFactor();
296 140274361 storage *= curVolVars.extrusionFactor();
297
298 161603656 storage -= prevStorage;
299 323207312 storage *= Extrusion::volume(fvGeometry, scv);
300 161603656 storage /= timeLoop_->timeStepSize();
301
302 363853420 residual[scv.localDofIndex()] += storage;
303 161603656 }
304
305 /*!
306 * \brief Compute the source local residual, i.e. the deviation of the
307 * source term from zero.
308 *
309 * \param residual The residual vector to fill
310 * \param problem The problem to solve
311 * \param element The DUNE Codim<0> entity for which the residual
312 * ought to be calculated
313 * \param fvGeometry The finite-volume geometry of the element
314 * \param curElemVolVars The volume averaged variables for all
315 * sub-control volumes of the element at the current time level
316 * \param scv The sub control volume the source term is integrated over
317 */
318 743761860 void evalSource(ElementResidualVector& residual,
319 const Problem& problem,
320 const Element& element,
321 const FVElementGeometry& fvGeometry,
322 const ElementVolumeVariables& curElemVolVars,
323 const SubControlVolume& scv) const
324 {
325 //! Compute source with the model specific storage residual
326
2/2
✓ Branch 0 taken 13616464 times.
✓ Branch 1 taken 28671254 times.
 743761860 const auto& curVolVars = curElemVolVars[scv];
327
2/3
✓ Branch 0 taken 18993351 times.
✓ Branch 1 taken 55135473 times.
✗ Branch 2 not taken.
 743761860 NumEqVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
328 1757552610 source *= Extrusion::volume(fvGeometry, scv)*curVolVars.extrusionFactor();
329
330 //! subtract source from local rate (sign convention in user interface)
331 1685485464 residual[scv.localDofIndex()] -= source;
332 743761858 }
333
334 /*!
335 * \brief Compute the flux local residual, i.e. the deviation of the
336 * flux term from zero.
337 * \param residual The residual vector to fill
338 * \param problem The problem to solve
339 * \param element The DUNE Codim<0> entity for which the residual
340 * ought to be calculated
341 * \param fvGeometry The finite-volume geometry of the element
342 * \param elemVolVars The volume averaged variables for all
343 * sub-control volumes of the element at the current time level
344 * \param elemBcTypes the boundary types for the boundary entities of an elements
345 * \param elemFluxVarsCache The flux variable caches for the element stencil
346 * \param scvf The sub control volume face the flux term is integrated over
347 */
348 void evalFlux(ElementResidualVector& residual,
349 const Problem& problem,
350 const Element& element,
351 const FVElementGeometry& fvGeometry,
352 const ElementVolumeVariables& elemVolVars,
353 const ElementBoundaryTypes& elemBcTypes,
354 const ElementFluxVariablesCache& elemFluxVarsCache,
355 const SubControlVolumeFace& scvf) const {}
356
357 /*!
358 * \brief Compute the flux local residual, i.e. the deviation of the
359 * flux term from zero.
360 *
361 * \param problem The problem to solve
362 * \param element The DUNE Codim<0> entity for which the residual
363 * ought to be calculated
364 * \param fvGeometry The finite-volume geometry of the element
365 * \param elemVolVars The volume averaged variables for all
366 * sub-control volumes of the element at the current time level
367 * \param elemFluxVarsCache The flux variable caches for the element stencil
368 * \param scvf The sub control volume face the flux term is integrated over
369 */
370 NumEqVector evalFlux(const Problem& problem,
371 const Element& element,
372 const FVElementGeometry& fvGeometry,
373 const ElementVolumeVariables& elemVolVars,
374 const ElementFluxVariablesCache& elemFluxVarsCache,
375 const SubControlVolumeFace& scvf) const
376 {
377 return asImp().evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
378 }
379
380 //\}
381
382 /*!
383 * \name Interfaces for analytic Jacobian computation
384 */
385 // \{
386
387 //! Compute the derivative of the storage residual
388 template<class PartialDerivativeMatrix>
389 void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
390 const Problem& problem,
391 const Element& element,
392 const FVElementGeometry& fvGeometry,
393 const VolumeVariables& curVolVars,
394 const SubControlVolume& scv) const
395 {
396 DUNE_THROW(Dune::NotImplemented, "analytic storage derivative");
397 }
398
399 //! Compute the derivative of the source residual
400 template<class PartialDerivativeMatrix>
401 void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
402 const Problem& problem,
403 const Element& element,
404 const FVElementGeometry& fvGeometry,
405 const VolumeVariables& curVolVars,
406 const SubControlVolume& scv) const
407 {
408 DUNE_THROW(Dune::NotImplemented, "analytic source derivative");
409 }
410
411 //! Compute the derivative of the flux residual
412 template<class PartialDerivativeMatrices, class T = TypeTag>
413 std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod != DiscretizationMethods::box, void>
414 addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
415 const Problem& problem,
416 const Element& element,
417 const FVElementGeometry& fvGeometry,
418 const ElementVolumeVariables& curElemVolVars,
419 const ElementFluxVariablesCache& elemFluxVarsCache,
420 const SubControlVolumeFace& scvf) const
421 {
422 DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for cell-centered models");
423 }
424
425 //! Compute the derivative of the flux residual for the box method
426 template<class JacobianMatrix, class T = TypeTag>
427 std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::box, void>
428 addFluxDerivatives(JacobianMatrix& A,
429 const Problem& problem,
430 const Element& element,
431 const FVElementGeometry& fvGeometry,
432 const ElementVolumeVariables& curElemVolVars,
433 const ElementFluxVariablesCache& elemFluxVarsCache,
434 const SubControlVolumeFace& scvf) const
435 {
436 DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for box models");
437 }
438
439 //! Compute the derivative of the Dirichlet flux residual for cell-centered schemes
440 template<class PartialDerivativeMatrices>
441 void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
442 const Problem& problem,
443 const Element& element,
444 const FVElementGeometry& fvGeometry,
445 const ElementVolumeVariables& curElemVolVars,
446 const ElementFluxVariablesCache& elemFluxVarsCache,
447 const SubControlVolumeFace& scvf) const
448 {
449 DUNE_THROW(Dune::NotImplemented, "analytic Dirichlet flux derivative");
450 }
451
452 //! Compute the derivative of Robin type boundary conditions ("solution dependent Neumann")
453 template<class PartialDerivativeMatrices>
454 void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
455 const Problem& problem,
456 const Element& element,
457 const FVElementGeometry& fvGeometry,
458 const ElementVolumeVariables& curElemVolVars,
459 const ElementFluxVariablesCache& elemFluxVarsCache,
460 const SubControlVolumeFace& scvf) const
461 {
462 DUNE_THROW(Dune::NotImplemented, "analytic Robin flux derivative");
463 }
464
465 //\}
466
467 /*!
468 * \name Interfaces accessed by local residual implementations
469 */
470 // \{
471
472 //! the problem
473 ✗ const Problem& problem() const
474 ✗ { return *problem_; }
475
476 //! the timeloop for instationary problems
477 //! calling this for stationary leads to undefined behaviour
478 ✗ const TimeLoop& timeLoop() const
479 ✗ { return *timeLoop_; }
480
481 //! returns true if the residual is stationary
482 ✗ bool isStationary() const
483
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 9 { return !timeLoop_; }
484
485 // \}
486 protected:
487
488 Implementation &asImp()
489 { return *static_cast<Implementation*>(this); }
490
491 const Implementation &asImp() const
492 { return *static_cast<const Implementation*>(this); }
493
494 private:
495 const Problem* problem_; //!< the problem we are assembling this residual for
496 const TimeLoop* timeLoop_; //!< the timeloop for instationary problems
497 };
498
499 } // end namespace Dumux
500
501 #endif
502