GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/assembly/cvfelocalresidual.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 39 40 97.5%
Functions: 347 395 87.8%
Branches: 36 42 85.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 * \ingroup CVFEDiscretization
11 * \brief Calculates the element-wise residual for control-volume finite element schemes
12 */
13 #ifndef DUMUX_CVFE_LOCAL_RESIDUAL_HH
14 #define DUMUX_CVFE_LOCAL_RESIDUAL_HH
15
16 #include <dune/common/std/type_traits.hh>
17 #include <dune/geometry/type.hh>
18 #include <dune/istl/matrix.hh>
19
20 #include <dumux/common/properties.hh>
21 #include <dumux/common/numeqvector.hh>
22 #include <dumux/assembly/fvlocalresidual.hh>
23 #include <dumux/discretization/extrusion.hh>
24
25 namespace Dumux::Detail {
26
27 template<class Imp, class P, class G, class S, class V>
28 using TimeInfoInterfaceCVFEDetector = decltype(
29 std::declval<Imp>().computeStorage(
30 std::declval<P>(), std::declval<G>(), std::declval<S>(), std::declval<V>(), true
31 )
32 );
33
34 template<class Imp, class P, class G, class S, class V>
35 constexpr inline bool hasTimeInfoInterfaceCVFE()
36 { return Dune::Std::is_detected<TimeInfoInterfaceCVFEDetector, Imp, P, G, S, V>::value; }
37
38 template<class Imp>
39 using SCVFIsOverlappingDetector = decltype(
40 std::declval<Imp>().isOverlapping()
41 );
42
43 template<class Imp>
44 constexpr inline bool hasScvfIsOverlapping()
45 { return Dune::Std::is_detected<SCVFIsOverlappingDetector, Imp>::value; }
46
47 } // end namespace Dumux::Detail
48
49
50 namespace Dumux {
51
52 /*!
53 * \ingroup Assembly
54 * \ingroup CVFEDiscretization
55 * \brief The element-wise residual for control-volume finite element schemes
56 * \tparam TypeTag the TypeTag
57 */
58 template<class TypeTag>
59 class CVFELocalResidual : public FVLocalResidual<TypeTag>
60 {
61 using ParentType = FVLocalResidual<TypeTag>;
62 using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
63 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
64 using Problem = GetPropType<TypeTag, Properties::Problem>;
65 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
66 using GridView = typename GridGeometry::GridView;
67 using Element = typename GridView::template Codim<0>::Entity;
68 using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
69 using FVElementGeometry = typename GridGeometry::LocalView;
70 using GridVolumeVariables = GetPropType<TypeTag, Properties::GridVolumeVariables>;
71 using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
72 using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
73 using SubControlVolume = typename FVElementGeometry::SubControlVolume;
74 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
75 using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
76 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
77 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
78 using Extrusion = Extrusion_t<GridGeometry>;
79
80 public:
81 using ElementResidualVector = typename ParentType::ElementResidualVector;
82
2/4
✓ Branch 1 taken 8653695 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 50002 times.
✗ Branch 5 not taken.
 14461971 using ParentType::ParentType;
83
84 //! evaluate flux residuals for one sub control volume face and add to residual
85 517173535 void evalFlux(ElementResidualVector& residual,
86 const Problem& problem,
87 const Element& element,
88 const FVElementGeometry& fvGeometry,
89 const ElementVolumeVariables& elemVolVars,
90 const ElementBoundaryTypes& elemBcTypes,
91 const ElementFluxVariablesCache& elemFluxVarsCache,
92 const SubControlVolumeFace& scvf) const
93 {
94 524524345 const auto flux = evalFlux(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
95
2/2
✓ Branch 0 taken 417768916 times.
✓ Branch 1 taken 39019355 times.
 517173535 if (!scvf.boundary())
96 {
97
4/4
✓ Branch 0 taken 7400 times.
✓ Branch 1 taken 7400 times.
✓ Branch 2 taken 7400 times.
✓ Branch 3 taken 7400 times.
 961587138 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
98
4/4
✓ Branch 0 taken 7400 times.
✓ Branch 1 taken 7400 times.
✓ Branch 2 taken 7400 times.
✓ Branch 3 taken 7400 times.
 961587138 const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
99
4/4
✓ Branch 0 taken 7400 times.
✓ Branch 1 taken 7400 times.
✓ Branch 2 taken 7400 times.
✓ Branch 3 taken 7400 times.
 961587138 residual[insideScv.localDofIndex()] += flux;
100
101 // for control-volume finite element schemes with overlapping control volumes
102 if constexpr (Detail::hasScvfIsOverlapping<SubControlVolumeFace>())
103 {
104
2/2
✓ Branch 0 taken 18452304 times.
✓ Branch 1 taken 18231704 times.
 69030016 if (!scvf.isOverlapping())
105 69471216 residual[outsideScv.localDofIndex()] -= flux;
106 }
107 else
108 816176296 residual[outsideScv.localDofIndex()] -= flux;
109 }
110 else
111 {
112 80110742 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
113 80110742 residual[insideScv.localDofIndex()] += flux;
114 }
115 517173535 }
116
117 //! evaluate flux residuals for one sub control volume face
118 509799564 NumEqVector evalFlux(const Problem& problem,
119 const Element& element,
120 const FVElementGeometry& fvGeometry,
121 const ElementVolumeVariables& elemVolVars,
122 const ElementBoundaryTypes& elemBcTypes,
123 const ElementFluxVariablesCache& elemFluxVarsCache,
124 const SubControlVolumeFace& scvf) const
125 {
126
2/2
✓ Branch 0 taken 40819720 times.
✓ Branch 1 taken 3492893 times.
 513533142 NumEqVector flux(0.0);
127
128 // inner faces
129
2/2
✓ Branch 0 taken 417768916 times.
✓ Branch 1 taken 39019355 times.
 513533142 if (!scvf.boundary())
130 883022575 flux += this->asImp().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
131
132 // boundary faces
133 else
134 {
135 80110742 const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
136 40055371 const auto& bcTypes = elemBcTypes.get(fvGeometry, scv);
137
138 // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions.
139 // For Dirichlet there is no addition to the residual here but they
140 // are enforced strongly by replacing the residual entry afterwards.
141
4/4
✓ Branch 0 taken 31470912 times.
✓ Branch 1 taken 7548443 times.
✓ Branch 2 taken 31470912 times.
✓ Branch 3 taken 7548443 times.
 80110742 if (bcTypes.hasNeumann())
142 {
143
2/2
✓ Branch 0 taken 56222 times.
✓ Branch 1 taken 56089 times.
 31981876 auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
144
145 // multiply neumann fluxes with the area and the extrusion factor
146 95945628 neumannFluxes *= Extrusion::area(fvGeometry, scvf)*elemVolVars[scv].extrusionFactor();
147
148 // only add fluxes to equations for which Neumann is set
149
2/2
✓ Branch 0 taken 94743863 times.
✓ Branch 1 taken 31470912 times.
 127272979 for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx)
150
2/2
✓ Branch 0 taken 44048 times.
✓ Branch 1 taken 94699815 times.
 95291103 if (bcTypes.isNeumann(eqIdx))
151 278469675 flux[eqIdx] += neumannFluxes[eqIdx];
152 }
153 }
154
155 509799564 return flux;
156 }
157
158 using ParentType::evalStorage;
159 /*!
160 * \brief Compute the storage local residual, i.e. the deviation of the
161 * storage term from zero for instationary problems.
162 *
163 * \param residual The residual vector to fill
164 * \param problem The problem to solve
165 * \param element The DUNE Codim<0> entity for which the residual
166 * ought to be calculated
167 * \param fvGeometry The finite-volume geometry of the element
168 * \param prevElemVolVars The volume averaged variables for all
169 * sub-control volumes of the element at the previous time level
170 * \param curElemVolVars The volume averaged variables for all
171 * sub-control volumes of the element at the current time level
172 * \param scv The sub control volume the storage term is integrated over
173 */
174 337090030 void evalStorage(ElementResidualVector& residual,
175 const Problem& problem,
176 const Element& element,
177 const FVElementGeometry& fvGeometry,
178 const ElementVolumeVariables& prevElemVolVars,
179 const ElementVolumeVariables& curElemVolVars,
180 const SubControlVolume& scv) const
181 {
182
1/2
✓ Branch 0 taken 400000 times.
✗ Branch 1 not taken.
 337090030 const auto& curVolVars = curElemVolVars[scv];
183
1/2
✓ Branch 0 taken 400000 times.
✗ Branch 1 not taken.
 337090030 const auto& prevVolVars = prevElemVolVars[scv];
184
185 // Compute storage with the model specific storage residual
186 // This addresses issues #792/#940 in ad-hoc way by additionally providing crude time level information (previous or current)
187 // to the low-level interfaces if this is supported by the LocalResidual implementation
188 674180060 NumEqVector prevStorage = computeStorageImpl_(problem, fvGeometry, scv, prevVolVars, /*previous time level?*/true);
189 674180060 NumEqVector storage = computeStorageImpl_(problem, fvGeometry, scv, curVolVars, /*previous time level?*/false);
190
191 461209008 prevStorage *= prevVolVars.extrusionFactor();
192 333401390 storage *= curVolVars.extrusionFactor();
193
194 337090030 storage -= prevStorage;
195 674180060 storage *= Extrusion::volume(fvGeometry, scv);
196 337090030 storage /= this->timeLoop().timeStepSize();
197
198 690176956 residual[scv.localDofIndex()] += storage;
199 337090030 }
200
201 private:
202 NumEqVector computeStorageImpl_(const Problem& problem,
203 const FVElementGeometry& fvGeometry,
204 const SubControlVolume& scv,
205 const VolumeVariables& volVars,
206 [[maybe_unused]] bool isPreviousTimeStep) const
207 {
208 if constexpr (Detail::hasTimeInfoInterfaceCVFE<Implementation, Problem, FVElementGeometry, SubControlVolume, VolumeVariables>())
209 ✗ return this->asImp().computeStorage(problem, fvGeometry, scv, volVars, isPreviousTimeStep);
210 else
211
2/4
✓ Branch 3 taken 400000 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 400000 times.
✗ Branch 6 not taken.
 591505266 return this->asImp().computeStorage(problem, scv, volVars);
212 }
213 };
214
215 } // end namespace Dumux
216
217 #endif
218