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