blackoilboundaryratevector.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_BOUNDARY_RATE_VECTOR_HH
29#define EWOMS_BLACK_OIL_BOUNDARY_RATE_VECTOR_HH
30
31#include <opm/material/common/Valgrind.hpp>
32#include <opm/material/constraintsolvers/NcpFlash.hpp>
33
36
37namespace Opm {
38
44template <class TypeTag>
45class BlackOilBoundaryRateVector : public GetPropType<TypeTag, Properties::RateVector>
46{
55
56 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
57 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
58 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
59 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
60 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
61 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
62 enum { conti0EqIdx = Indices::conti0EqIdx };
63 enum { contiEnergyEqIdx = Indices::contiEnergyEqIdx };
64 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
65 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
66
67 static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
68
70
71public:
76 {}
77
81 BlackOilBoundaryRateVector(Scalar value) : ParentType(value)
82 {}
83
89
93 template <class Context, class FluidState>
94 void setFreeFlow(const Context& context,
95 unsigned bfIdx,
96 unsigned timeIdx,
97 const FluidState& fluidState)
98 {
99 ExtensiveQuantities extQuants;
100 extQuants.updateBoundary(context, bfIdx, timeIdx, fluidState);
101 const auto& insideIntQuants = context.intensiveQuantities(bfIdx, timeIdx);
102 unsigned focusDofIdx = context.focusDofIndex();
103 unsigned interiorDofIdx = context.interiorScvIndex(bfIdx, timeIdx);
104
106 // advective fluxes of all components in all phases
108 (*this) = 0.0;
109 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
110 if (!FluidSystem::phaseIsActive(phaseIdx)) {
111 continue;
112 }
113 const auto& pBoundary = fluidState.pressure(phaseIdx);
114 const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx);
115
116 RateVector tmp;
117
118 // mass conservation
119 if (pBoundary < pInside)
120 // outflux
121 LocalResidual::template evalPhaseFluxes_<Evaluation>(tmp,
122 phaseIdx,
123 insideIntQuants.pvtRegionIndex(),
124 extQuants,
125 insideIntQuants.fluidState());
126 else if (pBoundary > pInside) {
127 using RhsEval = typename std::conditional<std::is_same<typename FluidState::Scalar, Evaluation>::value,
128 Evaluation, Scalar>::type;
129 // influx
130 LocalResidual::template evalPhaseFluxes_<RhsEval>(tmp,
131 phaseIdx,
132 insideIntQuants.pvtRegionIndex(),
133 extQuants,
134 fluidState);
135 }
136
137 for (unsigned i = 0; i < tmp.size(); ++i)
138 (*this)[i] += tmp[i];
139
140 // energy conservation
141 if constexpr (enableEnergy) {
142 Evaluation density;
143 Evaluation specificEnthalpy;
144 if (pBoundary > pInside) {
145 if (focusDofIdx == interiorDofIdx) {
146 density = fluidState.density(phaseIdx);
147 specificEnthalpy = fluidState.enthalpy(phaseIdx);
148 }
149 else {
150 density = getValue(fluidState.density(phaseIdx));
151 specificEnthalpy = getValue(fluidState.enthalpy(phaseIdx));
152 }
153 }
154 else if (focusDofIdx == interiorDofIdx) {
155 density = insideIntQuants.fluidState().density(phaseIdx);
156 specificEnthalpy = insideIntQuants.fluidState().enthalpy(phaseIdx);
157 }
158 else {
159 density = getValue(insideIntQuants.fluidState().density(phaseIdx));
160 specificEnthalpy = getValue(insideIntQuants.fluidState().enthalpy(phaseIdx));
161 }
162
163 Evaluation enthalpyRate = density*extQuants.volumeFlux(phaseIdx)*specificEnthalpy;
164 EnergyModule::addToEnthalpyRate(*this, enthalpyRate*getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>());
165 }
166 }
167
168 if constexpr (enableSolvent) {
169 (*this)[Indices::contiSolventEqIdx] = extQuants.solventVolumeFlux();
170 if (blackoilConserveSurfaceVolume)
171 (*this)[Indices::contiSolventEqIdx] *= insideIntQuants.solventInverseFormationVolumeFactor();
172 else
173 (*this)[Indices::contiSolventEqIdx] *= insideIntQuants.solventDensity();
174
175 }
176
177 if constexpr (enablePolymer) {
178 (*this)[Indices::contiPolymerEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.polymerConcentration();
179 }
180
181 if constexpr (enableMICP) {
182 (*this)[Indices::contiMicrobialEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.microbialConcentration();
183 (*this)[Indices::contiOxygenEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.oxygenConcentration();
184 (*this)[Indices::contiUreaEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.ureaConcentration();
185 }
186
187 // make sure that the right mass conservation quantities are used
188 LocalResidual::adaptMassConservationQuantities_(*this, insideIntQuants.pvtRegionIndex());
189
190 // heat conduction
191 if constexpr (enableEnergy)
192 EnergyModule::addToEnthalpyRate(*this, extQuants.energyFlux()*getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>());
193
194#ifndef NDEBUG
195 for (unsigned i = 0; i < numEq; ++i) {
196 Valgrind::CheckDefined((*this)[i]);
197 }
198 Valgrind::CheckDefined(*this);
199#endif
200 }
201
205 template <class Context, class FluidState>
206 void setInFlow(const Context& context,
207 unsigned bfIdx,
208 unsigned timeIdx,
209 const FluidState& fluidState)
210 {
211 this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
212
213 // we only allow fluxes in the direction opposite to the outer
214 // unit normal
215 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
216 Scalar& val = this->operator[](eqIdx);
217 val = std::min<Scalar>(0.0, val);
218 }
219 }
220
224 template <class Context, class FluidState>
225 void setOutFlow(const Context& context,
226 unsigned bfIdx,
227 unsigned timeIdx,
228 const FluidState& fluidState)
229 {
230 this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
231
232 // we only allow fluxes in the same direction as the outer
233 // unit normal
234 for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
235 Scalar& val = this->operator[](eqIdx);
236 val = std::max( Scalar(0), val);
237 }
238 }
239
244 { (*this) = Scalar(0); }
245
253 template <class Context, class FluidState>
254 void setThermalFlow([[maybe_unused]] const Context& context,
255 [[maybe_unused]] unsigned bfIdx,
256 [[maybe_unused]] unsigned timeIdx,
257 [[maybe_unused]] const FluidState& boundaryFluidState)
258 {
259 // set the mass no-flow condition
260 setNoFlow();
261
262 // if we do not conserve energy there is nothing we should do in addition
263 if constexpr (enableEnergy) {
264 ExtensiveQuantities extQuants;
265 extQuants.updateBoundary(context, bfIdx, timeIdx, boundaryFluidState);
266
267 (*this)[contiEnergyEqIdx] += extQuants.energyFlux();
268
269#ifndef NDEBUG
270 for (unsigned i = 0; i < numEq; ++i)
271 Valgrind::CheckDefined((*this)[i]);
272 Valgrind::CheckDefined(*this);
273#endif
274 }
275 }
276};
277
278} // namespace Opm
279
280#endif
Contains the classes required to extend the black-oil model by energy.
Implements a boundary vector for the fully implicit black-oil model.
Definition: blackoilboundaryratevector.hh:46
BlackOilBoundaryRateVector()
Default constructor.
Definition: blackoilboundaryratevector.hh:75
BlackOilBoundaryRateVector(Scalar value)
Definition: blackoilboundaryratevector.hh:81
void setOutFlow(const Context &context, unsigned bfIdx, unsigned timeIdx, const FluidState &fluidState)
Specify an outflow boundary.
Definition: blackoilboundaryratevector.hh:225
BlackOilBoundaryRateVector & operator=(const BlackOilBoundaryRateVector &value)=default
BlackOilBoundaryRateVector(const BlackOilBoundaryRateVector &value)=default
Copy constructor.
void setFreeFlow(const Context &context, unsigned bfIdx, unsigned timeIdx, const FluidState &fluidState)
Specify a free-flow boundary.
Definition: blackoilboundaryratevector.hh:94
void setNoFlow()
Specify a no-flow boundary for all conserved quantities.
Definition: blackoilboundaryratevector.hh:243
void setInFlow(const Context &context, unsigned bfIdx, unsigned timeIdx, const FluidState &fluidState)
Specify an inflow boundary.
Definition: blackoilboundaryratevector.hh:206
void setThermalFlow(const Context &context, unsigned bfIdx, unsigned timeIdx, const FluidState &boundaryFluidState)
an energy flux that corresponds to the thermal conduction from
Definition: blackoilboundaryratevector.hh:254
Contains the high level supplements required to extend the black oil model by energy.
Definition: blackoilenergymodules.hh:52
static void addToEnthalpyRate(RateVector &flux, const Evaluation &hRate)
Definition: blackoilenergymodules.hh:237
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