GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/test/porousmediumflow/3pwateroil/problem.hh
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 50 54 92.6%
Functions: 3 4 75.0%
Branches: 73 118 61.9%
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 ThreePWaterOilTests
10 * \brief Non-isothermal steam-assisted gravity
11 * drainage (SAGD) problem.
12 */
13
14 #ifndef DUMUX_SAGDPROBLEM_HH
15 #define DUMUX_SAGDPROBLEM_HH
16
17 #include <dumux/common/properties.hh>
18 #include <dumux/common/parameters.hh>
19 #include <dumux/common/boundarytypes.hh>
20 #include <dumux/common/numeqvector.hh>
21
22 #include <dumux/porousmediumflow/problem.hh>
23
24 namespace Dumux {
25
26 /*!
27 * \file
28 * \ingroup ThreePWaterOilTests
29 * \brief Non-isothermal steam-assisted gravity
30 * drainage (SAGD) problem.
31 */
32 template <class TypeTag >
33 class SagdProblem : public PorousMediumFlowProblem<TypeTag>
34 {
35 using ParentType = PorousMediumFlowProblem<TypeTag>;
36 using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
37 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
38 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
39 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
40 using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
41 enum {
42 pressureIdx = Indices::pressureIdx,
43 switch1Idx = Indices::switch1Idx,
44 switch2Idx = Indices::switch2Idx,
45
46 contiWEqIdx = Indices::conti0EqIdx + FluidSystem::wCompIdx,
47 contiNEqIdx = Indices::conti0EqIdx + FluidSystem::nCompIdx,
48 energyEqIdx = Indices::energyEqIdx,
49
50 // phase indices
51 wPhaseIdx = FluidSystem::wPhaseIdx,
52 nPhaseIdx = FluidSystem::nPhaseIdx,
53
54 // phase state
55 wnPhaseOnly = Indices::wnPhaseOnly,
56
57 // world dimension
58 dimWorld = GridView::dimensionworld
59 };
60
61 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
62 using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
63 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
64 using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
65 using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
66 using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
67 using Element = typename GridView::template Codim<0>::Entity;
68 using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
69 using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
70 using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
71
72 public:
73
74 1 SagdProblem(std::shared_ptr<const GridGeometry> gridGeometry)
75
7/20
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
 3 : ParentType(gridGeometry), pOut_(4e6)
76 {
77 1 maxDepth_ = 400.0; // [m]
78
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 FluidSystem::init();
79 1 }
80
81 /*!
82 * \name Boundary conditions
83 */
84 // \{
85
86 /*!
87 * \brief Specifies which kind of boundary condition should be
88 * used for which equation on a given boundary control volume.
89 *
90 * \param globalPos The position for which the boundary types are evaluated
91 */
92 45874 BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
93 {
94 45874 BoundaryTypes bcTypes;
95 // on bottom
96
4/4
✓ Branch 0 taken 13264 times.
✓ Branch 1 taken 32610 times.
✓ Branch 2 taken 13264 times.
✓ Branch 3 taken 32610 times.
 91748 if (globalPos[1] < eps_)
97 bcTypes.setAllNeumann();
98
99 // on top
100
4/4
✓ Branch 0 taken 13142 times.
✓ Branch 1 taken 19468 times.
✓ Branch 2 taken 13142 times.
✓ Branch 3 taken 19468 times.
 65220 else if (globalPos[1] > 40.0 - eps_)
101 bcTypes.setAllNeumann();
102
103 // on bottom other than corners
104
4/4
✓ Branch 0 taken 8518 times.
✓ Branch 1 taken 10950 times.
✓ Branch 2 taken 8518 times.
✓ Branch 3 taken 10950 times.
 38936 else if (globalPos[0] > 60 - eps_ )
105 bcTypes.setAllDirichlet();
106
107 // on Left
108 else
109 bcTypes.setAllNeumann();
110 45874 return bcTypes;
111 }
112
113 /*!
114 * \brief Evaluates the boundary conditions for a Dirichlet control volume.
115 *
116 * \param globalPos The center of the finite volume which ought to be set.
117 */
118 PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
119 {
120 12168 return initial_(globalPos);
121 }
122
123 /*!
124 * \brief Evaluates the boundary conditions for a Neumann boundary segment.
125 *
126 * \param element The finite element
127 * \param fvGeometry The finite-volume geometry in the box scheme
128 * \param elemVolVars The element volume variables
129 * \param elemFluxVarsCache Flux variables caches for all faces in stencil
130 * \param scvf The sub-control volume face
131 *
132 * Negative values mean influx.
133 */
134 267904 NumEqVector neumann(const Element& element,
135 const FVElementGeometry& fvGeometry,
136 const ElementVolumeVariables& elemVolVars,
137 const ElementFluxVariablesCache& elemFluxVarsCache,
138 const SubControlVolumeFace& scvf) const
139 {
140 267904 NumEqVector values(0.0);
141
142
4/4
✓ Branch 0 taken 153088 times.
✓ Branch 1 taken 114816 times.
✓ Branch 2 taken 153088 times.
✓ Branch 3 taken 114816 times.
 535808 const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
143
2/2
✓ Branch 0 taken 153088 times.
✓ Branch 1 taken 114816 times.
 267904 const auto& globalPos = insideScv.dofPosition();
144
145 // negative values for injection at injection well
146
8/8
✓ Branch 0 taken 153088 times.
✓ Branch 1 taken 114816 times.
✓ Branch 2 taken 153088 times.
✓ Branch 3 taken 114816 times.
✓ Branch 4 taken 1664 times.
✓ Branch 5 taken 151424 times.
✓ Branch 6 taken 1664 times.
✓ Branch 7 taken 151424 times.
 535808 if (globalPos[1] > 8.5 - eps_ && globalPos[1] < 9.5 + eps_)
147 {
148 1664 values[contiNEqIdx] = -0.0;
149 1664 values[contiWEqIdx] = -0.193;//*0.5; // (55.5 mol*12.5)/3600 mol/s m = 0.193
150 3328 values[energyEqIdx] = -9132;//*0.5; // J/sec m 9132
151 }
152
8/8
✓ Branch 0 taken 161408 times.
✓ Branch 1 taken 104832 times.
✓ Branch 2 taken 161408 times.
✓ Branch 3 taken 104832 times.
✓ Branch 4 taken 1664 times.
✓ Branch 5 taken 159744 times.
✓ Branch 6 taken 1664 times.
✓ Branch 7 taken 159744 times.
 532480 else if (globalPos[1] > 2.5 - eps_ && globalPos[1] < 3.5 + eps_) // production well
153 {
154
155
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
 4992 const Scalar elemPressW = elemVolVars[scvf.insideScvIdx()].pressure(wPhaseIdx); //Pressures
156
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
 4992 const Scalar elemPressN = elemVolVars[scvf.insideScvIdx()].pressure(nPhaseIdx);
157
158
5/10
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1664 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1664 times.
 8320 const Scalar densityW = elemVolVars[scvf.insideScvIdx()].fluidState().density(wPhaseIdx); //Densities
159
4/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1664 times.
 6656 const Scalar densityN = elemVolVars[scvf.insideScvIdx()].fluidState().density(nPhaseIdx);
160
161
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
 4992 const Scalar elemMobW = elemVolVars[scvf.insideScvIdx()].mobility(wPhaseIdx); //Mobilities
162
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
 4992 const Scalar elemMobN = elemVolVars[scvf.insideScvIdx()].mobility(nPhaseIdx);
163
164
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
 4992 const Scalar enthW = elemVolVars[scvf.insideScvIdx()].enthalpy(wPhaseIdx); //Enthalpies
165
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
 4992 const Scalar enthN = elemVolVars[scvf.insideScvIdx()].enthalpy(nPhaseIdx);
166
167 1664 const Scalar wellRadius = 0.50 * 0.3048; // 0.50 ft as specified by SPE9
168
169
170 1664 const Scalar gridHeight_ = 0.5;
171 1664 const Scalar effectiveRadius_ = 0.208 * gridHeight_; //Peaceman's Well Model
172
173 using std::log;
174 // divided by molarMass() of water to convert from kg/m s to mol/m s
175 3328 const Scalar qW = (((2*3.1415*0.5*4e-14)/(log(effectiveRadius_/wellRadius))) *
176 1664 densityW * elemMobW * ( elemPressW-pOut_))/0.01801528;
177 // divided by molarMass() of HeavyOil to convert from kg/m s to mol/m s
178 3328 const Scalar qN = (((2*3.1415*0.5*4e-14)/(log(effectiveRadius_/wellRadius))) *
179 1664 densityN * elemMobN * (elemPressN-pOut_))/0.35;
180
181 Scalar qE;
182 // without cooling:
183 // qE = qW*0.018*enthW + qN*enthN*0.350;
184
185 // with cooling: see Diplomarbeit Stefan Roll, Sept. 2015
186
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1664 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1664 times.
 4992 Scalar wT = elemVolVars[scvf.insideScvIdx()].temperature(); // well temperature
187
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1664 times.
 1664 if ( wT > 495. )
188 {
189 ✗ qE = qW*0.018*enthW + qN*enthN*0.350 + (wT-495.)*5000.; // ~3x injected enthalpy
190 ✗ std::cout<< "Cooling now! Extracted enthalpy: " << qE << std::endl;
191 } else {
192 1664 qE = qW*0.018*enthW + qN*enthN*0.350;
193 }
194
195
196 1664 values[contiWEqIdx] = qW;
197 1664 values[contiNEqIdx] = qN;
198 3328 values[energyEqIdx] = qE;
199 }
200 267904 return values;
201 }
202
203 // \}
204
205 /*!
206 * \name Volume terms
207 */
208 // \{
209
210 /*!
211 * \brief Evaluates the initial value for a control volume.
212 *
213 * \param globalPos The position for which the initial condition should be evaluated
214 */
215 PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
216 {
217 5002 return initial_(globalPos);
218 }
219
220 private:
221 // internal method for the initial condition (reused for the
222 // dirichlet conditions!)
223 ✗ PrimaryVariables initial_(const GlobalPosition &globalPos) const
224 {
225 8585 PrimaryVariables values(0.0);
226 6084 values.setState(wnPhaseOnly);
227 8585 Scalar densityW = 1000.0;
228 12168 values[pressureIdx] = 101300.0 + (maxDepth_ - globalPos[1])*densityW*9.81;
229
230 6084 values[switch1Idx] = 295.13; // temperature
231 12168 values[switch2Idx] = 0.3; // NAPL saturation
232 ✗ return values;
233 }
234
235 Scalar maxDepth_;
236 static constexpr Scalar eps_ = 1e-6;
237 Scalar pIn_;
238 Scalar pOut_;
239
240 std::string name_;
241 };
242 } // end namespace Dumux
243
244 #endif
245