73 using Toolbox = MathToolbox<Evaluation>;
75 using TabulatedFunction =
typename BlackOilExtboParams<Scalar>::TabulatedFunction;
76 using Tabulated2DFunction =
typename BlackOilExtboParams<Scalar>::Tabulated2DFunction;
78 static constexpr unsigned zFractionIdx = Indices::zFractionIdx;
79 static constexpr unsigned contiZfracEqIdx = Indices::contiZfracEqIdx;
80 static constexpr unsigned enableExtbo = enableExtboV;
82 static constexpr unsigned numPhases = FluidSystem::numPhases;
83 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
84 static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
85 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
110 static bool primaryVarApplies(
unsigned pvIdx)
112 if constexpr (enableExtbo)
113 return pvIdx == zFractionIdx;
118 static std::string primaryVarName([[maybe_unused]]
unsigned pvIdx)
120 assert(primaryVarApplies(pvIdx));
125 static Scalar primaryVarWeight([[maybe_unused]]
unsigned pvIdx)
127 assert(primaryVarApplies(pvIdx));
130 return static_cast<Scalar
>(1.0);
133 static bool eqApplies(
unsigned eqIdx)
135 if constexpr (enableExtbo)
136 return eqIdx == contiZfracEqIdx;
141 static std::string eqName([[maybe_unused]]
unsigned eqIdx)
143 assert(eqApplies(eqIdx));
145 return "conti^solvent";
148 static Scalar eqWeight([[maybe_unused]]
unsigned eqIdx)
150 assert(eqApplies(eqIdx));
153 return static_cast<Scalar
>(1.0);
156 template <
class LhsEval>
157 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
158 const IntensiveQuantities& intQuants)
160 if constexpr (enableExtbo) {
161 if constexpr (blackoilConserveSurfaceVolume) {
162 storage[contiZfracEqIdx] =
163 Toolbox::template decay<LhsEval>(intQuants.porosity())
164 * Toolbox::template decay<LhsEval>(intQuants.yVolume())
165 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(gasPhaseIdx))
166 * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(gasPhaseIdx));
167 if (FluidSystem::enableDissolvedGas()) {
168 storage[contiZfracEqIdx] +=
169 Toolbox::template decay<LhsEval>(intQuants.porosity())
170 * Toolbox::template decay<LhsEval>(intQuants.xVolume())
171 * Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs())
172 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(oilPhaseIdx))
173 * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(oilPhaseIdx));
178 const Scalar regWghtFactor = 1.0e-6;
179 storage[contiZfracEqIdx] += regWghtFactor*(1.0-Toolbox::template decay<LhsEval>(intQuants.zFraction()))
180 + regWghtFactor*Toolbox::template decay<LhsEval>(intQuants.porosity())
181 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(gasPhaseIdx))
182 * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(gasPhaseIdx));
183 storage[contiZfracEqIdx-1] += regWghtFactor*Toolbox::template decay<LhsEval>(intQuants.zFraction());
186 throw std::runtime_error(
"Only component conservation in terms of surface volumes is implemented. ");
191 static void computeFlux([[maybe_unused]] RateVector& flux,
192 [[maybe_unused]]
const ElementContext& elemCtx,
193 [[maybe_unused]]
unsigned scvfIdx,
194 [[maybe_unused]]
unsigned timeIdx)
197 if constexpr (enableExtbo) {
198 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
200 if constexpr (blackoilConserveSurfaceVolume) {
201 unsigned inIdx = extQuants.interiorIndex();
203 unsigned upIdxGas =
static_cast<unsigned>(extQuants.upstreamIndex(gasPhaseIdx));
204 const auto& upGas = elemCtx.intensiveQuantities(upIdxGas, timeIdx);
205 const auto& fsGas = upGas.fluidState();
206 if (upIdxGas == inIdx) {
207 flux[contiZfracEqIdx] =
208 extQuants.volumeFlux(gasPhaseIdx)
210 * fsGas.invB(gasPhaseIdx);
213 flux[contiZfracEqIdx] =
214 extQuants.volumeFlux(gasPhaseIdx)
215 * (decay<Scalar>(upGas.yVolume()))
216 * decay<Scalar>(fsGas.invB(gasPhaseIdx));
218 if (FluidSystem::enableDissolvedGas()) {
219 unsigned upIdxOil =
static_cast<unsigned>(extQuants.upstreamIndex(oilPhaseIdx));
220 const auto& upOil = elemCtx.intensiveQuantities(upIdxOil, timeIdx);
221 const auto& fsOil = upOil.fluidState();
222 if (upIdxOil == inIdx) {
223 flux[contiZfracEqIdx] +=
224 extQuants.volumeFlux(oilPhaseIdx)
227 * fsOil.invB(oilPhaseIdx);
230 flux[contiZfracEqIdx] +=
231 extQuants.volumeFlux(oilPhaseIdx)
232 * decay<Scalar>(upOil.xVolume())
233 * decay<Scalar>(fsOil.Rs())
234 * decay<Scalar>(fsOil.invB(oilPhaseIdx));
239 throw std::runtime_error(
"Only component conservation in terms of surface volumes is implemented. ");
250 if constexpr (enableExtbo)
251 priVars[zFractionIdx] = zFraction;
258 const PrimaryVariables& oldPv,
259 const EqVector& delta)
261 if constexpr (enableExtbo)
263 newPv[zFractionIdx] = oldPv[zFractionIdx] - delta[zFractionIdx];
275 return static_cast<Scalar
>(0.0);
284 return std::abs(Toolbox::scalarValue(resid[contiZfracEqIdx]));
287 template <
class DofEntity>
288 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
290 if constexpr (enableExtbo) {
291 unsigned dofIdx = model.dofMapper().index(dof);
293 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
294 outstream << priVars[zFractionIdx];
298 template <
class DofEntity>
299 static void deserializeEntity(Model& model, std::istream& instream,
const DofEntity& dof)
301 if constexpr (enableExtbo) {
302 unsigned dofIdx = model.dofMapper().index(dof);
304 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
305 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
307 instream >> priVars0[zFractionIdx];
310 priVars1 = priVars0[zFractionIdx];
314 template <
typename Value>
315 static Value xVolume(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
316 return params_.X_[pvtRegionIdx].eval(z, pressure,
true);
319 template <
typename Value>
320 static Value yVolume(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
321 return params_.Y_[pvtRegionIdx].eval(z, pressure,
true);
324 template <
typename Value>
325 static Value pbubRs(
unsigned pvtRegionIdx,
const Value& z,
const Value& rs) {
326 return params_.PBUB_RS_[pvtRegionIdx].eval(z, rs,
true);
329 template <
typename Value>
330 static Value pbubRv(
unsigned pvtRegionIdx,
const Value& z,
const Value& rv) {
331 return params_.PBUB_RV_[pvtRegionIdx].eval(z, rv,
true);
334 template <
typename Value>
335 static Value oilViscosity(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
336 return params_.VISCO_[pvtRegionIdx].eval(z, pressure,
true);
339 template <
typename Value>
340 static Value gasViscosity(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
341 return params_.VISCG_[pvtRegionIdx].eval(z, pressure,
true);
344 template <
typename Value>
345 static Value bo(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
346 return params_.BO_[pvtRegionIdx].eval(z, pressure,
true);
349 template <
typename Value>
350 static Value bg(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
351 return params_.BG_[pvtRegionIdx].eval(z, pressure,
true);
354 template <
typename Value>
355 static Value rs(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
356 return params_.RS_[pvtRegionIdx].eval(z, pressure,
true);
359 template <
typename Value>
360 static Value rv(
unsigned pvtRegionIdx,
const Value& pressure,
const Value& z) {
361 return params_.RV_[pvtRegionIdx].eval(z, pressure,
true);
364 static Scalar referenceDensity(
unsigned regionIdx) {
365 return params_.zReferenceDensity_[regionIdx];
368 static Scalar zLim(
unsigned regionIdx) {
369 return params_.zLim_[regionIdx];
372 template <
typename Value>
373 static Value oilCmp(
unsigned pvtRegionIdx,
const Value& z) {
374 return params_.oilCmp_[pvtRegionIdx].eval(z,
true);
377 template <
typename Value>
378 static Value gasCmp(
unsigned pvtRegionIdx,
const Value& z) {
379 return params_.gasCmp_[pvtRegionIdx].eval(z,
true);
383 static BlackOilExtboParams<Scalar> params_;
430 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
431 unsigned pvtRegionIdx = priVars.pvtRegionIndex();
432 auto& fs = asImp_().fluidState_;
434 zFraction_ = priVars.makeEvaluation(zFractionIdx, timeIdx);
436 oilViscosity_ = ExtboModule::oilViscosity(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
437 gasViscosity_ = ExtboModule::gasViscosity(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
439 bo_ = ExtboModule::bo(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
440 bg_ = ExtboModule::bg(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
442 bz_ = ExtboModule::bg(pvtRegionIdx, fs.pressure(oilPhaseIdx), Evaluation{0.99});
444 if (FluidSystem::enableDissolvedGas())
445 rs_ = ExtboModule::rs(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
449 if (FluidSystem::enableVaporizedOil())
450 rv_ = ExtboModule::rv(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
454 xVolume_ = ExtboModule::xVolume(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
455 yVolume_ = ExtboModule::yVolume(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
457 Evaluation pbub = fs.pressure(oilPhaseIdx);
459 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
460 static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
461 Scalar sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx).value();
462 if (sw >= thresholdWaterFilledCell)
466 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
467 rs_ = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
468 const Evaluation zLim = ExtboModule::zLim(pvtRegionIdx);
469 if (zFraction_ > zLim) {
470 pbub = ExtboModule::pbubRs(pvtRegionIdx, zLim, rs_);
472 pbub = ExtboModule::pbubRs(pvtRegionIdx, zFraction_, rs_);
474 bo_ = ExtboModule::bo(pvtRegionIdx, pbub, zFraction_) + ExtboModule::oilCmp(pvtRegionIdx, zFraction_)*(fs.pressure(oilPhaseIdx)-pbub);
476 xVolume_ = ExtboModule::xVolume(pvtRegionIdx, pbub, zFraction_);
479 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
480 rv_ = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
481 Evaluation rvsat = ExtboModule::rv(pvtRegionIdx, pbub, zFraction_);
482 bg_ = ExtboModule::bg(pvtRegionIdx, pbub, zFraction_) + ExtboModule::gasCmp(pvtRegionIdx, zFraction_)*(rv_-rvsat);
484 yVolume_ = ExtboModule::yVolume(pvtRegionIdx, pbub, zFraction_);