blackoillocalresidualtpfa.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_TPFA_RESIDUAL_HH
29#define EWOMS_BLACK_OIL_LOCAL_TPFA_RESIDUAL_HH
30
31#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
32#include <opm/input/eclipse/Schedule/BCProp.hpp>
33
34#include <opm/material/common/MathToolbox.hpp>
35#include <opm/material/fluidstates/BlackOilFluidState.hpp>
36#include <opm/material/common/ConditionalStorage.hpp>
37
49
50#include <array>
51#include <cassert>
52#include <stdexcept>
53#include <string>
54
55namespace Opm {
61template <class TypeTag>
62class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLocalResidual>
63{
74 using FluidState = typename IntensiveQuantities::FluidState;
75
76 enum { conti0EqIdx = Indices::conti0EqIdx };
77 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
78 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
79 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
80
81 enum { dimWorld = GridView::dimensionworld };
82 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
83 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
84 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
85
86 enum { gasCompIdx = FluidSystem::gasCompIdx };
87 enum { oilCompIdx = FluidSystem::oilCompIdx };
88 enum { waterCompIdx = FluidSystem::waterCompIdx };
89 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
90
91 static constexpr bool waterEnabled = Indices::waterEnabled;
92 static constexpr bool gasEnabled = Indices::gasEnabled;
93 static constexpr bool oilEnabled = Indices::oilEnabled;
94 static constexpr bool compositionSwitchEnabled = compositionSwitchIdx >= 0;
95
96 static constexpr bool blackoilConserveSurfaceVolume =
97 getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
98
99 static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
100 static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
101 static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
102 static constexpr bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
103 static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
104 static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
105 static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
106 static constexpr bool enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>();
107 static constexpr bool enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>();
108 static constexpr bool enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();
109 static constexpr bool enableMICP = Indices::enableMICP;
110
119 using ConvectiveMixingModuleParam = typename ConvectiveMixingModule::ConvectiveMixingModuleParam;
120
123
124 using Toolbox = MathToolbox<Evaluation>;
125
126public:
128 {
129 double trans;
130 double faceArea;
131 double thpres;
132 double dZg;
133 FaceDir::DirEnum faceDir;
134 double Vin;
135 double Vex;
136 ConditionalStorage<enableEnergy, double> inAlpha;
137 ConditionalStorage<enableEnergy, double> outAlpha;
138 ConditionalStorage<enableDiffusion, double> diffusivity;
139 ConditionalStorage<enableDispersion, double> dispersivity;
140 };
141
143 {
144 ConvectiveMixingModuleParam convectiveMixingModuleParam;
145 };
146
150 template <class LhsEval>
151 void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
152 const ElementContext& elemCtx,
153 unsigned dofIdx,
154 unsigned timeIdx) const
155 {
156 const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
157 computeStorage(storage, intQuants);
158 }
159
160 template <class LhsEval>
161 static void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
162 const IntensiveQuantities& intQuants)
163 {
164 OPM_TIMEBLOCK_LOCAL(computeStorage, Subsystem::Assembly);
165 // retrieve the intensive quantities for the SCV at the specified point in time
166 const auto& fs = intQuants.fluidState();
167 storage = 0.0;
168
169 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
170 if (!FluidSystem::phaseIsActive(phaseIdx)) {
171 continue;
172 }
173 unsigned activeCompIdx =
174 FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
175 LhsEval surfaceVolume =
176 Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx)) *
177 Toolbox::template decay<LhsEval>(fs.invB(phaseIdx)) *
178 Toolbox::template decay<LhsEval>(intQuants.porosity());
179
180 storage[conti0EqIdx + activeCompIdx] += surfaceVolume;
181
182 // account for dissolved gas
183 if (phaseIdx == oilPhaseIdx && FluidSystem::enableDissolvedGas()) {
184 unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
185 storage[conti0EqIdx + activeGasCompIdx] +=
186 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs()) *
187 surfaceVolume;
188 }
189
190 // account for dissolved gas in water
191 if (phaseIdx == waterPhaseIdx && FluidSystem::enableDissolvedGasInWater()) {
192 unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
193 storage[conti0EqIdx + activeGasCompIdx] +=
194 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rsw()) *
195 surfaceVolume;
196 }
197
198 // account for vaporized oil
199 if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedOil()) {
200 unsigned activeOilCompIdx = FluidSystem::canonicalToActiveCompIdx(oilCompIdx);
201 storage[conti0EqIdx + activeOilCompIdx] +=
202 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rv()) *
203 surfaceVolume;
204 }
205
206 // account for vaporized water
207 if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedWater()) {
208 unsigned activeWaterCompIdx = FluidSystem::canonicalToActiveCompIdx(waterCompIdx);
209 storage[conti0EqIdx + activeWaterCompIdx] +=
210 Toolbox::template decay<LhsEval>(intQuants.fluidState().Rvw()) *
211 surfaceVolume;
212 }
213 }
214
215 adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex());
216
217 // deal with solvents (if present)
218 SolventModule::addStorage(storage, intQuants);
219
220 // deal with zFracton (if present)
221 ExtboModule::addStorage(storage, intQuants);
222
223 // deal with polymer (if present)
224 PolymerModule::addStorage(storage, intQuants);
225
226 // deal with energy (if present)
227 EnergyModule::addStorage(storage, intQuants);
228
229 // deal with foam (if present)
230 FoamModule::addStorage(storage, intQuants);
231
232 // deal with salt (if present)
233 BrineModule::addStorage(storage, intQuants);
234
235 // deal with bioeffects (if present)
236 BioeffectsModule::addStorage(storage, intQuants);
237 }
238
244 static void computeFlux(RateVector& flux,
245 RateVector& darcy,
246 const unsigned globalIndexIn,
247 const unsigned globalIndexEx,
248 const IntensiveQuantities& intQuantsIn,
249 const IntensiveQuantities& intQuantsEx,
250 const ResidualNBInfo& nbInfo,
251 const ModuleParams& moduleParams)
252 {
253 OPM_TIMEBLOCK_LOCAL(computeFlux, Subsystem::Assembly);
254 flux = 0.0;
255 darcy = 0.0;
256
257 calculateFluxes_(flux,
258 darcy,
259 intQuantsIn,
260 intQuantsEx,
261 globalIndexIn,
262 globalIndexEx,
263 nbInfo,
264 moduleParams);
265 }
266
267 // This function demonstrates compatibility with the ElementContext-based interface.
268 // Actually using it will lead to double work since the element context already contains
269 // fluxes through its stored ExtensiveQuantities.
270 static void computeFlux(RateVector& flux,
271 const ElementContext& elemCtx,
272 unsigned scvfIdx,
273 unsigned timeIdx)
274 {
275 OPM_TIMEBLOCK_LOCAL(computeFlux, Subsystem::Assembly);
276 assert(timeIdx == 0);
277
278 flux = 0.0;
279 RateVector darcy = 0.0;
280 // need for dary flux calculation
281 const auto& problem = elemCtx.problem();
282 const auto& stencil = elemCtx.stencil(timeIdx);
283 const auto& scvf = stencil.interiorFace(scvfIdx);
284
285 unsigned interiorDofIdx = scvf.interiorIndex();
286 unsigned exteriorDofIdx = scvf.exteriorIndex();
287 assert(interiorDofIdx != exteriorDofIdx);
288
289 // unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
290 // unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
291 Scalar Vin = elemCtx.dofVolume(interiorDofIdx, /*timeIdx=*/0);
292 Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, /*timeIdx=*/0);
293 const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx);
294 const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx);
295 Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
296 Scalar faceArea = scvf.area();
297 const auto faceDir = faceDirFromDirId(scvf.dirId());
298 Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx);
299
300 // estimate the gravity correction: for performance reasons we use a simplified
301 // approach for this flux module that assumes that gravity is constant and always
302 // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
303 const Scalar g = problem.gravity()[dimWorld - 1];
304 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
305 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
306
307 // this is quite hacky because the dune grid interface does not provide a
308 // cellCenterDepth() method (so we ask the problem to provide it). The "good"
309 // solution would be to take the Z coordinate of the element centroids, but since
310 // ECL seems to like to be inconsistent on that front, it needs to be done like
311 // here...
312 const Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
313 const Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
314 // the distances from the DOF's depths. (i.e., the additional depth of the
315 // exterior DOF)
316 const Scalar distZ = zIn - zEx;
317 // for thermal harmonic mean of half trans
318 const Scalar inAlpha = problem.thermalHalfTransmissibility(globalIndexIn, globalIndexEx);
319 const Scalar outAlpha = problem.thermalHalfTransmissibility(globalIndexEx, globalIndexIn);
320 const Scalar diffusivity = problem.diffusivity(globalIndexEx, globalIndexIn);
321 const Scalar dispersivity = problem.dispersivity(globalIndexEx, globalIndexIn);
322
323 const ResidualNBInfo res_nbinfo {
324 trans, faceArea, thpres, distZ * g, faceDir, Vin, Vex,
325 inAlpha, outAlpha, diffusivity, dispersivity
326 };
327
328 calculateFluxes_(flux,
329 darcy,
330 intQuantsIn,
331 intQuantsEx,
332 globalIndexIn,
333 globalIndexEx,
334 res_nbinfo,
335 problem.moduleParams());
336 }
337
338 static void calculateFluxes_(RateVector& flux,
339 RateVector& darcy,
340 const IntensiveQuantities& intQuantsIn,
341 const IntensiveQuantities& intQuantsEx,
342 const unsigned& globalIndexIn,
343 const unsigned& globalIndexEx,
344 const ResidualNBInfo& nbInfo,
345 const ModuleParams& moduleParams)
346 {
347 OPM_TIMEBLOCK_LOCAL(calculateFluxes, Subsystem::Assembly);
348 const Scalar Vin = nbInfo.Vin;
349 const Scalar Vex = nbInfo.Vex;
350 const Scalar distZg = nbInfo.dZg;
351 const Scalar thpres = nbInfo.thpres;
352 const Scalar trans = nbInfo.trans;
353 const Scalar faceArea = nbInfo.faceArea;
354 FaceDir::DirEnum facedir = nbInfo.faceDir;
355
356 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
357 if (!FluidSystem::phaseIsActive(phaseIdx)) {
358 continue;
359 }
360 // darcy flux calculation
361 short dnIdx;
362 //
363 short upIdx;
364 // fake intices should only be used to get upwind anc compatibility with old functions
365 short interiorDofIdx = 0; // NB
366 short exteriorDofIdx = 1; // NB
367 Evaluation pressureDifference;
368 ExtensiveQuantities::calculatePhasePressureDiff_(upIdx,
369 dnIdx,
370 pressureDifference,
371 intQuantsIn,
372 intQuantsEx,
373 phaseIdx, // input
374 interiorDofIdx, // input
375 exteriorDofIdx, // input
376 Vin,
377 Vex,
378 globalIndexIn,
379 globalIndexEx,
380 distZg,
381 thpres,
382 moduleParams);
383
384 const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
385 unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
386 // Use arithmetic average (more accurate with harmonic, but that requires recomputing the transmissbility)
387 Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() +
388 Toolbox::value(intQuantsEx.rockCompTransMultiplier())) / 2;
389 if constexpr (enableBioeffects) {
390 transMult *= (intQuantsIn.permFactor() + Toolbox::value(intQuantsEx.permFactor())) / 2;
391 }
392 Evaluation darcyFlux;
393 if (globalUpIndex == globalIndexIn) {
394 darcyFlux = pressureDifference * up.mobility(phaseIdx, facedir) * transMult * (-trans / faceArea);
395 } else {
396 darcyFlux = pressureDifference *
397 (Toolbox::value(up.mobility(phaseIdx, facedir)) * transMult * (-trans / faceArea));
398 }
399
400 unsigned activeCompIdx =
401 FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
402 // NB! For the FLORES fluxes without derivatives
403 darcy[conti0EqIdx + activeCompIdx] = darcyFlux.value() * faceArea;
404
405 unsigned pvtRegionIdx = up.pvtRegionIndex();
406 // if (upIdx == globalFocusDofIdx){
407 if (globalUpIndex == globalIndexIn) {
408 const auto& invB
409 = getInvB_<FluidSystem, FluidState, Evaluation>(up.fluidState(), phaseIdx, pvtRegionIdx);
410 const auto& surfaceVolumeFlux = invB * darcyFlux;
411 evalPhaseFluxes_<Evaluation>(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
412 if constexpr (enableEnergy) {
413 EnergyModule::template
414 addPhaseEnthalpyFluxes_<Evaluation>(flux, phaseIdx, darcyFlux, up.fluidState());
415 }
416 if constexpr (enableBioeffects) {
417 BioeffectsModule::template
418 addBioeffectsFluxes_<Evaluation>(flux, phaseIdx, darcyFlux, up);
419 }
420 if constexpr (enableBrine) {
421 BrineModule::template
422 addBrineFluxes_<Evaluation, FluidState>(flux, phaseIdx, darcyFlux, up.fluidState());
423 }
424 } else {
425 const auto& invB = getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), phaseIdx, pvtRegionIdx);
426 const auto& surfaceVolumeFlux = invB * darcyFlux;
427 evalPhaseFluxes_<Scalar>(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
428 if constexpr (enableEnergy) {
429 EnergyModule::template
430 addPhaseEnthalpyFluxes_<Scalar>(flux, phaseIdx, darcyFlux, up.fluidState());
431 }
432 if constexpr (enableBioeffects) {
433 BioeffectsModule::template
434 addBioeffectsFluxes_<Scalar>(flux, phaseIdx, darcyFlux, up);
435 }
436 if constexpr (enableBrine) {
437 BrineModule::template
438 addBrineFluxes_<Scalar, FluidState>(flux, phaseIdx, darcyFlux, up.fluidState());
439 }
440 }
441 }
442
443 // deal with solvents (if present)
444 static_assert(!enableSolvent,
445 "Relevant computeFlux() method must be implemented for this module before enabling.");
446 // SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
447
448 // deal with zFracton (if present)
449 static_assert(!enableExtbo,
450 "Relevant computeFlux() method must be implemented for this module before enabling.");
451 // ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
452
453 // deal with polymer (if present)
454 static_assert(!enablePolymer,
455 "Relevant computeFlux() method must be implemented for this module before enabling.");
456 // PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
457
458 // deal with convective mixing
459 if constexpr (enableConvectiveMixing) {
460 ConvectiveMixingModule::addConvectiveMixingFlux(flux,
461 intQuantsIn,
462 intQuantsEx,
463 globalIndexIn,
464 globalIndexEx,
465 nbInfo.dZg,
466 nbInfo.trans,
467 nbInfo.faceArea,
468 moduleParams.convectiveMixingModuleParam);
469 }
470
471 // deal with energy (if present)
472 if constexpr (enableEnergy) {
473 const Scalar inAlpha = nbInfo.inAlpha;
474 const Scalar outAlpha = nbInfo.outAlpha;
475 Evaluation heatFlux;
476
477 short interiorDofIdx = 0; // NB
478 short exteriorDofIdx = 1; // NB
479
480 EnergyModule::ExtensiveQuantities::updateEnergy(heatFlux,
481 interiorDofIdx, // focusDofIndex,
482 interiorDofIdx,
483 exteriorDofIdx,
484 intQuantsIn,
485 intQuantsEx,
486 intQuantsIn.fluidState(),
487 intQuantsEx.fluidState(),
488 inAlpha,
489 outAlpha,
490 faceArea);
491 EnergyModule::addHeatFlux(flux, heatFlux);
492 }
493 // NB need to be tha last energy call since it does scaling
494 // EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
495
496 // deal with foam (if present)
497 static_assert(!enableFoam,
498 "Relevant computeFlux() method must be implemented for this module before enabling.");
499 // FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
500
501 // deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity).
502 if constexpr (enableDiffusion) {
503 typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
504 DiffusionModule::ExtensiveQuantities::update(effectiveDiffusionCoefficient, intQuantsIn, intQuantsEx);
505 const Scalar diffusivity = nbInfo.diffusivity;
506 const Scalar tmpdiffusivity = diffusivity / faceArea;
507 DiffusionModule::addDiffusiveFlux(flux,
508 intQuantsIn,
509 intQuantsEx,
510 tmpdiffusivity,
511 effectiveDiffusionCoefficient);
512 }
513
514 // deal with dispersion (if present). opm-models expects per area flux (added in the tmpdispersivity).
515 if constexpr (enableDispersion) {
516 typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg;
517 DispersionModule::ExtensiveQuantities::update(normVelocityAvg, intQuantsIn, intQuantsEx);
518 const Scalar dispersivity = nbInfo.dispersivity;
519 const Scalar tmpdispersivity = dispersivity / faceArea;
520 DispersionModule::addDispersiveFlux(flux,
521 intQuantsIn,
522 intQuantsEx,
523 tmpdispersivity,
524 normVelocityAvg);
525 }
526
527 // apply the scaling for the urea equation in MICP
528 if constexpr (enableMICP) {
530 }
531 }
532
533 template <class BoundaryConditionData>
534 static void computeBoundaryFlux(RateVector& bdyFlux,
535 const Problem& problem,
536 const BoundaryConditionData& bdyInfo,
537 const IntensiveQuantities& insideIntQuants,
538 unsigned globalSpaceIdx)
539 {
540 switch (bdyInfo.type) {
541 case BCType::NONE:
542 bdyFlux = 0.0;
543 break;
544 case BCType::RATE:
545 computeBoundaryFluxRate(bdyFlux, bdyInfo);
546 break;
547 case BCType::FREE:
548 case BCType::DIRICHLET:
549 computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
550 break;
551 case BCType::THERMAL:
552 computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
553 break;
554 default:
555 throw std::logic_error("Unknown boundary condition type " +
556 std::to_string(static_cast<int>(bdyInfo.type)) +
557 " in computeBoundaryFlux()." );
558 }
559 }
560
561 template <class BoundaryConditionData>
562 static void computeBoundaryFluxRate(RateVector& bdyFlux,
563 const BoundaryConditionData& bdyInfo)
564 {
565 bdyFlux.setMassRate(bdyInfo.massRate, bdyInfo.pvtRegionIdx);
566 }
567
568 template <class BoundaryConditionData>
569 static void computeBoundaryFluxFree(const Problem& problem,
570 RateVector& bdyFlux,
571 const BoundaryConditionData& bdyInfo,
572 const IntensiveQuantities& insideIntQuants,
573 unsigned globalSpaceIdx)
574 {
575 OPM_TIMEBLOCK_LOCAL(computeBoundaryFluxFree, Subsystem::Assembly);
576 std::array<short, numPhases> upIdx;
577 std::array<short, numPhases> dnIdx;
578 std::array<Evaluation, numPhases> volumeFlux;
579 std::array<Evaluation, numPhases> pressureDifference;
580
581 ExtensiveQuantities::calculateBoundaryGradients_(problem,
582 globalSpaceIdx,
583 insideIntQuants,
584 bdyInfo.boundaryFaceIndex,
585 bdyInfo.faceArea,
586 bdyInfo.faceZCoord,
587 bdyInfo.exFluidState,
588 upIdx,
589 dnIdx,
590 volumeFlux,
591 pressureDifference);
592
594 // advective fluxes of all components in all phases
596 bdyFlux = 0.0;
597 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
598 if (!FluidSystem::phaseIsActive(phaseIdx)) {
599 continue;
600 }
601 const auto& pBoundary = bdyInfo.exFluidState.pressure(phaseIdx);
602 const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx);
603 const unsigned pvtRegionIdx = insideIntQuants.pvtRegionIndex();
604
605 RateVector tmp(0.0);
606 const auto& darcyFlux = volumeFlux[phaseIdx];
607 // mass conservation
608 if (pBoundary < pInside) {
609 // outflux
610 const auto& invB =
611 getInvB_<FluidSystem, FluidState, Evaluation>(insideIntQuants.fluidState(), phaseIdx, pvtRegionIdx);
612 Evaluation surfaceVolumeFlux = invB * darcyFlux;
613 evalPhaseFluxes_<Evaluation>(tmp,
614 phaseIdx,
615 insideIntQuants.pvtRegionIndex(),
616 surfaceVolumeFlux,
617 insideIntQuants.fluidState());
618 if constexpr (enableEnergy) {
619 EnergyModule::template
620 addPhaseEnthalpyFluxes_<Evaluation>(tmp, phaseIdx, darcyFlux, insideIntQuants.fluidState());
621 }
622 } else if (pBoundary > pInside) {
623 // influx
624 using ScalarFluidState = decltype(bdyInfo.exFluidState);
625 const auto& invB =
626 getInvB_<FluidSystem, ScalarFluidState, Scalar>(bdyInfo.exFluidState, phaseIdx, pvtRegionIdx);
627 Evaluation surfaceVolumeFlux = invB * darcyFlux;
628 evalPhaseFluxes_<Scalar>(tmp,
629 phaseIdx,
630 insideIntQuants.pvtRegionIndex(),
631 surfaceVolumeFlux,
632 bdyInfo.exFluidState);
633 if constexpr (enableEnergy) {
634 EnergyModule::template
635 addPhaseEnthalpyFluxes_<Scalar>(tmp, phaseIdx, darcyFlux, bdyInfo.exFluidState);
636 }
637 }
638
639 for (unsigned i = 0; i < tmp.size(); ++i) {
640 bdyFlux[i] += tmp[i];
641 }
642 }
643
644 // conductive heat flux from boundary
645 if constexpr (enableEnergy) {
646 Evaluation heatFlux;
647 // avoid overload of functions with same number of elements in eclproblem
648 Scalar alpha =
649 problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
650 unsigned inIdx = 0; // dummy
651 // always calculated with derivatives of this cell
652 EnergyModule::ExtensiveQuantities::updateEnergyBoundary(heatFlux,
653 insideIntQuants,
654 /*focusDofIndex*/ inIdx,
655 inIdx,
656 alpha,
657 bdyInfo.exFluidState);
658 EnergyModule::addHeatFlux(bdyFlux, heatFlux);
659 }
660
661 static_assert(!enableSolvent,
662 "Relevant treatment of boundary conditions must be implemented before enabling.");
663 static_assert(!enablePolymer,
664 "Relevant treatment of boundary conditions must be implemented before enabling.");
665
666 // make sure that the right mass conservation quantities are used
667 adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex());
668
669#ifndef NDEBUG
670 for (unsigned i = 0; i < numEq; ++i) {
671 Valgrind::CheckDefined(bdyFlux[i]);
672 }
673 Valgrind::CheckDefined(bdyFlux);
674#endif
675 }
676
677 template <class BoundaryConditionData>
678 static void computeBoundaryThermal(const Problem& problem,
679 RateVector& bdyFlux,
680 const BoundaryConditionData& bdyInfo,
681 const IntensiveQuantities& insideIntQuants,
682 [[maybe_unused]] unsigned globalSpaceIdx)
683 {
684 OPM_TIMEBLOCK_LOCAL(computeBoundaryThermal, Subsystem::Assembly);
685 // only heat is allowed to flow through this boundary
686 bdyFlux = 0.0;
687
688 // conductive heat flux from boundary
689 if constexpr (enableEnergy) {
690 Evaluation heatFlux;
691 // avoid overload of functions with same numeber of elements in eclproblem
692 Scalar alpha =
693 problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
694 unsigned inIdx = 0;//dummy
695 // always calculated with derivatives of this cell
696 EnergyModule::ExtensiveQuantities::updateEnergyBoundary(heatFlux,
697 insideIntQuants,
698 /*focusDofIndex*/ inIdx,
699 inIdx,
700 alpha,
701 bdyInfo.exFluidState);
702 EnergyModule::addHeatFlux(bdyFlux, heatFlux);
703 }
704
705#ifndef NDEBUG
706 for (unsigned i = 0; i < numEq; ++i) {
707 Valgrind::CheckDefined(bdyFlux[i]);
708 }
709 Valgrind::CheckDefined(bdyFlux);
710#endif
711 }
712
713 static void computeSource(RateVector& source,
714 const Problem& problem,
715 const IntensiveQuantities& insideIntQuants,
716 unsigned globalSpaceIdex,
717 unsigned timeIdx)
718 {
719 OPM_TIMEBLOCK_LOCAL(computeSource, Subsystem::Assembly);
720 // retrieve the source term intrinsic to the problem
721 problem.source(source, globalSpaceIdex, timeIdx);
722
723 // deal with bioeffects (if present)
724 BioeffectsModule::addSource(source, problem, insideIntQuants, globalSpaceIdex);
725
726 // scale the source term of the energy equation
727 if constexpr (enableEnergy) {
728 source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
729 }
730 }
731
732 static void computeSourceDense(RateVector& source,
733 const Problem& problem,
734 const IntensiveQuantities& insideIntQuants,
735 unsigned globalSpaceIdex,
736 unsigned timeIdx)
737 {
738 source = 0.0;
739 problem.addToSourceDense(source, globalSpaceIdex, timeIdx);
740
741 // deal with bioeffects (if present)
742 BioeffectsModule::addSource(source, problem, insideIntQuants, globalSpaceIdex);
743
744 // scale the source term of the energy equation
745 if constexpr (enableEnergy) {
746 source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
747 }
748 }
749
753 void computeSource(RateVector& source,
754 const ElementContext& elemCtx,
755 unsigned dofIdx,
756 unsigned timeIdx) const
757 {
758 OPM_TIMEBLOCK_LOCAL(computeSource, Subsystem::Assembly);
759 // retrieve the source term intrinsic to the problem
760 elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
761
762 // deal with bioeffects (if present)
763 BioeffectsModule::addSource(source, elemCtx, dofIdx, timeIdx);
764
765 // scale the source term of the energy equation
766 if constexpr (enableEnergy) {
767 source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
768 }
769 }
770
771 template <class UpEval, class FluidState>
772 static void evalPhaseFluxes_(RateVector& flux,
773 unsigned phaseIdx,
774 unsigned pvtRegionIdx,
775 const ExtensiveQuantities& extQuants,
776 const FluidState& upFs)
777 {
778 const auto& invB = getInvB_<FluidSystem, FluidState, UpEval>(upFs, phaseIdx, pvtRegionIdx);
779 const auto& surfaceVolumeFlux = invB * extQuants.volumeFlux(phaseIdx);
780 evalPhaseFluxes_<UpEval>(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, upFs);
781 }
782
787 template <class UpEval, class Eval,class FluidState>
788 static void evalPhaseFluxes_(RateVector& flux,
789 unsigned phaseIdx,
790 unsigned pvtRegionIdx,
791 const Eval& surfaceVolumeFlux,
792 const FluidState& upFs)
793 {
794 unsigned activeCompIdx =
795 FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
796
797 if constexpr (blackoilConserveSurfaceVolume) {
798 flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux;
799 }
800 else {
801 flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux *
802 FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
803 }
804
805 if (phaseIdx == oilPhaseIdx) {
806 // dissolved gas (in the oil phase).
807 if (FluidSystem::enableDissolvedGas()) {
808 const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
809
810 const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
811 if constexpr (blackoilConserveSurfaceVolume) {
812 flux[conti0EqIdx + activeGasCompIdx] += Rs * surfaceVolumeFlux;
813 }
814 else {
815 flux[conti0EqIdx + activeGasCompIdx] +=
816 Rs * surfaceVolumeFlux *
817 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
818 }
819 }
820 }
821 else if (phaseIdx == waterPhaseIdx) {
822 // dissolved gas (in the water phase).
823 if (FluidSystem::enableDissolvedGasInWater()) {
824 const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
825
826 const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
827 if constexpr (blackoilConserveSurfaceVolume) {
828 flux[conti0EqIdx + activeGasCompIdx] += Rsw * surfaceVolumeFlux;
829 }
830 else {
831 flux[conti0EqIdx + activeGasCompIdx] +=
832 Rsw * surfaceVolumeFlux *
833 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
834 }
835 }
836 }
837 else if (phaseIdx == gasPhaseIdx) {
838 // vaporized oil (in the gas phase).
839 if (FluidSystem::enableVaporizedOil()) {
840 const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
841
842 const unsigned activeOilCompIdx = FluidSystem::canonicalToActiveCompIdx(oilCompIdx);
843 if constexpr (blackoilConserveSurfaceVolume) {
844 flux[conti0EqIdx + activeOilCompIdx] += Rv * surfaceVolumeFlux;
845 }
846 else {
847 flux[conti0EqIdx + activeOilCompIdx] +=
848 Rv * surfaceVolumeFlux *
849 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
850 }
851 }
852 // vaporized water (in the gas phase).
853 if (FluidSystem::enableVaporizedWater()) {
854 const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
855
856 const unsigned activeWaterCompIdx = FluidSystem::canonicalToActiveCompIdx(waterCompIdx);
857 if constexpr (blackoilConserveSurfaceVolume) {
858 flux[conti0EqIdx + activeWaterCompIdx] += Rvw * surfaceVolumeFlux;
859 }
860 else {
861 flux[conti0EqIdx + activeWaterCompIdx] +=
862 Rvw * surfaceVolumeFlux *
863 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
864 }
865 }
866 }
867 }
868
880 template <class Scalar>
881 static void adaptMassConservationQuantities_(Dune::FieldVector<Scalar, numEq>& container,
882 unsigned pvtRegionIdx)
883 {
884 if constexpr (!blackoilConserveSurfaceVolume) {
885 // convert "surface volume" to mass. this is complicated a bit by the fact that
886 // not all phases are necessarily enabled. (we here assume that if a fluid phase
887 // is disabled, its respective "main" component is not considered as well.)
888
889 if constexpr (waterEnabled) {
890 const unsigned activeWaterCompIdx = FluidSystem::canonicalToActiveCompIdx(waterCompIdx);
891 container[conti0EqIdx + activeWaterCompIdx] *=
892 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
893 }
894
895 if constexpr (gasEnabled) {
896 const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
897 container[conti0EqIdx + activeGasCompIdx] *=
898 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
899 }
900
901 if constexpr (oilEnabled) {
902 const unsigned activeOilCompIdx = FluidSystem::canonicalToActiveCompIdx(oilCompIdx);
903 container[conti0EqIdx + activeOilCompIdx] *=
904 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
905 }
906 }
907 }
908
909 // NNC does not have a direction
910 static FaceDir::DirEnum faceDirFromDirId(const int dirId)
911 { return dirId < 0 ? FaceDir::DirEnum::Unknown : FaceDir::FromIntersectionIndex(dirId); }
912};
913
914} // namespace Opm
915
916#endif
Contains the classes required to extend the black-oil model by bioeffects.
Contains the classes required to extend the black-oil model by brine.
Classes required for dynamic convective mixing.
Classes required for molecular diffusion.
Classes required for mechanical dispersion.
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 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 bioeffects.
Definition: blackoilbioeffectsmodules.hh:93
static void applyScaling(RateVector &flux)
Definition: blackoilbioeffectsmodules.hh:235
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilbioeffectsmodules.hh:177
static void addSource(RateVector &source, const Problem &problem, const IntensiveQuantities &intQuants, unsigned globalSpaceIdex)
Definition: blackoilbioeffectsmodules.hh:276
Contains the high level supplements required to extend the black oil model by brine.
Definition: blackoilbrinemodules.hh:56
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilbrinemodules.hh:163
Definition: blackoilconvectivemixingmodule.hh:64
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: blackoildiffusionmodule.hh:50
Provides the auxiliary methods required for consideration of the dispersion equation.
Definition: blackoildispersionmodule.hh:58
Contains the high level supplements required to extend the black oil model by energy.
Definition: blackoilenergymodules.hh:58
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilenergymodules.hh:153
static void addHeatFlux(RateVector &flux, const Evaluation &heatFlux)
Definition: blackoilenergymodules.hh:211
Contains the high level supplements required to extend the black oil model.
Definition: blackoilextbomodules.hh:62
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilextbomodules.hh:156
Contains the high level supplements required to extend the black oil model to include the effects of ...
Definition: blackoilfoammodules.hh:58
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilfoammodules.hh:164
Calculates the local residual of the black oil model.
Definition: blackoillocalresidualtpfa.hh:63
void computeSource(RateVector &source, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx) const
Calculate the source term of the equation.
Definition: blackoillocalresidualtpfa.hh:753
static void computeBoundaryFlux(RateVector &bdyFlux, const Problem &problem, const BoundaryConditionData &bdyInfo, const IntensiveQuantities &insideIntQuants, unsigned globalSpaceIdx)
Definition: blackoillocalresidualtpfa.hh:534
static void computeBoundaryFluxFree(const Problem &problem, RateVector &bdyFlux, const BoundaryConditionData &bdyInfo, const IntensiveQuantities &insideIntQuants, unsigned globalSpaceIdx)
Definition: blackoillocalresidualtpfa.hh:569
static void computeBoundaryFluxRate(RateVector &bdyFlux, const BoundaryConditionData &bdyInfo)
Definition: blackoillocalresidualtpfa.hh:562
static FaceDir::DirEnum faceDirFromDirId(const int dirId)
Definition: blackoillocalresidualtpfa.hh:910
static void evalPhaseFluxes_(RateVector &flux, unsigned phaseIdx, unsigned pvtRegionIdx, const ExtensiveQuantities &extQuants, const FluidState &upFs)
Definition: blackoillocalresidualtpfa.hh:772
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoillocalresidualtpfa.hh:270
static void computeStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoillocalresidualtpfa.hh:161
static void computeFlux(RateVector &flux, RateVector &darcy, const unsigned globalIndexIn, const unsigned globalIndexEx, const IntensiveQuantities &intQuantsIn, const IntensiveQuantities &intQuantsEx, const ResidualNBInfo &nbInfo, const ModuleParams &moduleParams)
Definition: blackoillocalresidualtpfa.hh:244
static void computeBoundaryThermal(const Problem &problem, RateVector &bdyFlux, const BoundaryConditionData &bdyInfo, const IntensiveQuantities &insideIntQuants, unsigned globalSpaceIdx)
Definition: blackoillocalresidualtpfa.hh:678
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: blackoillocalresidualtpfa.hh:881
static void evalPhaseFluxes_(RateVector &flux, unsigned phaseIdx, unsigned pvtRegionIdx, const Eval &surfaceVolumeFlux, const FluidState &upFs)
Helper function to calculate the flux of mass in terms of conservation quantities via specific fluid ...
Definition: blackoillocalresidualtpfa.hh:788
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: blackoillocalresidualtpfa.hh:151
static void calculateFluxes_(RateVector &flux, RateVector &darcy, const IntensiveQuantities &intQuantsIn, const IntensiveQuantities &intQuantsEx, const unsigned &globalIndexIn, const unsigned &globalIndexEx, const ResidualNBInfo &nbInfo, const ModuleParams &moduleParams)
Definition: blackoillocalresidualtpfa.hh:338
static void computeSource(RateVector &source, const Problem &problem, const IntensiveQuantities &insideIntQuants, unsigned globalSpaceIdex, unsigned timeIdx)
Definition: blackoillocalresidualtpfa.hh:713
static void computeSourceDense(RateVector &source, const Problem &problem, const IntensiveQuantities &insideIntQuants, unsigned globalSpaceIdex, unsigned timeIdx)
Definition: blackoillocalresidualtpfa.hh:732
Contains the high level supplements required to extend the black oil model by polymer.
Definition: blackoilpolymermodules.hh:64
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilpolymermodules.hh:236
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilsolventmodules.hh:68
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilsolventmodules.hh:183
@ NONE
Definition: DeferredLogger.hpp:46
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
std::string to_string(const ConvergenceReport::ReservoirFailure::Type t)
Definition: blackoillocalresidualtpfa.hh:143
ConvectiveMixingModuleParam convectiveMixingModuleParam
Definition: blackoillocalresidualtpfa.hh:144
Definition: blackoillocalresidualtpfa.hh:128
FaceDir::DirEnum faceDir
Definition: blackoillocalresidualtpfa.hh:133
double faceArea
Definition: blackoillocalresidualtpfa.hh:130
ConditionalStorage< enableEnergy, double > inAlpha
Definition: blackoillocalresidualtpfa.hh:136
ConditionalStorage< enableDiffusion, double > diffusivity
Definition: blackoillocalresidualtpfa.hh:138
double dZg
Definition: blackoillocalresidualtpfa.hh:132
double Vin
Definition: blackoillocalresidualtpfa.hh:134
ConditionalStorage< enableDispersion, double > dispersivity
Definition: blackoillocalresidualtpfa.hh:139
double thpres
Definition: blackoillocalresidualtpfa.hh:131
ConditionalStorage< enableEnergy, double > outAlpha
Definition: blackoillocalresidualtpfa.hh:137
double trans
Definition: blackoillocalresidualtpfa.hh:129
double Vex
Definition: blackoillocalresidualtpfa.hh:135