GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: dumux/dumux/assembly/fclocalresidual.hh
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 29 30 96.7%
Functions: 86 102 84.3%
Branches: 12 17 70.6%

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 Assembly
10 * \ingroup FaceCenteredStaggeredDiscretization
11 * \brief Calculates the element-wise residual for the box scheme
12 */
13 #ifndef DUMUX_FACECENTERED_LOCAL_RESIDUAL_HH
14 #define DUMUX_FACECENTERED_LOCAL_RESIDUAL_HH
15
16 #include <dune/geometry/type.hh>
17 #include <dune/istl/matrix.hh>
18
19 #include <dumux/common/numeqvector.hh>
20 #include <dumux/common/properties.hh>
21 #include <dumux/assembly/fvlocalresidual.hh>
22 #include <dumux/discretization/extrusion.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \ingroup Assembly
28 * \ingroup FaceCenteredStaggeredDiscretization
29 * \brief The element-wise residual for the box scheme
30 * \tparam TypeTag the TypeTag
31 */
32 template<class TypeTag>
33 class FaceCenteredLocalResidual : public FVLocalResidual<TypeTag>
34 {
35 using ParentType = FVLocalResidual<TypeTag>;
36 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
37 using Problem = GetPropType<TypeTag, Properties::Problem>;
38 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
39 using GridView = typename GridGeometry::GridView;
40 using Element = typename GridView::template Codim<0>::Entity;
41 using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
42 using FVElementGeometry = typename GridGeometry::LocalView;
43 using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
44 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
45 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
46 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
47 using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
48 using Extrusion = Extrusion_t<GridGeometry>;
49
50 public:
51 using ElementResidualVector = typename ParentType::ElementResidualVector;
52
1/2
✓ Branch 1 taken 2730761 times.
✗ Branch 2 not taken.
2761797 using ParentType::ParentType;
53
54 /*!
55 * \brief evaluate flux residuals for one sub control volume face and add to residual
56 *
57 */
58 315842860 void evalFlux(ElementResidualVector& residual,
59 const Problem& problem,
60 const Element& element,
61 const FVElementGeometry& fvGeometry,
62 const ElementVolumeVariables& elemVolVars,
63 const ElementBoundaryTypes& elemBcTypes,
64 const ElementFluxVariablesCache& elemFluxVarsCache,
65 const SubControlVolumeFace& scvf) const
66 {
67 315842860 const auto flux = evalFlux(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
68
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 309484660 times.
315842860 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
69 315842860 residual[insideScv.localDofIndex()] += flux;
70 315842860 }
71
72 //! evaluate flux residuals for one sub control volume face
73
2/2
✓ Branch 0 taken 21094298 times.
✓ Branch 1 taken 294748562 times.
315842860 NumEqVector evalFlux(const Problem& problem,
74 const Element& element,
75 const FVElementGeometry& fvGeometry,
76 const ElementVolumeVariables& elemVolVars,
77 const ElementBoundaryTypes& elemBcTypes,
78 const ElementFluxVariablesCache& elemFluxVarsCache,
79 const SubControlVolumeFace& scvf) const
80 {
81
2/2
✓ Branch 0 taken 21094298 times.
✓ Branch 1 taken 294748562 times.
315842860 if (elemBcTypes.hasDirichlet())
82 21094298 return this->asImp().maybeHandleDirichletBoundary(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
83
84
2/2
✓ Branch 0 taken 4687524 times.
✓ Branch 1 taken 290061038 times.
294748562 if (elemBcTypes.hasNeumann())
85 4687524 return this->asImp().maybeHandleNeumannBoundary(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
86
87 290061038 return this->asImp().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache, elemBcTypes);
88 }
89
90 /*!
91 * \brief Calculate the source term of the equation
92 *
93 * \param problem The problem to solve
94 * \param element The DUNE Codim<0> entity for which the residual
95 * ought to be calculated
96 * \param fvGeometry The finite-volume geometry of the element
97 * \param elemVolVars The volume variables associated with the element stencil
98 * \param scv The sub-control volume over which we integrate the source term
99 * \note This is the default implementation for all models as sources are computed
100 * in the user interface of the problem
101 *
102 */
103 104220774 NumEqVector computeSource(const Problem& problem,
104 const Element& element,
105 const FVElementGeometry& fvGeometry,
106 const ElementVolumeVariables& elemVolVars,
107 const SubControlVolume& scv) const
108 {
109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 401232 times.
104220774 NumEqVector source(0.0);
110
111 // add contributions from volume flux sources
112
2/3
✗ Branch 0 not taken.
✓ Branch 1 taken 86030398 times.
✓ Branch 2 taken 18190376 times.
150070644 source += problem.source(element, fvGeometry, elemVolVars, scv)[scv.dofAxis()];
113
114 // add contribution from possible point sources
115
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 104220774 times.
104220774 if (!problem.pointSourceMap().empty())
116 source += problem.scvPointSources(element, fvGeometry, elemVolVars, scv)[scv.dofAxis()];
117
118 103819542 return source;
119 }
120
121 using ParentType::evalStorage;
122
123 /*!
124 * \brief Compute the storage local residual, i.e. the deviation of the
125 * storage term from zero for instationary problems.
126 *
127 * \param residual The residual vector to fill
128 * \param problem The problem to solve
129 * \param element The DUNE Codim<0> entity for which the residual
130 * ought to be calculated
131 * \param fvGeometry The finite-volume geometry of the element
132 * \param prevElemVolVars The volume averaged variables for all
133 * sub-control volumes of the element at the previous time level
134 * \param curElemVolVars The volume averaged variables for all
135 * sub-control volumes of the element at the current time level
136 * \param scv The sub control volume the storage term is integrated over
137 */
138 59953376 void evalStorage(ElementResidualVector& residual,
139 const Problem& problem,
140 const Element& element,
141 const FVElementGeometry& fvGeometry,
142 const ElementVolumeVariables& prevElemVolVars,
143 const ElementVolumeVariables& curElemVolVars,
144 const SubControlVolume& scv) const
145 {
146 59953376 const auto& curVolVars = curElemVolVars[scv];
147 59953376 const auto& prevVolVars = prevElemVolVars[scv];
148
149 // mass balance within the element. this is the
150 // \f$\frac{m}{\partial t}\f$ term if using implicit or explicit
151 // euler as time discretization.
152 //
153
154 //! Compute storage with the model specific storage residual
155 59953376 NumEqVector prevStorage = this->asImp().computeStorage(problem, scv, prevVolVars, true/*isPreviousStorage*/);
156 59953376 NumEqVector storage = this->asImp().computeStorage(problem, scv, curVolVars, false/*isPreviousStorage*/);
157
158 59953376 prevStorage *= prevVolVars.extrusionFactor();
159 119906752 storage *= curVolVars.extrusionFactor();
160
161 59953376 storage -= prevStorage;
162 59953376 storage *= Extrusion::volume(fvGeometry, scv);
163 59953376 storage /= this->timeLoop().timeStepSize();
164
165 59953376 residual[scv.localDofIndex()] += storage;
166 59953376 }
167 };
168
169 } // end namespace Dumux
170
171 #endif
172