GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/flux/upwindscheme.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 39 41 95.1%
Functions: 473 589 80.3%
Branches: 39 46 84.8%
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 Flux
10 * \brief Base class for the upwind scheme
11 */
12 #ifndef DUMUX_DISCRETIZATION_UPWINDSCHEME_HH
13 #define DUMUX_DISCRETIZATION_UPWINDSCHEME_HH
14
15 #include <dumux/common/parameters.hh>
16 #include <dumux/discretization/method.hh>
17
18 namespace Dumux {
19
20 //! Forward declaration of the upwind scheme implementation
21 template<class GridGeometry, class DiscretizationMethod>
22 class UpwindSchemeImpl;
23
24 /*!
25 * \ingroup Flux
26 * \brief The upwind scheme used for the advective fluxes.
27 * This depends on the chosen discretization method.
28 */
29 template<class GridGeometry>
30 using UpwindScheme = UpwindSchemeImpl<GridGeometry, typename GridGeometry::DiscretizationMethod>;
31
32 namespace Detail {
33
34 //! returns the upwind factor which is multiplied to the advective flux across the given scvf
35 template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
36 6229258166 Scalar upwindSchemeMultiplier(const ElemVolVars& elemVolVars,
37 const SubControlVolumeFace& scvf,
38 const UpwindTermFunction& upwindTerm,
39 Scalar flux, int phaseIdx)
40 {
41 // TODO: pass this from outside?
42
6/8
✓ Branch 0 taken 1727 times.
✓ Branch 1 taken 3242079895 times.
✓ Branch 3 taken 753 times.
✓ Branch 4 taken 974 times.
✓ Branch 6 taken 753 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 753 times.
✗ Branch 10 not taken.
 6229258166 static const Scalar upwindWeight = getParamFromGroup<Scalar>(elemVolVars.gridVolVars().problem().paramGroup(), "Flux.UpwindWeight");
43
44
4/4
✓ Branch 0 taken 16831789 times.
✓ Branch 1 taken 1427222373 times.
✓ Branch 2 taken 16831789 times.
✓ Branch 3 taken 1427222373 times.
 10819295760 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
45
4/4
✓ Branch 0 taken 16831789 times.
✓ Branch 1 taken 1427222373 times.
✓ Branch 2 taken 632598114 times.
✓ Branch 3 taken 811456048 times.
 9821147512 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
46
47 using std::signbit;
48
4/4
✓ Branch 0 taken 1462835897 times.
✓ Branch 1 taken 1779244995 times.
✓ Branch 2 taken 1462836627 times.
✓ Branch 3 taken 1779244585 times.
 12458514052 if (signbit(flux)) // if sign of flux is negative
49
1/2
✓ Branch 0 taken 18431989 times.
✗ Branch 1 not taken.
 2807182830 return (upwindWeight*upwindTerm(outsideVolVars)
50
1/2
✓ Branch 0 taken 18431989 times.
✗ Branch 1 not taken.
 4395995540 + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
51 else
52
1/2
✓ Branch 0 taken 13906739 times.
✗ Branch 1 not taken.
 3422075380 return (upwindWeight*upwindTerm(insideVolVars)
53
1/2
✓ Branch 0 taken 13906739 times.
✗ Branch 1 not taken.
 5518895840 + (1.0 - upwindWeight)*upwindTerm(outsideVolVars));
54 }
55
56 } // end namespace Detail
57
58 //! Upwind scheme for control-volume finite element methods (uses the standard scheme)
59 template<class GridGeometry, class DM>
60 class UpwindSchemeImpl<GridGeometry, DiscretizationMethods::CVFE<DM>>
61 {
62 public:
63 //! applies a simple upwind scheme to the precalculated advective flux
64 template<class FluxVariables, class UpwindTermFunction, class Scalar>
65 static Scalar apply(const FluxVariables& fluxVars,
66 const UpwindTermFunction& upwindTerm,
67 Scalar flux, int phaseIdx)
68 {
69 2888062244 return apply(fluxVars.elemVolVars(), fluxVars.scvFace(), upwindTerm, flux, phaseIdx);
70 }
71
72 //! applies a simple upwind scheme to the precalculated advective flux across the given scvf
73 template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
74 static Scalar apply(const ElemVolVars& elemVolVars,
75 const SubControlVolumeFace& scvf,
76 const UpwindTermFunction& upwindTerm,
77 Scalar flux, int phaseIdx)
78 {
79 2888062244 return flux*multiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
80 }
81
82 //! returns the upwind factor which is multiplied to the advective flux across the given scvf
83 template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
84 static Scalar multiplier(const ElemVolVars& elemVolVars,
85 const SubControlVolumeFace& scvf,
86 const UpwindTermFunction& upwindTerm,
87 Scalar flux, int phaseIdx)
88 {
89 1444054162 return Detail::upwindSchemeMultiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
90 }
91 };
92
93 //! Upwind scheme for the cell-centered tpfa scheme
94 template<class GridGeometry>
95 class UpwindSchemeImpl<GridGeometry, DiscretizationMethods::CCTpfa>
96 {
97 using GridView = typename GridGeometry::GridView;
98 static constexpr int dim = GridView::dimension;
99 static constexpr int dimWorld = GridView::dimensionworld;
100
101 public:
102 //! returns the upwind factor which is multiplied to the advective flux across the given scvf
103 template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
104 ✗ static Scalar multiplier(const ElemVolVars& elemVolVars,
105 const SubControlVolumeFace& scvf,
106 const UpwindTermFunction& upwindTerm,
107 Scalar flux, int phaseIdx)
108 {
109 static_assert(dim == dimWorld, "Multiplier cannot be computed on surface grids using cell-centered scheme!");
110 1468881570 return Detail::upwindSchemeMultiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
111 }
112
113 //! applies a simple upwind scheme to the precalculated advective flux across the given scvf
114 template<class ElemVolVars, class SubControlVolumeFace, class UpwindTermFunction, class Scalar>
115 ✗ static Scalar apply(const ElemVolVars& elemVolVars,
116 const SubControlVolumeFace& scvf,
117 const UpwindTermFunction& upwindTerm,
118 Scalar flux, int phaseIdx)
119 {
120 static_assert(dim == dimWorld, "This upwind scheme cannot be used on surface grids using cell-centered schemes!");
121 2937707112 return flux*multiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
122 }
123
124 // For surface and network grids (dim < dimWorld) we have to do a special upwind scheme
125 template<class FluxVariables, class UpwindTermFunction, class Scalar>
126 320742313 static Scalar apply(const FluxVariables& fluxVars,
127 const UpwindTermFunction& upwindTerm,
128 Scalar flux, int phaseIdx)
129 {
130 1803588943 const auto& scvf = fluxVars.scvFace();
131 1803588943 const auto& elemVolVars = fluxVars.elemVolVars();
132
133 // standard scheme
134 if constexpr (dim == dimWorld)
135 1482846630 return apply(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
136
137 // on non-branching points the standard scheme works
138
4/4
✓ Branch 0 taken 8952144 times.
✓ Branch 1 taken 523679 times.
✓ Branch 2 taken 8952144 times.
✓ Branch 3 taken 523679 times.
 641484626 if (scvf.numOutsideScvs() == 1)
139 315146726 return flux*Detail::upwindSchemeMultiplier(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
140
141 // handle branching points with a more complicated upwind scheme
142 // we compute a flux-weighted average of all inflowing branches
143 6842602 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
144
145 5595587 Scalar branchingPointUpwindTerm = 0.0;
146 5595587 Scalar sumUpwindFluxes = 0.0;
147
148 // if the inside flux is positive (outflow) do fully upwind and return flux
149 using std::signbit;
150
4/4
✓ Branch 0 taken 386380 times.
✓ Branch 1 taken 137299 times.
✓ Branch 2 taken 386380 times.
✓ Branch 3 taken 137299 times.
 11191174 if (!signbit(flux))
151 4798918 return upwindTerm(insideVolVars)*flux;
152 else
153 3159984 sumUpwindFluxes += flux;
154
155
4/4
✓ Branch 0 taken 426066 times.
✓ Branch 1 taken 137299 times.
✓ Branch 2 taken 426066 times.
✓ Branch 3 taken 137299 times.
 12639866 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
156 {
157 // compute the outside flux
158 9479882 const auto& fvGeometry = fluxVars.fvGeometry();
159 9479882 const auto outsideScvIdx = scvf.outsideScvIdx(i);
160 9479882 const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx);
161 9479882 const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
162
163 using AdvectionType = typename FluxVariables::AdvectionType;
164 9479882 const auto outsideFlux = AdvectionType::flux(fluxVars.problem(),
165 outsideElement,
166 fvGeometry,
167 elemVolVars,
168 flippedScvf,
169 phaseIdx,
170 fluxVars.elemFluxVarsCache());
171
172
4/4
✓ Branch 0 taken 400936 times.
✓ Branch 1 taken 25130 times.
✓ Branch 2 taken 400936 times.
✓ Branch 3 taken 25130 times.
 18959764 if (!signbit(outsideFlux))
173 4791808 branchingPointUpwindTerm += upwindTerm(elemVolVars[outsideScvIdx])*outsideFlux;
174 else
175 4688074 sumUpwindFluxes += outsideFlux;
176 }
177
178 // the flux might be zero
179
1/2
✓ Branch 0 taken 137299 times.
✗ Branch 1 not taken.
 3159984 if (sumUpwindFluxes != 0.0)
180 3159860 branchingPointUpwindTerm /= -sumUpwindFluxes;
181 else
182 branchingPointUpwindTerm = 0.0;
183
184 // upwind scheme (always do fully upwind at branching points)
185 // a weighting here would lead to an error since the derivation is based on a fully upwind scheme
186 // TODO How to implement a weight of e.g. 0.5
187 3159984 if (signbit(flux))
188 3159984 return flux*branchingPointUpwindTerm;
189 else
190 return flux*upwindTerm(insideVolVars);
191 }
192 };
193
194 //! Upwind scheme for cell-centered mpfa schemes
195 template<class GridGeometry>
196 class UpwindSchemeImpl<GridGeometry, DiscretizationMethods::CCMpfa>
197 : public UpwindSchemeImpl<GridGeometry, DiscretizationMethods::CCTpfa> {};
198
199 } // end namespace Dumux
200
201 #endif
202