GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/flux/upwindscheme.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 39 41 95.1%
Functions: 467 583 80.1%
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 6499678710 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 1739 times.
✓ Branch 1 taken 3378631463 times.
✓ Branch 3 taken 755 times.
✓ Branch 4 taken 984 times.
✓ Branch 6 taken 755 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 755 times.
✗ Branch 10 not taken.
 6499678710 static const Scalar upwindWeight = getParamFromGroup<Scalar>(elemVolVars.gridVolVars().problem().paramGroup(), "Flux.UpwindWeight");
43
44
4/4
✓ Branch 0 taken 16831741 times.
✓ Branch 1 taken 1433307689 times.
✓ Branch 2 taken 16831741 times.
✓ Branch 3 taken 1433307689 times.
 11360136848 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
45
4/4
✓ Branch 0 taken 16831741 times.
✓ Branch 1 taken 1433307689 times.
✓ Branch 2 taken 634236055 times.
✓ Branch 3 taken 815903375 times.
 10351294412 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
46
47 using std::signbit;
48
4/4
✓ Branch 0 taken 1529335429 times.
✓ Branch 1 taken 1849297043 times.
✓ Branch 2 taken 1529336159 times.
✓ Branch 3 taken 1849296633 times.
 12999355140 if (signbit(flux)) // if sign of flux is negative
49
1/2
✓ Branch 0 taken 18439011 times.
✗ Branch 1 not taken.
 2939520194 return (upwindWeight*upwindTerm(outsideVolVars)
50
1/2
✓ Branch 0 taken 18439011 times.
✗ Branch 1 not taken.
 4656325796 + (1.0 - upwindWeight)*upwindTerm(insideVolVars));
51 else
52
1/2
✓ Branch 0 taken 13942692 times.
✗ Branch 1 not taken.
 3560158560 return (upwindWeight*upwindTerm(insideVolVars)
53
1/2
✓ Branch 0 taken 13942692 times.
✗ Branch 1 not taken.
 5788357552 + (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 2900232780 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 2900232780 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 1450139430 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 1593009670 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 3185963312 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 327415837 static Scalar apply(const FluxVariables& fluxVars,
127 const UpwindTermFunction& upwindTerm,
128 Scalar flux, int phaseIdx)
129 {
130 1934390567 const auto& scvf = fluxVars.scvFace();
131 1934390567 const auto& elemVolVars = fluxVars.elemVolVars();
132
133 // standard scheme
134 if constexpr (dim == dimWorld)
135 1606974730 return apply(elemVolVars, scvf, upwindTerm, flux, phaseIdx);
136
137 // on non-branching points the standard scheme works
138
4/4
✓ Branch 0 taken 15290356 times.
✓ Branch 1 taken 858991 times.
✓ Branch 2 taken 15290356 times.
✓ Branch 3 taken 858991 times.
 654831674 if (scvf.numOutsideScvs() == 1)
139 321484938 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 7513226 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
144
145 5930899 Scalar branchingPointUpwindTerm = 0.0;
146 5930899 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 638946 times.
✓ Branch 1 taken 220045 times.
✓ Branch 2 taken 638946 times.
✓ Branch 3 taken 220045 times.
 11861798 if (!signbit(flux))
151 5304050 return upwindTerm(insideVolVars)*flux;
152 else
153 3242730 sumUpwindFluxes += flux;
154
155
4/4
✓ Branch 0 taken 710563 times.
✓ Branch 1 taken 220045 times.
✓ Branch 2 taken 710563 times.
✓ Branch 3 taken 220045 times.
 13007109 for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
156 {
157 // compute the outside flux
158 9764379 const auto& fvGeometry = fluxVars.fvGeometry();
159 9764379 const auto outsideScvIdx = scvf.outsideScvIdx(i);
160 9764379 const auto outsideElement = fvGeometry.gridGeometry().element(outsideScvIdx);
161 9764379 const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
162
163 using AdvectionType = typename FluxVariables::AdvectionType;
164 9764379 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 667991 times.
✓ Branch 1 taken 42572 times.
✓ Branch 2 taken 667991 times.
✓ Branch 3 taken 42572 times.
 19528758 if (!signbit(outsideFlux))
173 5058863 branchingPointUpwindTerm += upwindTerm(elemVolVars[outsideScvIdx])*outsideFlux;
174 else
175 4705516 sumUpwindFluxes += outsideFlux;
176 }
177
178 // the flux might be zero
179
1/2
✓ Branch 0 taken 220045 times.
✗ Branch 1 not taken.
 3242730 if (sumUpwindFluxes != 0.0)
180 3242606 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 3242730 if (signbit(flux))
188 3242730 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