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);
147 computeStorage(storage,
148 intQuants);
149 }
150
151 template <class LhsEval>
152 static void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
153 const IntensiveQuantities& intQuants)
154 {
155 OPM_TIMEBLOCK_LOCAL(computeStorage);
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));
165 LhsEval surfaceVolume =
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())
177 * surfaceVolume;
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())
185 * surfaceVolume;
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())
193 * surfaceVolume;
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())
201 * surfaceVolume;
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 {
243 OPM_TIMEBLOCK_LOCAL(computeFlux);
244 flux = 0.0;
245 darcy = 0.0;
246
247 calculateFluxes_(flux,
248 darcy,
249 intQuantsIn,
250 intQuantsEx,
251 globalIndexIn,
252 globalIndexEx,
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 {
265 OPM_TIMEBLOCK_LOCAL(computeFlux);
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();
277 assert(interiorDofIdx != exteriorDofIdx);
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,
317 intQuantsIn,
318 intQuantsEx,
319 globalIndexIn,
320 globalIndexEx,
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 {
334 OPM_TIMEBLOCK_LOCAL(calculateFluxes);
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,
357 intQuantsIn,
358 intQuantsEx,
359 phaseIdx, // input
360 interiorDofIdx, // input
361 exteriorDofIdx, // input
362 Vin,
363 Vex,
364 globalIndexIn,
365 globalIndexEx,
366 distZg,
367 thpres,
368 moduleParams);
369
370
371
372 const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
373 unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
374 // Use arithmetic average (more accurate with harmonic, but that requires recomputing the transmissbility)
375 const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;
376 Evaluation darcyFlux;
377 if (globalUpIndex == globalIndexIn) {
378 darcyFlux = pressureDifference * up.mobility(phaseIdx, facedir) * transMult * (-trans / faceArea);
379 } else {
380 darcyFlux = pressureDifference *
381 (Toolbox::value(up.mobility(phaseIdx, facedir)) * transMult * (-trans / faceArea));
382 }
383
384 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
385 darcy[conti0EqIdx + activeCompIdx] = darcyFlux.value() * faceArea; // NB! For the FLORES fluxes without derivatives
386
387 unsigned pvtRegionIdx = up.pvtRegionIndex();
388 // if (upIdx == globalFocusDofIdx){
389 if (globalUpIndex == globalIndexIn) {
390 const auto& invB
391 = getInvB_<FluidSystem, FluidState, Evaluation>(up.fluidState(), phaseIdx, pvtRegionIdx);
392 const auto& surfaceVolumeFlux = invB * darcyFlux;
393 evalPhaseFluxes_<Evaluation, Evaluation, FluidState>(
394 flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
395 if constexpr (enableEnergy) {
396 EnergyModule::template addPhaseEnthalpyFluxes_<Evaluation, Evaluation, FluidState>(
397 flux, phaseIdx, darcyFlux, up.fluidState());
398 }
399 } else {
400 const auto& invB = getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), phaseIdx, pvtRegionIdx);
401 const auto& surfaceVolumeFlux = invB * darcyFlux;
402 evalPhaseFluxes_<Scalar, Evaluation, FluidState>(
403 flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
404 if constexpr (enableEnergy) {
405 EnergyModule::template
406 addPhaseEnthalpyFluxes_<Scalar, Evaluation, FluidState>
407 (flux,phaseIdx,darcyFlux, up.fluidState());
408 }
409 }
410
411 }
412
413 // deal with solvents (if present)
414 static_assert(!enableSolvent, "Relevant computeFlux() method must be implemented for this module before enabling.");
415 // SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
416
417 // deal with zFracton (if present)
418 static_assert(!enableExtbo, "Relevant computeFlux() method must be implemented for this module before enabling.");
419 // ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
420
421 // deal with polymer (if present)
422 static_assert(!enablePolymer, "Relevant computeFlux() method must be implemented for this module before enabling.");
423 // PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
424
425 // deal with convective mixing
426 if constexpr(enableConvectiveMixing) {
427 ConvectiveMixingModule::addConvectiveMixingFlux(flux,
428 intQuantsIn,
429 intQuantsEx,
430 globalIndexIn,
431 globalIndexEx,
432 nbInfo.dZg,
433 nbInfo.trans,
434 nbInfo.faceArea,
435 moduleParams.convectiveMixingModuleParam);
436 }
437
438 // deal with energy (if present)
439 if constexpr(enableEnergy){
440 const Scalar inAlpha = nbInfo.inAlpha;
441 const Scalar outAlpha = nbInfo.outAlpha;
442 Evaluation heatFlux;
443 {
444 short interiorDofIdx = 0; // NB
445 short exteriorDofIdx = 1; // NB
446
447 EnergyModule::ExtensiveQuantities::template updateEnergy(heatFlux,
448 interiorDofIdx, // focusDofIndex,
449 interiorDofIdx,
450 exteriorDofIdx,
451 intQuantsIn,
452 intQuantsEx,
453 intQuantsIn.fluidState(),
454 intQuantsEx.fluidState(),
455 inAlpha,
456 outAlpha,
457 faceArea);
458 }
459 EnergyModule::addHeatFlux(flux, heatFlux);
460 }
461 // NB need to be tha last energy call since it does scaling
462 // EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
463
464 // deal with foam (if present)
465 static_assert(!enableFoam, "Relevant computeFlux() method must be implemented for this module before enabling.");
466 // FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
467
468 // deal with salt (if present)
469 static_assert(!enableBrine, "Relevant computeFlux() method must be implemented for this module before enabling.");
470 // BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
471
472
473
474 // deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity).
475 if constexpr(enableDiffusion){
476 typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
477 DiffusionModule::ExtensiveQuantities::update(effectiveDiffusionCoefficient, intQuantsIn, intQuantsEx);
478 const Scalar diffusivity = nbInfo.diffusivity;
479 const Scalar tmpdiffusivity = diffusivity / faceArea;
480 DiffusionModule::addDiffusiveFlux(flux,
481 intQuantsIn.fluidState(),
482 intQuantsEx.fluidState(),
483 tmpdiffusivity,
484 effectiveDiffusionCoefficient);
485
486 }
487 // deal with dispersion (if present). opm-models expects per area flux (added in the tmpdispersivity).
488 if constexpr(enableDispersion){
489 typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg;
490 DispersionModule::ExtensiveQuantities::update(normVelocityAvg, intQuantsIn, intQuantsEx);
491 const Scalar dispersivity = nbInfo.dispersivity;
492 const Scalar tmpdispersivity = dispersivity / faceArea;
493 DispersionModule::addDispersiveFlux(flux,
494 intQuantsIn.fluidState(),
495 intQuantsEx.fluidState(),
496 tmpdispersivity,
497 normVelocityAvg);
498
499 }
500 // deal with micp (if present)
501 static_assert(!enableMICP, "Relevant computeFlux() method must be implemented for this module before enabling.");
502 // MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
503
504 }
505
506 template <class BoundaryConditionData>
507 static void computeBoundaryFlux(RateVector& bdyFlux,
508 const Problem& problem,
509 const BoundaryConditionData& bdyInfo,
510 const IntensiveQuantities& insideIntQuants,
511 unsigned globalSpaceIdx)
512 {
513 if (bdyInfo.type == BCType::NONE) {
514 bdyFlux = 0.0;
515 } else if (bdyInfo.type == BCType::RATE) {
516 computeBoundaryFluxRate(bdyFlux, bdyInfo);
517 } else if (bdyInfo.type == BCType::FREE || bdyInfo.type == BCType::DIRICHLET) {
518 computeBoundaryFluxFree(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
519 } else if (bdyInfo.type == BCType::THERMAL) {
520 computeBoundaryThermal(problem, bdyFlux, bdyInfo, insideIntQuants, globalSpaceIdx);
521 } else {
522 throw std::logic_error("Unknown boundary condition type " + std::to_string(static_cast<int>(bdyInfo.type)) + " in computeBoundaryFlux()." );
523 }
524 }
525
526 template <class BoundaryConditionData>
527 static void computeBoundaryFluxRate(RateVector& bdyFlux,
528 const BoundaryConditionData& bdyInfo)
529 {
530 bdyFlux.setMassRate(bdyInfo.massRate, bdyInfo.pvtRegionIdx);
531 }
532
533 template <class BoundaryConditionData>
534 static void computeBoundaryFluxFree(const Problem& problem,
535 RateVector& bdyFlux,
536 const BoundaryConditionData& bdyInfo,
537 const IntensiveQuantities& insideIntQuants,
538 unsigned globalSpaceIdx)
539 {
540 OPM_TIMEBLOCK_LOCAL(computeBoundaryFluxFree);
541 std::array<short, numPhases> upIdx;
542 std::array<short, numPhases> dnIdx;
543 std::array<Evaluation, numPhases> volumeFlux;
544 std::array<Evaluation, numPhases> pressureDifference;
545
546 ExtensiveQuantities::calculateBoundaryGradients_(problem,
547 globalSpaceIdx,
548 insideIntQuants,
549 bdyInfo.boundaryFaceIndex,
550 bdyInfo.faceArea,
551 bdyInfo.faceZCoord,
552 bdyInfo.exFluidState,
553 upIdx,
554 dnIdx,
555 volumeFlux,
556 pressureDifference);
557
559 // advective fluxes of all components in all phases
561 bdyFlux = 0.0;
562 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
563 if (!FluidSystem::phaseIsActive(phaseIdx)) {
564 continue;
565 }
566 const auto& pBoundary = bdyInfo.exFluidState.pressure(phaseIdx);
567 const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx);
568 const unsigned pvtRegionIdx = insideIntQuants.pvtRegionIndex();
569
570 RateVector tmp(0.0);
571 const auto& darcyFlux = volumeFlux[phaseIdx];
572 // mass conservation
573 if (pBoundary < pInside) {
574 // outflux
575 const auto& invB = getInvB_<FluidSystem, FluidState, Evaluation>(insideIntQuants.fluidState(), phaseIdx, pvtRegionIdx);
576 Evaluation surfaceVolumeFlux = invB * darcyFlux;
577 evalPhaseFluxes_<Evaluation>(tmp,
578 phaseIdx,
579 insideIntQuants.pvtRegionIndex(),
580 surfaceVolumeFlux,
581 insideIntQuants.fluidState());
582 if constexpr (enableEnergy) {
583 EnergyModule::template
584 addPhaseEnthalpyFluxes_<Evaluation, Evaluation, FluidState>
585 (tmp, phaseIdx, darcyFlux, insideIntQuants.fluidState());
586 }
587 } else if (pBoundary > pInside) {
588 // influx
589 using ScalarFluidState = decltype(bdyInfo.exFluidState);
590 const auto& invB = getInvB_<FluidSystem, ScalarFluidState, Scalar>(bdyInfo.exFluidState, phaseIdx, pvtRegionIdx);
591 Evaluation surfaceVolumeFlux = invB * darcyFlux;
592 evalPhaseFluxes_<Scalar>(tmp,
593 phaseIdx,
594 insideIntQuants.pvtRegionIndex(),
595 surfaceVolumeFlux,
596 bdyInfo.exFluidState);
597 if constexpr (enableEnergy) {
598 EnergyModule::template
599 addPhaseEnthalpyFluxes_<Scalar, Evaluation, ScalarFluidState>
600 (tmp,
601 phaseIdx,
602 darcyFlux,
603 bdyInfo.exFluidState);
604 }
605 }
606
607 for (unsigned i = 0; i < tmp.size(); ++i) {
608 bdyFlux[i] += tmp[i];
609 }
610 }
611
612 // conductive heat flux from boundary
613 if constexpr(enableEnergy){
614 Evaluation heatFlux;
615 // avoid overload of functions with same numeber of elements in eclproblem
616 Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
617 unsigned inIdx = 0;//dummy
618 // always calculated with derivatives of this cell
619 EnergyModule::ExtensiveQuantities::template updateEnergyBoundary(heatFlux,
620 insideIntQuants,
621 /*focusDofIndex*/ inIdx,
622 inIdx,
623 alpha,
624 bdyInfo.exFluidState);
625 EnergyModule::addHeatFlux(bdyFlux, heatFlux);
626 }
627
628 static_assert(!enableSolvent, "Relevant treatment of boundary conditions must be implemented before enabling.");
629 static_assert(!enablePolymer, "Relevant treatment of boundary conditions must be implemented before enabling.");
630 static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling.");
631
632 // make sure that the right mass conservation quantities are used
633 adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex());
634
635#ifndef NDEBUG
636 for (unsigned i = 0; i < numEq; ++i) {
637 Valgrind::CheckDefined(bdyFlux[i]);
638 }
639 Valgrind::CheckDefined(bdyFlux);
640#endif
641 }
642
643 template <class BoundaryConditionData>
644 static void computeBoundaryThermal(const Problem& problem,
645 RateVector& bdyFlux,
646 const BoundaryConditionData& bdyInfo,
647 const IntensiveQuantities& insideIntQuants,
648 [[maybe_unused]] unsigned globalSpaceIdx)
649 {
650 OPM_TIMEBLOCK_LOCAL(computeBoundaryThermal);
651 // only heat is allowed to flow through this boundary
652 bdyFlux = 0.0;
653
654 // conductive heat flux from boundary
655 if constexpr(enableEnergy){
656 Evaluation heatFlux;
657 // avoid overload of functions with same numeber of elements in eclproblem
658 Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
659 unsigned inIdx = 0;//dummy
660 // always calculated with derivatives of this cell
661 EnergyModule::ExtensiveQuantities::template updateEnergyBoundary(heatFlux,
662 insideIntQuants,
663 /*focusDofIndex*/ inIdx,
664 inIdx,
665 alpha,
666 bdyInfo.exFluidState);
667 EnergyModule::addHeatFlux(bdyFlux, heatFlux);
668 }
669
670#ifndef NDEBUG
671 for (unsigned i = 0; i < numEq; ++i) {
672 Valgrind::CheckDefined(bdyFlux[i]);
673 }
674 Valgrind::CheckDefined(bdyFlux);
675#endif
676 }
677
678 static void computeSource(RateVector& source,
679 const Problem& problem,
680 unsigned globalSpaceIdex,
681 unsigned timeIdx)
682 {
683 OPM_TIMEBLOCK_LOCAL(computeSource);
684 // retrieve the source term intrinsic to the problem
685 problem.source(source, globalSpaceIdex, timeIdx);
686
687 // deal with MICP (if present)
688 // deal with micp (if present)
689 static_assert(!enableMICP, "Relevant addSource() method must be implemented for this module before enabling.");
690 // MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
691
692 // scale the source term of the energy equation
693 if (enableEnergy)
694 source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
695 }
696
697 static void computeSourceDense(RateVector& source,
698 const Problem& problem,
699 unsigned globalSpaceIdex,
700 unsigned timeIdx)
701 {
702 source = 0.0;
703 problem.addToSourceDense(source, globalSpaceIdex, timeIdx);
704
705 // deal with MICP (if present)
706 // deal with micp (if present)
707 static_assert(!enableMICP, "Relevant addSource() method must be implemented for this module before enabling.");
708 // MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
709
710 // scale the source term of the energy equation
711 if (enableEnergy)
712 source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
713 }
714
718 void computeSource(RateVector& source,
719 const ElementContext& elemCtx,
720 unsigned dofIdx,
721 unsigned timeIdx) const
722 {
723 OPM_TIMEBLOCK_LOCAL(computeSource);
724 // retrieve the source term intrinsic to the problem
725 elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
726
727 // deal with MICP (if present)
728 MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
729
730 // scale the source term of the energy equation
731 if constexpr(enableEnergy)
732 source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
733 }
734
735 template <class UpEval, class FluidState>
736 static void evalPhaseFluxes_(RateVector& flux,
737 unsigned phaseIdx,
738 unsigned pvtRegionIdx,
739 const ExtensiveQuantities& extQuants,
740 const FluidState& upFs)
741 {
742
743 const auto& invB = getInvB_<FluidSystem, FluidState, UpEval>(upFs, phaseIdx, pvtRegionIdx);
744 const auto& surfaceVolumeFlux = invB * extQuants.volumeFlux(phaseIdx);
745 evalPhaseFluxes_<UpEval>(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, upFs);
746 }
747
752 template <class UpEval, class Eval,class FluidState>
753 static void evalPhaseFluxes_(RateVector& flux,
754 unsigned phaseIdx,
755 unsigned pvtRegionIdx,
756 const Eval& surfaceVolumeFlux,
757 const FluidState& upFs)
758 {
759 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
760
761 if (blackoilConserveSurfaceVolume)
762 flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux;
763 else
764 flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux*FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
765
766 if (phaseIdx == oilPhaseIdx) {
767 // dissolved gas (in the oil phase).
768 if (FluidSystem::enableDissolvedGas()) {
769 const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
770
771 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
772 if (blackoilConserveSurfaceVolume)
773 flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux;
774 else
775 flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
776 }
777 } else if (phaseIdx == waterPhaseIdx) {
778 // dissolved gas (in the water phase).
779 if (FluidSystem::enableDissolvedGasInWater()) {
780 const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
781
782 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
783 if (blackoilConserveSurfaceVolume)
784 flux[conti0EqIdx + activeGasCompIdx] += Rsw*surfaceVolumeFlux;
785 else
786 flux[conti0EqIdx + activeGasCompIdx] += Rsw*surfaceVolumeFlux*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
787 }
788 }
789 else if (phaseIdx == gasPhaseIdx) {
790 // vaporized oil (in the gas phase).
791 if (FluidSystem::enableVaporizedOil()) {
792 const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
793
794 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
795 if (blackoilConserveSurfaceVolume)
796 flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux;
797 else
798 flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux*FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
799 }
800 // vaporized water (in the gas phase).
801 if (FluidSystem::enableVaporizedWater()) {
802 const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
803
804 unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
805 if (blackoilConserveSurfaceVolume)
806 flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux;
807 else
808 flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux*FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
809 }
810 }
811 }
812
824 template <class Scalar>
825 static void adaptMassConservationQuantities_(Dune::FieldVector<Scalar, numEq>& container, unsigned pvtRegionIdx)
826 {
827 if (blackoilConserveSurfaceVolume)
828 return;
829
830 // convert "surface volume" to mass. this is complicated a bit by the fact that
831 // not all phases are necessarily enabled. (we here assume that if a fluid phase
832 // is disabled, its respective "main" component is not considered as well.)
833
834 if (waterEnabled) {
835 unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
836 container[conti0EqIdx + activeWaterCompIdx] *=
837 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
838 }
839
840 if (gasEnabled) {
841 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
842 container[conti0EqIdx + activeGasCompIdx] *=
843 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
844 }
845
846 if (oilEnabled) {
847 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
848 container[conti0EqIdx + activeOilCompIdx] *=
849 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
850 }
851 }
852
853
854 static FaceDir::DirEnum faceDirFromDirId(const int dirId)
855 {
856 // NNC does not have a direction
857 if (dirId < 0 ) {
858 return FaceDir::DirEnum::Unknown;
859 }
860 return FaceDir::FromIntersectionIndex(dirId);
861 }
862};
863
864} // namespace Opm
865
866#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.
Provides the auxiliary methods required for consideration of the dispersion equation.
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:718
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:825
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:753
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
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
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Definition blackoillocalresidualtpfa.hh:133
Definition blackoillocalresidualtpfa.hh:119