GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/assembly/fvlocalresidual.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 37 42 88.1%
Functions: 1240 1909 65.0%
Branches: 30 61 49.2%
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 277435652 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 213847769 times.
 277435652 assert(timeLoop_ && "no time loop set for storage term evaluation");
137
138 // initialize the residual vector for all scvs in this element
139 552386696 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 471906567 times.
✓ Branch 1 taken 217492181 times.
✓ Branch 2 taken 439773322 times.
✓ Branch 3 taken 185445156 times.
✓ Branch 4 taken 1214804 times.
 1120751664 for (auto&& scv : scvs(fvGeometry))
144 568309968 asImp().evalStorage(residual, this->problem(), element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
145
146 277435652 return residual;
147 }
148
149 341278805 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 677338702 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 554724951 times.
✓ Branch 1 taken 260516605 times.
✓ Branch 2 taken 518498976 times.
✓ Branch 3 taken 226667707 times.
✓ Branch 4 taken 2581782 times.
 1394568642 for (auto&& scv : scvs(fvGeometry))
161 717167870 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 1096332334 times.
✓ Branch 1 taken 281392071 times.
✓ Branch 2 taken 1093672121 times.
✓ Branch 3 taken 276495284 times.
✓ Branch 4 taken 2581782 times.
 3077605561 for (auto&& scvf : scvfs(fvGeometry))
