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
35
36#include <algorithm>
37#include <type_traits>
38
39namespace Opm {
40
46template <class TypeTag>
47class BlackOilBoundaryRateVector : public GetPropType<TypeTag, Properties::RateVector>
48{
57
58 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
59 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
60 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
61 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
62 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
63 enum { contiEnergyEqIdx = Indices::contiEnergyEqIdx };
64 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
65 enum { enableMICP = Indices::enableMICP };
66
67 static constexpr bool blackoilConserveSurfaceVolume =
68 getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
69
71
72public:
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 const unsigned focusDofIdx = context.focusDofIndex();
103 const 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 }
127 else if (pBoundary > pInside) {
128 using RhsEval = std::conditional_t<std::is_same_v<typename FluidState::Scalar, Evaluation>,
129 Evaluation, Scalar>;
130 // influx
131 LocalResidual::template evalPhaseFluxes_<RhsEval>(tmp,
132 phaseIdx,
133 insideIntQuants.pvtRegionIndex(),
134 extQuants,
135 fluidState);
136 }
137
138 for (unsigned i = 0; i < tmp.size(); ++i) {
139 (*this)[i] += tmp[i];
140 }
141
142 // energy conservation
143 if constexpr (enableEnergy) {
144 Evaluation density;
145 Evaluation specificEnthalpy;
146 if (pBoundary > pInside) {
147 if (focusDofIdx == interiorDofIdx) {
148 density = fluidState.density(phaseIdx);
149 specificEnthalpy = fluidState.enthalpy(phaseIdx);
150 }
151 else {
152 density = getValue(fluidState.density(phaseIdx));
153 specificEnthalpy = getValue(fluidState.enthalpy(phaseIdx));
154 }
155 }
156 else if (focusDofIdx == interiorDofIdx) {
157 density = insideIntQuants.fluidState().density(phaseIdx);
158 specificEnthalpy = insideIntQuants.fluidState().enthalpy(phaseIdx);
159 }
160 else {
161 density = getValue(insideIntQuants.fluidState().density(phaseIdx));
162 specificEnthalpy = getValue(insideIntQuants.fluidState().enthalpy(phaseIdx));
163 }
164
165 const Evaluation enthalpyRate = density * extQuants.volumeFlux(phaseIdx) * specificEnthalpy;
166 EnergyModule::addToEnthalpyRate(*this, enthalpyRate *
167 getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>());
168 }
169 }
170
171 if constexpr (enableSolvent) {
172 (*this)[Indices::contiSolventEqIdx] = extQuants.solventVolumeFlux();
173 if (blackoilConserveSurfaceVolume) {
174 (*this)[Indices::contiSolventEqIdx] *= insideIntQuants.solventInverseFormationVolumeFactor();
175 }
176 else {
177 (*this)[Indices::contiSolventEqIdx] *= insideIntQuants.solventDensity();
178 }
179 }
180
181 if constexpr (enablePolymer) {
182 (*this)[Indices::contiPolymerEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) *
183 insideIntQuants.polymerConcentration();
184 }
185
186 if constexpr (enableMICP) {
187 (*this)[Indices::contiMicrobialEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) *
188 insideIntQuants.microbialConcentration();
189 (*this)[Indices::contiOxygenEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) *
190 insideIntQuants.oxygenConcentration();
191 (*this)[Indices::contiUreaEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) *
192 insideIntQuants.ureaConcentration();
193 // since the urea concentration can be much larger than 1, then we apply a scaling factor
194 (*this)[Indices::contiUreaEqIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
195 }
196
197 // make sure that the right mass conservation quantities are used
198 LocalResidual::adaptMassConservationQuantities_(*this, insideIntQuants.pvtRegionIndex());
199
200 // heat conduction
201 if constexpr (enableEnergy) {
202 EnergyModule::addToEnthalpyRate(*this, extQuants.energyFlux() *
203 getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>());
204 }
205
206#ifndef NDEBUG
207 for (unsigned i = 0; i < numEq; ++i) {
208 Valgrind::CheckDefined((*this)[i]);
209 }
210 Valgrind::CheckDefined(*this);
211#endif
212 }
213
217 template <class Context, class FluidState>
218 void setInFlow(const Context& context,
219 unsigned bfIdx,
220 unsigned timeIdx,
221 const FluidState& fluidState)
222 {
223 this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
224
225 // we only allow fluxes in the direction opposite to the outer
226 // unit normal
227 std::for_each(this->begin(), this->end(),
228 [](auto& val) { val = std::min(Scalar(0), val); });
229 }
230
234 template <class Context, class FluidState>
235 void setOutFlow(const Context& context,
236 unsigned bfIdx,
237 unsigned timeIdx,
238 const FluidState& fluidState)
239 {
240 this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
241
242 // we only allow fluxes in the same direction as the outer
243 // unit normal
244 std::for_each(this->begin(), this->end(),
245 [](auto& val) { val = std::max(Scalar(0), val); });
246 }
247
252 { (*this) = Scalar(0); }
253
261 template <class Context, class FluidState>
262 void setThermalFlow([[maybe_unused]] const Context& context,
263 [[maybe_unused]] unsigned bfIdx,
264 [[maybe_unused]] unsigned timeIdx,
265 [[maybe_unused]] const FluidState& boundaryFluidState)
266 {
267 // set the mass no-flow condition
268 setNoFlow();
269
270 // if we do not conserve energy there is nothing we should do in addition
271 if constexpr (enableEnergy) {
272 ExtensiveQuantities extQuants;
273 extQuants.updateBoundary(context, bfIdx, timeIdx, boundaryFluidState);
274
275 (*this)[contiEnergyEqIdx] += extQuants.energyFlux();
276
277#ifndef NDEBUG
278 for (unsigned i = 0; i < numEq; ++i) {
279 Valgrind::CheckDefined((*this)[i]);
280 }
281 Valgrind::CheckDefined(*this);
282#endif
283 }
284 }
285};
286
287} // namespace Opm
288
289#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:48
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:235
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
BlackOilBoundaryRateVector()=default
Default constructor.
void setNoFlow()
Specify a no-flow boundary for all conserved quantities.
Definition: blackoilboundaryratevector.hh:251
void setInFlow(const Context &context, unsigned bfIdx, unsigned timeIdx, const FluidState &fluidState)
Specify an inflow boundary.
Definition: blackoilboundaryratevector.hh:218
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:262
Contains the high level supplements required to extend the black oil model by energy.
Definition: blackoilenergymodules.hh:58
static void addToEnthalpyRate(RateVector &flux, const Evaluation &hRate)
Definition: blackoilenergymodules.hh:251
Definition: blackoilbioeffectsmodules.hh:43
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:233