blackoillocalresidual.hh
Go to the documentation of this file.
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 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_BLACK_OIL_LOCAL_RESIDUAL_HH
29#define EWOMS_BLACK_OIL_LOCAL_RESIDUAL_HH
30
31#include "blackoilproperties.hh"
41#include <opm/material/fluidstates/BlackOilFluidState.hpp>
42
43namespace Opm {
49template <class TypeTag>
50class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalResidual>
51{
61
62 enum { conti0EqIdx = Indices::conti0EqIdx };
63 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
64 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
65 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
66
67 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
68 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
69 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
70
71 enum { gasCompIdx = FluidSystem::gasCompIdx };
72 enum { oilCompIdx = FluidSystem::oilCompIdx };
73 enum { waterCompIdx = FluidSystem::waterCompIdx };
74 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
75
76 static const bool waterEnabled = Indices::waterEnabled;
77 static const bool gasEnabled = Indices::gasEnabled;
78 static const bool oilEnabled = Indices::oilEnabled;
79 static const bool compositionSwitchEnabled = (compositionSwitchIdx >= 0);
80
81 static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
82 static constexpr bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
83 static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
84 static constexpr bool enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>();
85
86 using Toolbox = MathToolbox<Evaluation>;
96
97public:
101 template <class LhsEval>
102 void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
103 const ElementContext& elemCtx,
104 unsigned dofIdx,
105 unsigned timeIdx) const
106 {
107 // retrieve the intensive quantities for the SCV at the specified point in time
108 const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
109 const auto& fs = intQuants.fluidState();
110
111 storage = 0.0;
112
113 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
114 if (!FluidSystem::phaseIsActive(phaseIdx)) {
115 if (Indices::numPhases == 3) { // add trivial equation for the pseudo phase
116 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
117 if (timeIdx == 0)
118 storage[conti0EqIdx + activeCompIdx] = variable<LhsEval>(0.0, conti0EqIdx + activeCompIdx);
119 else
120 storage[conti0EqIdx + activeCompIdx] = 0.0;
121 }
122 continue;
123 }
124
125 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
126 LhsEval surfaceVolume =
127 Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
128 * Toolbox::template decay<LhsEval>(fs.invB(phaseIdx))
129 * Toolbox::template decay<LhsEval>(intQuants.porosity());
130
131 storage[conti0EqIdx + activeCompIdx] += surfaceVolume;
132
133 // account for dissolved gas
134 if (phaseIdx == oilPhaseIdx && FluidSystem::enableDissolvedGas()) {
135 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
136 storage[conti0EqIdx + activeGasCompIdx] +=
137 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs())
138 * surfaceVolume;
139 }
140
141 // account for dissolved gas in water phase
142 if (phaseIdx == waterPhaseIdx && FluidSystem::enableDissolvedGasInWater()) {
143 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
144 storage[conti0EqIdx + activeGasCompIdx] +=
145 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rsw())
146 * surfaceVolume;
147 }
148
149 // account for vaporized oil
150 if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedOil()) {
151 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
152 storage[conti0EqIdx + activeOilCompIdx] +=
153 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rv())
154 * surfaceVolume;
155 }
156
157 // account for vaporized water
158 if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedWater()) {
159 unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
160 storage[conti0EqIdx + activeWaterCompIdx] +=
161 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rvw())
162 * surfaceVolume;
163 }
164 }
165
166 adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex());
167
168 // deal with solvents (if present)
169 SolventModule::addStorage(storage, intQuants);
170
171 // deal with zFracton (if present)
172 ExtboModule::addStorage(storage, intQuants);
173
174 // deal with polymer (if present)
175 PolymerModule::addStorage(storage, intQuants);
176
177 // deal with energy (if present)
178 EnergyModule::addStorage(storage, intQuants);
179
180 // deal with foam (if present)
181 FoamModule::addStorage(storage, intQuants);
182
183 // deal with salt (if present)
184 BrineModule::addStorage(storage, intQuants);
185
186 // deal with micp (if present)
187 MICPModule::addStorage(storage, intQuants);
188 }
189
193 void computeFlux(RateVector& flux,
194 const ElementContext& elemCtx,
195 unsigned scvfIdx,
196 unsigned timeIdx) const
197 {
198 assert(timeIdx == 0);
199
200 flux = 0.0;
201
202 const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
203 unsigned focusDofIdx = elemCtx.focusDofIndex();
204 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
205 if (!FluidSystem::phaseIsActive(phaseIdx))
206 continue;
207
208 unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
209 const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
210 unsigned pvtRegionIdx = up.pvtRegionIndex();
211 if (upIdx == focusDofIdx)
212 evalPhaseFluxes_<Evaluation>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
213 else
214 evalPhaseFluxes_<Scalar>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
215 }
216 // deal with solvents (if present)
217 SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
218
219 // deal with zFracton (if present)
220 ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
221
222 // deal with polymer (if present)
223 PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
224
225 // deal with energy (if present)
226 EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
227
228 // deal with foam (if present)
229 FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
230
231 // deal with salt (if present)
232 BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
233
234 // deal with micp (if present)
235 MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
236
237 DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
238
239 // deal with convective mixing (if enabled)
240 ConvectiveMixingModule::addConvectiveMixingFlux(flux,elemCtx, scvfIdx, timeIdx);
241 }
242
246 void computeSource(RateVector& source,
247 const ElementContext& elemCtx,
248 unsigned dofIdx,
249 unsigned timeIdx) const
250 {
251 // retrieve the source term intrinsic to the problem
252 elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
253
254 // deal with MICP (if present)
255 MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
256
257 // scale the source term of the energy equation
258 if (enableEnergy)
259 source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
260 }
261
266 template <class UpEval, class FluidState>
267 static void evalPhaseFluxes_(RateVector& flux,
268 unsigned phaseIdx,
269 unsigned pvtRegionIdx,
270 const ExtensiveQuantities& extQuants,
271 const FluidState& upFs)
272 {
273 const auto& invB = getInvB_<FluidSystem, FluidState, UpEval>(upFs, phaseIdx, pvtRegionIdx);
274 const auto& surfaceVolumeFlux = invB*extQuants.volumeFlux(phaseIdx);
275 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
276
277 if (blackoilConserveSurfaceVolume)
278 flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux;
279 else
280 flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux*FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
281
282 if (phaseIdx == oilPhaseIdx) {
283 // dissolved gas (in the oil phase).
284 if (FluidSystem::enableDissolvedGas()) {
285 const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
286
287 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
288 if (blackoilConserveSurfaceVolume)
289 flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux;
290 else
291 flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
292 }
293 } else if (phaseIdx == waterPhaseIdx) {
294 // dissolved gas (in the water phase).
295 if (FluidSystem::enableDissolvedGasInWater()) {
296 const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
297
298 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
299 if (blackoilConserveSurfaceVolume)
300 flux[conti0EqIdx + activeGasCompIdx] += Rsw*surfaceVolumeFlux;
301 else
302 flux[conti0EqIdx + activeGasCompIdx] += Rsw*surfaceVolumeFlux*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
303 }
304 }
305 else if (phaseIdx == gasPhaseIdx) {
306 // vaporized oil (in the gas phase).
307 if (FluidSystem::enableVaporizedOil()) {
308 const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
309
310 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
311 if (blackoilConserveSurfaceVolume)
312 flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux;
313 else
314 flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux*FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
315 }
316 // vaporized water (in the gas phase).
317 if (FluidSystem::enableVaporizedWater()) {
318 const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
319
320 unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
321 if (blackoilConserveSurfaceVolume)
322 flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux;
323 else
324 flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux*FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
325 }
326 }
327 }
328
340 template <class Scalar>
341 static void adaptMassConservationQuantities_(Dune::FieldVector<Scalar, numEq>& container, unsigned pvtRegionIdx)
342 {
343 if (blackoilConserveSurfaceVolume)
344 return;
345
346 // convert "surface volume" to mass. this is complicated a bit by the fact that
347 // not all phases are necessarily enabled. (we here assume that if a fluid phase
348 // is disabled, its respective "main" component is not considered as well.)
349
350 if (waterEnabled) {
351 unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
352 container[conti0EqIdx + activeWaterCompIdx] *=
353 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
354 }
355
356 if (gasEnabled) {
357 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
358 container[conti0EqIdx + activeGasCompIdx] *=
359 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
360 }
361
362 if (oilEnabled) {
363 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
364 container[conti0EqIdx + activeOilCompIdx] *=
365 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
366 }
367 }
368};
369
370} // namespace Opm
371
372#endif
Contains the classes required to extend the black-oil model by brine.
Classes required for dynamic convective mixing.
Classes required for molecular diffusion.
Contains the classes required to extend the black-oil model by energy.
Contains the classes required to extend the black-oil model by solvent component. For details,...
Contains the classes required to extend the black-oil model to include the effects of foam.
Contains the classes required to extend the black-oil model by MICP.
Contains the classes required to extend the black-oil model by polymer.
Declares the properties required by the black oil model.
Contains the classes required to extend the black-oil model by solvents.
Contains the high level supplements required to extend the black oil model by brine.
Definition: blackoilbrinemodules.hh:58
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilbrinemodules.hh:260
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilbrinemodules.hh:228
Definition: blackoilconvectivemixingmodule.hh:58
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: blackoildiffusionmodule.hh:48
Contains the high level supplements required to extend the black oil model by energy.
Definition: blackoilenergymodules.hh:52
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilenergymodules.hh:168
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilenergymodules.hh:141
Contains the high level supplements required to extend the black oil model.
Definition: blackoilextbomodules.hh:69
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilextbomodules.hh:320
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilextbomodules.hh:354
Contains the high level supplements required to extend the black oil model to include the effects of ...
Definition: blackoilfoammodules.hh:59
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilfoammodules.hh:247
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:287
Calculates the local residual of the black oil model.
Definition: blackoillocalresidual.hh:51
void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx) const
Evaluates the total mass flux of all conservation quantities over a face of a sub-control volume.
Definition: blackoillocalresidual.hh:193
void computeSource(RateVector &source, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx) const
Calculate the source term of the equation.
Definition: blackoillocalresidual.hh:246
static void evalPhaseFluxes_(RateVector &flux, unsigned phaseIdx, unsigned pvtRegionIdx, const ExtensiveQuantities &extQuants, const FluidState &upFs)
Helper function to calculate the flux of mass in terms of conservation quantities via specific fluid ...
Definition: blackoillocalresidual.hh:267
static void adaptMassConservationQuantities_(Dune::FieldVector< Scalar, numEq > &container, unsigned pvtRegionIdx)
Helper function to convert the mass-related parts of a Dune::FieldVector that stores conservation qua...
Definition: blackoillocalresidual.hh:341
void computeStorage(Dune::FieldVector< LhsEval, numEq > &storage, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx) const
Evaluate the amount all conservation quantities (e.g. phase mass) within a finite sub-control volume.
Definition: blackoillocalresidual.hh:102
Contains the high level supplements required to extend the black oil model by MICP.
Definition: blackoilmicpmodules.hh:56
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilmicpmodules.hh:271
static void addSource(RateVector &source, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: blackoilmicpmodules.hh:300
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilmicpmodules.hh:240
Contains the high level supplements required to extend the black oil model by polymer.
Definition: blackoilpolymermodules.hh:61
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilpolymermodules.hh:473
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilpolymermodules.hh:512
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilsolventmodules.hh:69
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilsolventmodules.hh:431
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilsolventmodules.hh:463
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235