165 1834345539 asImp().evalFlux(residual, this->problem(), element, fvGeometry, elemVolVars, bcTypes, elemFluxVarsCache, scvf);
166
167 341278803 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 523700340 NumEqVector computeSource(const Problem& problem,
209 const Element& element,
210 const FVElementGeometry& fvGeometry,
211 const ElementVolumeVariables& elemVolVars,
212 const SubControlVolume &scv) const
213 {
214
3/5
✓ Branch 0 taken 17400 times.
✓ Branch 1 taken 245400 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 148 times.
✗ Branch 5 not taken.
 593723509 NumEqVector source(0.0);
215
216 // add contributions from volume flux sources
217
4/16
✓ Branch 0 taken 17400 times.
✓ Branch 1 taken 245400 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 148 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 148 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 940925151 source += problem.source(element, fvGeometry, elemVolVars, scv);
218
219 // add contribution from possible point sources
220
1/7
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 148 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 602297165 source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv);
221
222 523700338 return source;
223 }
224
225 /*!
226 * \brief Calculate the flux term of the equation
227 *
228 * \param problem The problem to solve
229 * \param element The DUNE Codim<0> entity for which the residual
230 * ought to be calculated
231 * \param fvGeometry The finite-volume geometry of the element
232 * \param elemVolVars The volume variables associated with the element stencil
233 * \param scvf The sub-control volume over which we integrate the flux
234 * \param elemFluxVarsCache the flux variable caches for the element's flux stencils
235 *
236 * \note has to be implemented by the model specific residual class
237 *
238 */
239 NumEqVector computeFlux(const Problem& problem,
240 const Element& element,
241 const FVElementGeometry& fvGeometry,
242 const ElementVolumeVariables& elemVolVars,
243 const SubControlVolumeFace& scvf,
244 const ElementFluxVariablesCache& elemFluxVarsCache) const
245 {
246 DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!");
247 }
248
249 // \}
250
251 /*!
252 * \name Discretization specific interface
253 * \note The following method are the discretization specific wrapper methods
254 */
255 // \{
256
257 /*!
258 * \brief Compute the storage local residual, i.e. the deviation of the
259 * storage term from zero for instationary problems.
260 *
261 * \param residual The residual vector to fill
262 * \param problem The problem to solve
263 * \param element The DUNE Codim<0> entity for which the residual
264 * ought to be calculated
265 * \param fvGeometry The finite-volume geometry of the element
266 * \param prevElemVolVars The volume averaged variables for all
267 * sub-control volumes of the element at the previous time level
268 * \param curElemVolVars The volume averaged variables for all
269 * sub-control volumes of the element at the current time level
270 * \param scv The sub control volume the storage term is integrated over
271 */
272 191974622 void evalStorage(ElementResidualVector& residual,
273 const Problem& problem,
274 const Element& element,
275 const FVElementGeometry& fvGeometry,
276 const ElementVolumeVariables& prevElemVolVars,
277 const ElementVolumeVariables& curElemVolVars,
278 const SubControlVolume& scv) const
279 {
280
1/2
✓ Branch 0 taken 14460000 times.
✗ Branch 1 not taken.
 191974622 const auto& curVolVars = curElemVolVars[scv];
281
1/2
✓ Branch 0 taken 14460000 times.
✗ Branch 1 not taken.
 191974622 const auto& prevVolVars = prevElemVolVars[scv];
282
283 // mass balance within the element. this is the
284 // \f$\frac{m}{\partial t}\f$ term if using implicit or explicit
285 // euler as time discretization.
286 //
287 // TODO: We might need a more explicit way for
288 // doing the time discretization...
289
290 //! Compute storage with the model specific storage residual
291
1/2
✓ Branch 0 taken 21870800 times.
✗ Branch 1 not taken.
 191974622 NumEqVector prevStorage = asImp().computeStorage(problem, scv, prevVolVars);
292
1/2
✓ Branch 0 taken 21870800 times.
✗ Branch 1 not taken.
 191974622 NumEqVector storage = asImp().computeStorage(problem, scv, curVolVars);
293
294 191974622 prevStorage *= prevVolVars.extrusionFactor();
295 170645327 storage *= curVolVars.extrusionFactor();
296
297 191974622 storage -= prevStorage;
298 383949244 storage *= Extrusion::volume(fvGeometry, scv);
299 191974622 storage /= timeLoop_->timeStepSize();
300
301 454865794 residual[scv.localDofIndex()] += storage;
302 191974622 }
303
304 /*!
305 * \brief Compute the source local residual, i.e. the deviation of the
306 * source term from zero.
307 *
308 * \param residual The residual vector to fill
309 * \param problem The problem to solve
310 * \param element The DUNE Codim<0> entity for which the residual
311 * ought to be calculated
312 * \param fvGeometry The finite-volume geometry of the element
313 * \param curElemVolVars The volume averaged variables for all
314 * sub-control volumes of the element at the current time level
315 * \param scv The sub control volume the source term is integrated over
316 */
317 772920238 void evalSource(ElementResidualVector& residual,
318 const Problem& problem,
319 const Element& element,
320 const FVElementGeometry& fvGeometry,
321 const ElementVolumeVariables& curElemVolVars,
322 const SubControlVolume& scv) const
323 {
324 //! Compute source with the model specific storage residual
325
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 772936078 const auto& curVolVars = curElemVolVars[scv];
326
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 772936078 NumEqVector source = asImp().computeSource(problem, element, fvGeometry, curElemVolVars, scv);
327 1839232478 source *= Extrusion::volume(fvGeometry, scv)*curVolVars.extrusionFactor();
328
329 //! subtract source from local rate (sign convention in user interface)
330 1767195710 residual[scv.localDofIndex()] -= source;
331 772920236 }
332
333 /*!
334 * \brief Compute the flux local residual, i.e. the deviation of the
335 * flux term from zero.
336 * \param residual The residual vector to fill
337 * \param problem The problem to solve
338 * \param element The DUNE Codim<0> entity for which the residual
339 * ought to be calculated
340 * \param fvGeometry The finite-volume geometry of the element
341 * \param elemVolVars The volume averaged variables for all
342 * sub-control volumes of the element at the current time level
343 * \param elemBcTypes the boundary types for the boundary entities of an elements
344 * \param elemFluxVarsCache The flux variable caches for the element stencil
345 * \param scvf The sub control volume face the flux term is integrated over
346 */
347 void evalFlux(ElementResidualVector& residual,
348 const Problem& problem,
349 const Element& element,
350 const FVElementGeometry& fvGeometry,
351 const ElementVolumeVariables& elemVolVars,
352 const ElementBoundaryTypes& elemBcTypes,
353 const ElementFluxVariablesCache& elemFluxVarsCache,
354 const SubControlVolumeFace& scvf) const {}
355
356 /*!
357 * \brief Compute the flux local residual, i.e. the deviation of the
358 * flux term from zero.
359 *
360 * \param problem The problem to solve
361 * \param element The DUNE Codim<0> entity for which the residual
362 * ought to be calculated
363 * \param fvGeometry The finite-volume geometry of the element
364 * \param elemVolVars The volume averaged variables for all
365 * sub-control volumes of the element at the current time level
366 * \param elemFluxVarsCache The flux variable caches for the element stencil
367 * \param scvf The sub control volume face the flux term is integrated over
368 */
369 NumEqVector evalFlux(const Problem& problem,
370 const Element& element,
371 const FVElementGeometry& fvGeometry,
372 const ElementVolumeVariables& elemVolVars,
373 const ElementFluxVariablesCache& elemFluxVarsCache,
374 const SubControlVolumeFace& scvf) const
375 {
376 return asImp().evalFlux(problem, element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
377 }
378
379 //\}
380
381 /*!
382 * \name Interfaces for analytic Jacobian computation
383 */
384 // \{
385
386 //! Compute the derivative of the storage residual
387 template<class PartialDerivativeMatrix>
388 void addStorageDerivatives(PartialDerivativeMatrix& partialDerivatives,
389 const Problem& problem,
390 const Element& element,
391 const FVElementGeometry& fvGeometry,
392 const VolumeVariables& curVolVars,
393 const SubControlVolume& scv) const
394 {
395 DUNE_THROW(Dune::NotImplemented, "analytic storage derivative");
396 }
397
398 //! Compute the derivative of the source residual
399 template<class PartialDerivativeMatrix>
400 void addSourceDerivatives(PartialDerivativeMatrix& partialDerivatives,
401 const Problem& problem,
402 const Element& element,
403 const FVElementGeometry& fvGeometry,
404 const VolumeVariables& curVolVars,
405 const SubControlVolume& scv) const
406 {
407 DUNE_THROW(Dune::NotImplemented, "analytic source derivative");
408 }
409
410 //! Compute the derivative of the flux residual
411 template<class PartialDerivativeMatrices, class T = TypeTag>
412 std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod != DiscretizationMethods::box, void>
413 addFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
414 const Problem& problem,
415 const Element& element,
416 const FVElementGeometry& fvGeometry,
417 const ElementVolumeVariables& curElemVolVars,
418 const ElementFluxVariablesCache& elemFluxVarsCache,
419 const SubControlVolumeFace& scvf) const
420 {
421 DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for cell-centered models");
422 }
423
424 //! Compute the derivative of the flux residual for the box method
425 template<class JacobianMatrix, class T = TypeTag>
426 std::enable_if_t<GetPropType<T, Properties::GridGeometry>::discMethod == DiscretizationMethods::box, void>
427 addFluxDerivatives(JacobianMatrix& A,
428 const Problem& problem,
429 const Element& element,
430 const FVElementGeometry& fvGeometry,
431 const ElementVolumeVariables& curElemVolVars,
432 const ElementFluxVariablesCache& elemFluxVarsCache,
433 const SubControlVolumeFace& scvf) const
434 {
435 DUNE_THROW(Dune::NotImplemented, "analytic flux derivative for box models");
436 }
437
438 //! Compute the derivative of the Dirichlet flux residual for cell-centered schemes
439 template<class PartialDerivativeMatrices>
440 void addCCDirichletFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
441 const Problem& problem,
442 const Element& element,
443 const FVElementGeometry& fvGeometry,
444 const ElementVolumeVariables& curElemVolVars,
445 const ElementFluxVariablesCache& elemFluxVarsCache,
446 const SubControlVolumeFace& scvf) const
447 {
448 DUNE_THROW(Dune::NotImplemented, "analytic Dirichlet flux derivative");
449 }
450
451 //! Compute the derivative of Robin type boundary conditions ("solution dependent Neumann")
452 template<class PartialDerivativeMatrices>
453 void addRobinFluxDerivatives(PartialDerivativeMatrices& derivativeMatrices,
454 const Problem& problem,
455 const Element& element,
456 const FVElementGeometry& fvGeometry,
457 const ElementVolumeVariables& curElemVolVars,
458 const ElementFluxVariablesCache& elemFluxVarsCache,
459 const SubControlVolumeFace& scvf) const
460 {
461 DUNE_THROW(Dune::NotImplemented, "analytic Robin flux derivative");
462 }
463
464 //\}
465
466 /*!
467 * \name Interfaces accessed by local residual implementations
468 */
469 // \{
470
471 //! the problem
472 ✗ const Problem& problem() const
473 ✗ { return *problem_; }
474
475 //! the timeloop for instationary problems
476 //! calling this for stationary leads to undefined behaviour
477 ✗ const TimeLoop& timeLoop() const
478 ✗ { return *timeLoop_; }
479
480 //! returns true if the residual is stationary
481 ✗ bool isStationary() const
482
2/4
✓ Branch 1 taken 7 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 9 { return !timeLoop_; }
483
484 // \}
485 protected:
486
487 Implementation &asImp()
488 { return *static_cast<Implementation*>(this); }
489
490 const Implementation &asImp() const
491 { return *static_cast<const Implementation*>(this); }
492
493 private:
494 const Problem* problem_; //!< the problem we are assembling this residual for
495 const TimeLoop* timeLoop_; //!< the timeloop for instationary problems
496 };
497
498 } // end namespace Dumux
499
500 #endif
501