My Project
Loading...
Searching...
No Matches
OutputBlackoilModule.hpp
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*/
27#ifndef OPM_OUTPUT_BLACK_OIL_MODULE_HPP
28#define OPM_OUTPUT_BLACK_OIL_MODULE_HPP
29
30#include <dune/common/fvector.hh>
31
32#include <opm/simulators/utils/moduleVersion.hpp>
33
34#include <opm/common/Exceptions.hpp>
35#include <opm/common/TimingMacros.hpp>
36#include <opm/common/OpmLog/OpmLog.hpp>
37
38#include <opm/input/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
39
40#include <opm/material/common/Valgrind.hpp>
41#include <opm/material/fluidmatrixinteractions/EclEpsScalingPoints.hpp>
42#include <opm/material/fluidstates/BlackOilFluidState.hpp>
43#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
44
49
50#include <opm/output/data/Cells.hpp>
51#include <opm/output/eclipse/EclipseIO.hpp>
52#include <opm/output/eclipse/Inplace.hpp>
53
56
57#include <algorithm>
58#include <array>
59#include <cassert>
60#include <cstddef>
61#include <functional>
62#include <stdexcept>
63#include <string>
64#include <type_traits>
65#include <utility>
66#include <vector>
67
68namespace Opm::Parameters {
69
70struct ForceDisableFluidInPlaceOutput { static constexpr bool value = false; };
71struct ForceDisableResvFluidInPlaceOutput { static constexpr bool value = false; };
72
73} // namespace Opm::Parameters
74
75namespace Opm
76{
77
78// forward declaration
79template <class TypeTag>
81
88template <class TypeTag>
89class OutputBlackOilModule : public GenericOutputBlackoilModule<GetPropType<TypeTag, Properties::FluidSystem>>
90{
101 using Element = typename GridView::template Codim<0>::Entity;
102 using ElementIterator = typename GridView::template Codim<0>::Iterator;
105 using Dir = FaceDir::DirEnum;
106
107 enum { conti0EqIdx = Indices::conti0EqIdx };
108 enum { numPhases = FluidSystem::numPhases };
109 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
110 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
111 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
112 enum { gasCompIdx = FluidSystem::gasCompIdx };
113 enum { oilCompIdx = FluidSystem::oilCompIdx };
114 enum { waterCompIdx = FluidSystem::waterCompIdx };
115 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
116
117public:
118 template <class CollectDataToIORankType>
119 OutputBlackOilModule(const Simulator& simulator,
120 const SummaryConfig& smryCfg,
122 : BaseType(simulator.vanguard().eclState(),
123 simulator.vanguard().schedule(),
124 smryCfg,
125 simulator.vanguard().summaryState(),
137 , simulator_(simulator)
138 {
139 for (auto& region_pair : this->regions_) {
140 this->createLocalRegion_(region_pair.second);
141 }
142
143 auto isCartIdxOnThisRank = [&collectToIORank](const int idx) {
144 return collectToIORank.isCartIdxOnThisRank(idx);
145 };
146
147 this->setupBlockData(isCartIdxOnThisRank);
148
149 this->forceDisableFipOutput_ =
150 Parameters::Get<Parameters::ForceDisableFluidInPlaceOutput>();
151
152 this->forceDisableFipresvOutput_ =
153 Parameters::Get<Parameters::ForceDisableResvFluidInPlaceOutput>();
154
155 if (! Parameters::Get<Parameters::OwnerCellsFirst>()) {
156 const std::string msg = "The output code does not support --owner-cells-first=false.";
157 if (collectToIORank.isIORank()) {
158 OpmLog::error(msg);
159 }
160 OPM_THROW_NOLOG(std::runtime_error, msg);
161 }
162
163 if (smryCfg.match("[FB]PP[OGW]") || smryCfg.match("RPP[OGW]*")) {
164 auto rset = this->eclState_.fieldProps().fip_regions();
165 rset.push_back("PVTNUM");
166
167 // Note: We explicitly use decltype(auto) here because the
168 // default scheme (-> auto) will deduce an undesirable type. We
169 // need the "reference to vector" semantics in this instance.
170 this->regionAvgDensity_
171 .emplace(this->simulator_.gridView().comm(),
172 FluidSystem::numPhases, rset,
173 [fp = std::cref(this->eclState_.fieldProps())]
174 (const std::string& rsetName) -> decltype(auto)
175 { return fp.get().get_int(rsetName); });
176 }
177 }
178
182 static void registerParameters()
183 {
184 Parameters::Register<Parameters::ForceDisableFluidInPlaceOutput>
185 ("Do not print fluid-in-place values after each report step "
186 "even if requested by the deck.");
187 Parameters::Register<Parameters::ForceDisableResvFluidInPlaceOutput>
188 ("Do not print reservoir volumes values after each report step "
189 "even if requested by the deck.");
190 }
191
196 void
197 allocBuffers(const unsigned bufferSize,
198 const unsigned reportStepNum,
199 const bool substep,
200 const bool log,
201 const bool isRestart)
202 {
203 if (! std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value) {
204 return;
205 }
206
207 const auto& problem = this->simulator_.problem();
208
209 this->doAllocBuffers(bufferSize,
210 reportStepNum,
211 substep,
212 log,
213 isRestart,
214 problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)),
215 problem.materialLawManager()->enablePCHysteresis(),
216 problem.materialLawManager()->enableNonWettingHysteresis(),
217 problem.materialLawManager()->enableWettingHysteresis(),
218 problem.tracerModel().numTracers(),
219 problem.tracerModel().enableSolTracers(),
220 problem.eclWriter()->getOutputNnc().size());
221 }
222
223 void processElementMech(const ElementContext& elemCtx)
224 {
226 const auto& problem = elemCtx.simulator().problem();
227 const auto& model = problem.geoMechModel();
228 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
229 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
230 if (!this->mechPotentialForce_.empty()) {
231 // assume all mechanical things should be written
232 this->mechPotentialForce_[globalDofIdx] = model.mechPotentialForce(globalDofIdx);
233 this->mechPotentialPressForce_[globalDofIdx] = model.mechPotentialPressForce(globalDofIdx);
234 this->mechPotentialTempForce_[globalDofIdx] = model.mechPotentialTempForce(globalDofIdx);
235
236 this->dispX_[globalDofIdx] = model.disp(globalDofIdx, 0);
237 this->dispY_[globalDofIdx] = model.disp(globalDofIdx, 1);
238 this->dispZ_[globalDofIdx] = model.disp(globalDofIdx, 2);
239 this->stressXX_[globalDofIdx] = model.stress(globalDofIdx, 0);
240 this->stressYY_[globalDofIdx] = model.stress(globalDofIdx, 1);
241 this->stressZZ_[globalDofIdx] = model.stress(globalDofIdx, 2);
242 // voight notation
243 this->stressXY_[globalDofIdx] = model.stress(globalDofIdx, 5);
244 this->stressXZ_[globalDofIdx] = model.stress(globalDofIdx, 4);
245 this->stressYZ_[globalDofIdx] = model.stress(globalDofIdx, 3);
246
247 this->strainXX_[globalDofIdx] = model.strain(globalDofIdx, 0);
248 this->strainYY_[globalDofIdx] = model.strain(globalDofIdx, 1);
249 this->strainZZ_[globalDofIdx] = model.strain(globalDofIdx, 2);
250 // voight notation
251 this->strainXY_[globalDofIdx] = model.strain(globalDofIdx, 5);
252 this->strainXZ_[globalDofIdx] = model.strain(globalDofIdx, 4);
253 this->strainYZ_[globalDofIdx] = model.strain(globalDofIdx, 3);
254
255
256 this->delstressXX_[globalDofIdx] = model.delstress(globalDofIdx, 0);
257 this->delstressYY_[globalDofIdx] = model.delstress(globalDofIdx, 1);
258 this->delstressZZ_[globalDofIdx] = model.delstress(globalDofIdx, 2);
259 // voight notation
260 this->delstressXY_[globalDofIdx] = model.delstress(globalDofIdx, 5);
261 this->delstressXZ_[globalDofIdx] = model.delstress(globalDofIdx, 4);
262 this->delstressYZ_[globalDofIdx] = model.delstress(globalDofIdx, 3);
263 }
264 }
265 }
266 }
267
272 void processElement(const ElementContext& elemCtx)
273 {
275 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
276 return;
277
278 const auto& problem = elemCtx.simulator().problem();
279 const auto& modelResid = elemCtx.simulator().model().linearizer().residual();
280 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
281 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
282 const auto& fs = intQuants.fluidState();
283
284 using FluidState = std::remove_cv_t<std::remove_reference_t<decltype(fs)>>;
285
286 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
287 const unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex();
288
289 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
290 if (this->saturation_[phaseIdx].empty())
291 continue;
292
293 this->saturation_[phaseIdx][globalDofIdx] = getValue(fs.saturation(phaseIdx));
294 Valgrind::CheckDefined(this->saturation_[phaseIdx][globalDofIdx]);
295 }
296
297 if (this->regionAvgDensity_.has_value()) {
298 // Note: We intentionally exclude effects of rock
299 // compressibility by using referencePorosity() here.
300 const auto porv = intQuants.referencePorosity()
301 * elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
302
303 this->aggregateAverageDensityContributions_(fs, globalDofIdx,
304 static_cast<double>(porv));
305 }
306
307 if (!this->fluidPressure_.empty()) {
308 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
309 // Output oil pressure as default
310 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx));
311 } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
312 // Output gas if oil is not present
313 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx));
314 } else {
315 // Output water if neither oil nor gas is present
316 this->fluidPressure_[globalDofIdx] = getValue(fs.pressure(waterPhaseIdx));
317 }
318 Valgrind::CheckDefined(this->fluidPressure_[globalDofIdx]);
319 }
320
321 if (!this->temperature_.empty()) {
322 this->temperature_[globalDofIdx] = getValue(fs.temperature(oilPhaseIdx));
323 Valgrind::CheckDefined(this->temperature_[globalDofIdx]);
324 }
325 if (!this->gasDissolutionFactor_.empty()) {
326 Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
327 this->gasDissolutionFactor_[globalDofIdx]
329 fs, oilPhaseIdx, pvtRegionIdx, SoMax);
330 Valgrind::CheckDefined(this->gasDissolutionFactor_[globalDofIdx]);
331 }
332 if (!this->oilVaporizationFactor_.empty()) {
333 Scalar SoMax = elemCtx.problem().maxOilSaturation(globalDofIdx);
334 this->oilVaporizationFactor_[globalDofIdx]
336 fs, gasPhaseIdx, pvtRegionIdx, SoMax);
337 Valgrind::CheckDefined(this->oilVaporizationFactor_[globalDofIdx]);
338 }
339 if (!this->gasDissolutionFactorInWater_.empty()) {
340 Scalar SwMax = elemCtx.problem().maxWaterSaturation(globalDofIdx);
341 this->gasDissolutionFactorInWater_[globalDofIdx]
343 fs, waterPhaseIdx, pvtRegionIdx, SwMax);
344 Valgrind::CheckDefined(this->gasDissolutionFactorInWater_[globalDofIdx]);
345 }
346 if (!this->waterVaporizationFactor_.empty()) {
347 this->waterVaporizationFactor_[globalDofIdx]
349 fs, gasPhaseIdx, pvtRegionIdx);
350 Valgrind::CheckDefined(this->waterVaporizationFactor_[globalDofIdx]);
351 }
352 if (!this->gasFormationVolumeFactor_.empty()) {
353 this->gasFormationVolumeFactor_[globalDofIdx] = 1.0
355 fs, gasPhaseIdx, pvtRegionIdx);
356 Valgrind::CheckDefined(this->gasFormationVolumeFactor_[globalDofIdx]);
357 }
358 if (!this->saturatedOilFormationVolumeFactor_.empty()) {
359 this->saturatedOilFormationVolumeFactor_[globalDofIdx] = 1.0
361 fs, oilPhaseIdx, pvtRegionIdx);
362 Valgrind::CheckDefined(this->saturatedOilFormationVolumeFactor_[globalDofIdx]);
363 }
364 if (!this->oilSaturationPressure_.empty()) {
365 this->oilSaturationPressure_[globalDofIdx]
366 = FluidSystem::template saturationPressure<FluidState, Scalar>(fs, oilPhaseIdx, pvtRegionIdx);
367 Valgrind::CheckDefined(this->oilSaturationPressure_[globalDofIdx]);
368 }
369
370 if (!this->rs_.empty()) {
371 this->rs_[globalDofIdx] = getValue(fs.Rs());
372 Valgrind::CheckDefined(this->rs_[globalDofIdx]);
373 }
374 if (!this->rsw_.empty()) {
375 this->rsw_[globalDofIdx] = getValue(fs.Rsw());
376 Valgrind::CheckDefined(this->rsw_[globalDofIdx]);
377 }
378
379 if (!this->rv_.empty()) {
380 this->rv_[globalDofIdx] = getValue(fs.Rv());
381 Valgrind::CheckDefined(this->rv_[globalDofIdx]);
382 }
383 if (!this->pcgw_.empty()) {
384 this->pcgw_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
385 Valgrind::CheckDefined(this->pcgw_[globalDofIdx]);
386 }
387 if (!this->pcow_.empty()) {
388 this->pcow_[globalDofIdx] = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
389 Valgrind::CheckDefined(this->pcow_[globalDofIdx]);
390 }
391 if (!this->pcog_.empty()) {
392 this->pcog_[globalDofIdx] = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
393 Valgrind::CheckDefined(this->pcog_[globalDofIdx]);
394 }
395
396 if (!this->rvw_.empty()) {
397 this->rvw_[globalDofIdx] = getValue(fs.Rvw());
398 Valgrind::CheckDefined(this->rvw_[globalDofIdx]);
399 }
400
401 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
402 if (this->invB_[phaseIdx].empty())
403 continue;
404
405 this->invB_[phaseIdx][globalDofIdx] = getValue(fs.invB(phaseIdx));
406 Valgrind::CheckDefined(this->invB_[phaseIdx][globalDofIdx]);
407 }
408
409 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
410 if (this->density_[phaseIdx].empty())
411 continue;
412
413 this->density_[phaseIdx][globalDofIdx] = getValue(fs.density(phaseIdx));
414 Valgrind::CheckDefined(this->density_[phaseIdx][globalDofIdx]);
415 }
416
417 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
418 if (this->viscosity_[phaseIdx].empty())
419 continue;
420
421 if (!this->extboX_.empty() && phaseIdx == oilPhaseIdx)
422 this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.oilViscosity());
423 else if (!this->extboX_.empty() && phaseIdx == gasPhaseIdx)
424 this->viscosity_[phaseIdx][globalDofIdx] = getValue(intQuants.gasViscosity());
425 else
426 this->viscosity_[phaseIdx][globalDofIdx] = getValue(fs.viscosity(phaseIdx));
427 Valgrind::CheckDefined(this->viscosity_[phaseIdx][globalDofIdx]);
428 }
429
430 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
431 if (this->relativePermeability_[phaseIdx].empty())
432 continue;
433
434 this->relativePermeability_[phaseIdx][globalDofIdx]
435 = getValue(intQuants.relativePermeability(phaseIdx));
436 Valgrind::CheckDefined(this->relativePermeability_[phaseIdx][globalDofIdx]);
437 }
438
439 if (!this->drsdtcon_.empty()) {
440 this->drsdtcon_[globalDofIdx] = problem.drsdtcon(globalDofIdx, elemCtx.simulator().episodeIndex());
441 }
442
443 if (!this->sSol_.empty()) {
444 this->sSol_[globalDofIdx] = intQuants.solventSaturation().value();
445 }
446
447 if (!this->rswSol_.empty()) {
448 this->rswSol_[globalDofIdx] = intQuants.rsSolw().value();
449 }
450
451 if (!this->cPolymer_.empty()) {
452 this->cPolymer_[globalDofIdx] = intQuants.polymerConcentration().value();
453 }
454
455 if (!this->cFoam_.empty()) {
456 this->cFoam_[globalDofIdx] = intQuants.foamConcentration().value();
457 }
458
459 if (!this->cSalt_.empty()) {
460 this->cSalt_[globalDofIdx] = fs.saltConcentration().value();
461 }
462
463 if (!this->pSalt_.empty()) {
464 this->pSalt_[globalDofIdx] = intQuants.saltSaturation().value();
465 }
466
467 if (!this->permFact_.empty()) {
468 this->permFact_[globalDofIdx] = intQuants.permFactor().value();
469 }
470
471 if (!this->extboX_.empty()) {
472 this->extboX_[globalDofIdx] = intQuants.xVolume().value();
473 }
474
475 if (!this->extboY_.empty()) {
476 this->extboY_[globalDofIdx] = intQuants.yVolume().value();
477 }
478
479 if (!this->extboZ_.empty()) {
480 this->extboZ_[globalDofIdx] = intQuants.zFraction().value();
481 }
482
483 if (!this->rPorV_.empty()) {
484 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
485 this->rPorV_[globalDofIdx] = totVolume * intQuants.porosity().value();
486 }
487
488 if (!this->mFracCo2_.empty()) {
489 const Scalar stdVolOil = getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx))
490 + getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx)) * getValue(fs.Rv());
491 const Scalar stdVolGas = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
492 * (1.0 - intQuants.yVolume().value())
493 + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
494 * (1.0 - intQuants.xVolume().value());
495 const Scalar stdVolCo2 = getValue(fs.saturation(gasPhaseIdx)) * getValue(fs.invB(gasPhaseIdx))
496 * intQuants.yVolume().value()
497 + getValue(fs.saturation(oilPhaseIdx)) * getValue(fs.invB(oilPhaseIdx)) * getValue(fs.Rs())
498 * intQuants.xVolume().value();
499 const Scalar rhoO = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
500 const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
501 const Scalar rhoCO2 = intQuants.zRefDensity();
502 const Scalar stdMassTotal = 1.0e-10 + stdVolOil * rhoO + stdVolGas * rhoG + stdVolCo2 * rhoCO2;
503 this->mFracOil_[globalDofIdx] = stdVolOil * rhoO / stdMassTotal;
504 this->mFracGas_[globalDofIdx] = stdVolGas * rhoG / stdMassTotal;
505 this->mFracCo2_[globalDofIdx] = stdVolCo2 * rhoCO2 / stdMassTotal;
506 }
507
508 if (!this->cMicrobes_.empty()) {
509 this->cMicrobes_[globalDofIdx] = intQuants.microbialConcentration().value();
510 }
511
512 if (!this->cOxygen_.empty()) {
513 this->cOxygen_[globalDofIdx] = intQuants.oxygenConcentration().value();
514 }
515
516 if (!this->cUrea_.empty()) {
517 this->cUrea_[globalDofIdx] = 10
518 * intQuants.ureaConcentration()
519 .value(); // Reescaling back the urea concentration (see WellInterface_impl.hpp)
520 }
521
522 if (!this->cBiofilm_.empty()) {
523 this->cBiofilm_[globalDofIdx] = intQuants.biofilmConcentration().value();
524 }
525
526 if (!this->cCalcite_.empty()) {
527 this->cCalcite_[globalDofIdx] = intQuants.calciteConcentration().value();
528 }
529
530 if (!this->bubblePointPressure_.empty()) {
531 try {
532 this->bubblePointPressure_[globalDofIdx]
533 = getValue(FluidSystem::bubblePointPressure(fs, intQuants.pvtRegionIndex()));
534 } catch (const NumericalProblem&) {
535 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
536 this->failedCellsPb_.push_back(cartesianIdx);
537 }
538 }
539
540 if (!this->dewPointPressure_.empty()) {
541 try {
542 this->dewPointPressure_[globalDofIdx]
543 = getValue(FluidSystem::dewPointPressure(fs, intQuants.pvtRegionIndex()));
544 } catch (const NumericalProblem&) {
545 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
546 this->failedCellsPd_.push_back(cartesianIdx);
547 }
548 }
549
550 if (!this->minimumOilPressure_.empty())
551 this->minimumOilPressure_[globalDofIdx]
552 = std::min(getValue(fs.pressure(oilPhaseIdx)), problem.minOilPressure(globalDofIdx));
553
554 if (!this->overburdenPressure_.empty())
555 this->overburdenPressure_[globalDofIdx] = problem.overburdenPressure(globalDofIdx);
556
557 if (!this->rockCompPorvMultiplier_.empty())
558 this->rockCompPorvMultiplier_[globalDofIdx]
560
561 if (!this->rockCompTransMultiplier_.empty())
562 this->rockCompTransMultiplier_[globalDofIdx]
564
565 const auto& matLawManager = problem.materialLawManager();
566 if (matLawManager->enableHysteresis()) {
567 if (FluidSystem::phaseIsActive(oilPhaseIdx)
568 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
569 Scalar somax;
570 Scalar swmax;
571 Scalar swmin;
572
573 matLawManager->oilWaterHysteresisParams(
575
576 if (matLawManager->enableNonWettingHysteresis()) {
577 if (!this->soMax_.empty()) {
578 this->soMax_[globalDofIdx] = somax;
579 }
580 }
581 if (matLawManager->enableWettingHysteresis()) {
582 if (!this->swMax_.empty()) {
583 this->swMax_[globalDofIdx] = swmax;
584 }
585 }
586 if (matLawManager->enablePCHysteresis()) {
587 if (!this->swmin_.empty()) {
588 this->swmin_[globalDofIdx] = swmin;
589 }
590 }
591 }
592
593 if (FluidSystem::phaseIsActive(oilPhaseIdx)
594 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
595 Scalar sgmax;
596 Scalar shmax;
597 Scalar somin;
598 matLawManager->gasOilHysteresisParams(
600
601 if (matLawManager->enableNonWettingHysteresis()) {
602 if (!this->sgmax_.empty()) {
603 this->sgmax_[globalDofIdx] = sgmax;
604 }
605 }
606 if (matLawManager->enableWettingHysteresis()) {
607 if (!this->shmax_.empty()) {
608 this->shmax_[globalDofIdx] = shmax;
609 }
610 }
611 if (matLawManager->enablePCHysteresis()) {
612 if (!this->somin_.empty()) {
613 this->somin_[globalDofIdx] = somin;
614 }
615 }
616 }
617 } else {
618
619 if (!this->soMax_.empty())
620 this->soMax_[globalDofIdx]
621 = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx));
622
623 if (!this->swMax_.empty())
624 this->swMax_[globalDofIdx]
625 = std::max(getValue(fs.saturation(waterPhaseIdx)), problem.maxWaterSaturation(globalDofIdx));
626
627 }
628 if (!this->ppcw_.empty()) {
629 this->ppcw_[globalDofIdx] = matLawManager->oilWaterScaledEpsInfoDrainage(globalDofIdx).maxPcow;
630 // printf("ppcw_[%d] = %lg\n", globalDofIdx, ppcw_[globalDofIdx]);
631 }
632
633 // hack to make the intial output of rs and rv Ecl compatible.
634 // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
635 // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
636 // rs and rv with the values computed in the initially.
637 // Volume factors, densities and viscosities need to be recalculated with the updated rs and rv values.
638 if ((elemCtx.simulator().episodeIndex() < 0) &&
639 FluidSystem::phaseIsActive(oilPhaseIdx) &&
640 FluidSystem::phaseIsActive(gasPhaseIdx))
641 {
642 const auto& fsInitial = problem.initialFluidState(globalDofIdx);
643
644 // use initial rs and rv values
645 if (!this->rv_.empty())
646 this->rv_[globalDofIdx] = fsInitial.Rv();
647
648 if (!this->rs_.empty())
649 this->rs_[globalDofIdx] = fsInitial.Rs();
650
651 if (!this->rsw_.empty())
652 this->rsw_[globalDofIdx] = fsInitial.Rsw();
653
654 if (!this->rvw_.empty())
655 this->rvw_[globalDofIdx] = fsInitial.Rvw();
656
657 // re-compute the volume factors, viscosities and densities if asked for
658 if (!this->density_[oilPhaseIdx].empty())
659 this->density_[oilPhaseIdx][globalDofIdx]
660 = FluidSystem::density(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
661
662 if (!this->density_[gasPhaseIdx].empty())
663 this->density_[gasPhaseIdx][globalDofIdx]
664 = FluidSystem::density(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
665
666 if (!this->invB_[oilPhaseIdx].empty())
667 this->invB_[oilPhaseIdx][globalDofIdx]
668 = FluidSystem::inverseFormationVolumeFactor(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
669
670 if (!this->invB_[gasPhaseIdx].empty())
671 this->invB_[gasPhaseIdx][globalDofIdx]
672 = FluidSystem::inverseFormationVolumeFactor(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
673
674 if (!this->viscosity_[oilPhaseIdx].empty())
675 this->viscosity_[oilPhaseIdx][globalDofIdx]
676 = FluidSystem::viscosity(fsInitial, oilPhaseIdx, intQuants.pvtRegionIndex());
677
678 if (!this->viscosity_[gasPhaseIdx].empty())
679 this->viscosity_[gasPhaseIdx][globalDofIdx]
680 = FluidSystem::viscosity(fsInitial, gasPhaseIdx, intQuants.pvtRegionIndex());
681 }
682
683 // Adding Well RFT data
684 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
685 if (this->oilConnectionPressures_.count(cartesianIdx) > 0) {
686 this->oilConnectionPressures_[cartesianIdx] = getValue(fs.pressure(oilPhaseIdx));
687 }
688 if (this->waterConnectionSaturations_.count(cartesianIdx) > 0) {
689 this->waterConnectionSaturations_[cartesianIdx] = getValue(fs.saturation(waterPhaseIdx));
690 }
691 if (this->gasConnectionSaturations_.count(cartesianIdx) > 0) {
692 this->gasConnectionSaturations_[cartesianIdx] = getValue(fs.saturation(gasPhaseIdx));
693 }
694
695 // tracers
696 const auto& tracerModel = simulator_.problem().tracerModel();
697 if (! this->freeTracerConcentrations_.empty()) {
698 for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
699 if (this->freeTracerConcentrations_[tracerIdx].empty()) {
700 continue;
701 }
702 this->freeTracerConcentrations_[tracerIdx][globalDofIdx] =
703 tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx);
704 }
705 }
706 if (! this->solTracerConcentrations_.empty()) {
707 for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
708 if (this->solTracerConcentrations_[tracerIdx].empty()) {
709 continue;
710 }
711 this->solTracerConcentrations_[tracerIdx][globalDofIdx] =
712 tracerModel.solTracerConcentration(tracerIdx, globalDofIdx);
713
714 }
715 }
716
717 // output residual
718 for ( int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx )
719 {
720 if (!this->residual_[phaseIdx].empty()) {
721 const unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
723 }
724 }
725 }
726 }
727
728 void processElementFlows(const ElementContext& elemCtx)
729 {
730 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
731 if (!std::is_same_v<Discretization, EcfvDiscretization<TypeTag>>)
732 return;
733
734 const auto& problem = elemCtx.simulator().problem();
735 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
736
737 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
738 if (!problem.model().linearizer().getFlowsInfo().empty()) {
739 const auto& flowsInf = problem.model().linearizer().getFlowsInfo();
741 for (auto& flowsInfo : flowsInfos) {
742 if (flowsInfo.faceId >= 0) {
743 if (!this->flows_[flowsInfo.faceId][gasCompIdx].empty()) {
744 this->flows_[flowsInfo.faceId][gasCompIdx][globalDofIdx]
745 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
746 }
747 if (!this->flows_[flowsInfo.faceId][oilCompIdx].empty()) {
748 this->flows_[flowsInfo.faceId][oilCompIdx][globalDofIdx]
749 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
750 }
751 if (!this->flows_[flowsInfo.faceId][waterCompIdx].empty()) {
752 this->flows_[flowsInfo.faceId][waterCompIdx][globalDofIdx]
753 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
754 }
755 }
756 if (flowsInfo.faceId == -2) {
757 if (!this->flowsn_[gasCompIdx].indices.empty()) {
758 this->flowsn_[gasCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
759 this->flowsn_[gasCompIdx].values[flowsInfo.nncId]
760 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
761 }
762 if (!this->flowsn_[oilCompIdx].indices.empty()) {
763 this->flowsn_[oilCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
764 this->flowsn_[oilCompIdx].values[flowsInfo.nncId]
765 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
766 }
767 if (!this->flowsn_[waterCompIdx].indices.empty()) {
768 this->flowsn_[waterCompIdx].indices[flowsInfo.nncId] = flowsInfo.nncId;
769 this->flowsn_[waterCompIdx].values[flowsInfo.nncId]
770 = flowsInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
771 }
772 }
773 }
774 }
775
776 // flores
777 if (!problem.model().linearizer().getFloresInfo().empty()) {
778 const auto& floresInf = problem.model().linearizer().getFloresInfo();
780 for (auto& floresInfo : floresInfos) {
781 if (floresInfo.faceId >= 0) {
782 if (!this->flores_[floresInfo.faceId][gasCompIdx].empty()) {
783 this->flores_[floresInfo.faceId][gasCompIdx][globalDofIdx]
784 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
785 }
786 if (!this->flores_[floresInfo.faceId][oilCompIdx].empty()) {
787 this->flores_[floresInfo.faceId][oilCompIdx][globalDofIdx]
788 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
789 }
790 if (!this->flores_[floresInfo.faceId][waterCompIdx].empty()) {
791 this->flores_[floresInfo.faceId][waterCompIdx][globalDofIdx]
792 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
793 }
794 }
795
796 if (floresInfo.faceId == -2) {
797 if (!this->floresn_[gasCompIdx].indices.empty()) {
798 this->floresn_[gasCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
799 this->floresn_[gasCompIdx].values[floresInfo.nncId]
800 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(gasCompIdx)];
801 }
802 if (!this->floresn_[oilCompIdx].indices.empty()) {
803 this->floresn_[oilCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
804 this->floresn_[oilCompIdx].values[floresInfo.nncId]
805 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(oilCompIdx)];
806 }
807 if (!this->floresn_[waterCompIdx].indices.empty()) {
808 this->floresn_[waterCompIdx].indices[floresInfo.nncId] = floresInfo.nncId;
809 this->floresn_[waterCompIdx].values[floresInfo.nncId]
810 = floresInfo.flow[conti0EqIdx + Indices::canonicalToActiveComponentIndex(waterCompIdx)];
811 }
812 }
813 }
814 }
815 }
816 }
817
818 void processElementBlockData(const ElementContext& elemCtx)
819 {
820 OPM_TIMEBLOCK_LOCAL(processElementBlockData);
821 if (!std::is_same<Discretization, EcfvDiscretization<TypeTag>>::value)
822 return;
823
824 const auto& problem = elemCtx.simulator().problem();
825 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
826 // Adding block data
827 const auto globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
828 const auto cartesianIdx = elemCtx.simulator().vanguard().cartesianIndex(globalDofIdx);
829 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
830 const auto& fs = intQuants.fluidState();
831 for (auto& val : this->blockData_) {
832 const auto& key = val.first;
833 assert(key.second > 0);
834
835 const auto cartesianIdxBlock = static_cast<std::remove_cv_t<
836 std::remove_reference_t<decltype(cartesianIdx)>>>(key.second - 1);
837
839 if ((key.first == "BWSAT") || (key.first == "BSWAT"))
840 val.second = getValue(fs.saturation(waterPhaseIdx));
841 else if ((key.first == "BGSAT") || (key.first == "BSGAS"))
842 val.second = getValue(fs.saturation(gasPhaseIdx));
843 else if ((key.first == "BOSAT") || (key.first == "BSOIL"))
844 val.second = getValue(fs.saturation(oilPhaseIdx));
845 else if (key.first == "BNSAT")
846 val.second = intQuants.solventSaturation().value();
847 else if ((key.first == "BPR") || (key.first == "BPRESSUR")) {
848 if (FluidSystem::phaseIsActive(oilPhaseIdx))
849 val.second = getValue(fs.pressure(oilPhaseIdx));
850 else if (FluidSystem::phaseIsActive(gasPhaseIdx))
851 val.second = getValue(fs.pressure(gasPhaseIdx));
852 else if (FluidSystem::phaseIsActive(waterPhaseIdx))
853 val.second = getValue(fs.pressure(waterPhaseIdx));
854 }
855 else if ((key.first == "BTCNFHEA") || (key.first == "BTEMP")) {
856 if (FluidSystem::phaseIsActive(oilPhaseIdx))
857 val.second = getValue(fs.temperature(oilPhaseIdx));
858 else if (FluidSystem::phaseIsActive(gasPhaseIdx))
859 val.second = getValue(fs.temperature(gasPhaseIdx));
860 else if (FluidSystem::phaseIsActive(waterPhaseIdx))
861 val.second = getValue(fs.temperature(waterPhaseIdx));
862 }
863 else if (key.first == "BWKR" || key.first == "BKRW")
864 val.second = getValue(intQuants.relativePermeability(waterPhaseIdx));
865 else if (key.first == "BGKR" || key.first == "BKRG")
866 val.second = getValue(intQuants.relativePermeability(gasPhaseIdx));
867 else if (key.first == "BOKR" || key.first == "BKRO")
868 val.second = getValue(intQuants.relativePermeability(oilPhaseIdx));
869 else if (key.first == "BKROG") {
870 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0);
871 const auto krog
873 val.second = getValue(krog);
874 }
875 else if (key.first == "BKROW") {
876 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, /* timeIdx = */ 0);
877 const auto krow
879 val.second = getValue(krow);
880 }
881 else if (key.first == "BWPC")
882 val.second = getValue(fs.pressure(oilPhaseIdx)) - getValue(fs.pressure(waterPhaseIdx));
883 else if (key.first == "BGPC")
884 val.second = getValue(fs.pressure(gasPhaseIdx)) - getValue(fs.pressure(oilPhaseIdx));
885 else if (key.first == "BWPR")
886 val.second = getValue(fs.pressure(waterPhaseIdx));
887 else if (key.first == "BGPR")
888 val.second = getValue(fs.pressure(gasPhaseIdx));
889 else if (key.first == "BVWAT" || key.first == "BWVIS")
890 val.second = getValue(fs.viscosity(waterPhaseIdx));
891 else if (key.first == "BVGAS" || key.first == "BGVIS")
892 val.second = getValue(fs.viscosity(gasPhaseIdx));
893 else if (key.first == "BVOIL" || key.first == "BOVIS")
894 val.second = getValue(fs.viscosity(oilPhaseIdx));
895 else if ((key.first == "BODEN") || (key.first == "BDENO"))
896 val.second = getValue(fs.density(oilPhaseIdx));
897 else if ((key.first == "BGDEN") || (key.first == "BDENG"))
898 val.second = getValue(fs.density(gasPhaseIdx));
899 else if ((key.first == "BWDEN") || (key.first == "BDENW"))
900 val.second = getValue(fs.density(waterPhaseIdx));
901 else if ((key.first == "BRPV") ||
902 (key.first == "BOPV") ||
903 (key.first == "BWPV") ||
904 (key.first == "BGPV"))
905 {
906 if (key.first == "BRPV") {
907 val.second = 1.0;
908 }
909 else if (key.first == "BOPV") {
910 val.second = getValue(fs.saturation(oilPhaseIdx));
911 }
912 else if (key.first == "BWPV") {
913 val.second = getValue(fs.saturation(waterPhaseIdx));
914 }
915 else {
916 val.second = getValue(fs.saturation(gasPhaseIdx));
917 }
918
919 // Include active pore-volume.
920 val.second *= getValue(intQuants.porosity())
921 * elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
922 }
923 else if (key.first == "BRS")
924 val.second = getValue(fs.Rs());
925 else if (key.first == "BRV")
926 val.second = getValue(fs.Rv());
927 else if ((key.first == "BOIP") || (key.first == "BOIPL") || (key.first == "BOIPG") ||
928 (key.first == "BGIP") || (key.first == "BGIPL") || (key.first == "BGIPG") ||
929 (key.first == "BWIP"))
930 {
931 if ((key.first == "BOIP") || (key.first == "BOIPL")) {
932 val.second = getValue(fs.invB(oilPhaseIdx)) * getValue(fs.saturation(oilPhaseIdx));
933
934 if (key.first == "BOIP") {
935 val.second += getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
936 * getValue(fs.saturation(gasPhaseIdx));
937 }
938 }
939 else if (key.first == "BOIPG") {
940 val.second = getValue(fs.Rv()) * getValue(fs.invB(gasPhaseIdx))
941 * getValue(fs.saturation(gasPhaseIdx));
942 }
943 else if ((key.first == "BGIP") || (key.first == "BGIPG")) {
944 val.second = getValue(fs.invB(gasPhaseIdx)) * getValue(fs.saturation(gasPhaseIdx));
945
946 if (key.first == "BGIP") {
947 if (!FluidSystem::phaseIsActive(oilPhaseIdx)) {
948 val.second += getValue(fs.Rsw()) * getValue(fs.invB(waterPhaseIdx))
949 * getValue(fs.saturation(waterPhaseIdx));
950 }
951 else {
952 val.second += getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
953 * getValue(fs.saturation(oilPhaseIdx));
954 }
955 }
956 }
957 else if (key.first == "BGIPL") {
958 if (!FluidSystem::phaseIsActive(oilPhaseIdx)) {
959 val.second = getValue(fs.Rsw()) * getValue(fs.invB(waterPhaseIdx))
960 * getValue(fs.saturation(waterPhaseIdx));
961 }
962 else {
963 val.second = getValue(fs.Rs()) * getValue(fs.invB(oilPhaseIdx))
964 * getValue(fs.saturation(oilPhaseIdx));
965 }
966 }
967 else { // BWIP
968 val.second = getValue(fs.invB(waterPhaseIdx)) * getValue(fs.saturation(waterPhaseIdx));
969 }
970
971 // Include active pore-volume.
972 val.second *= elemCtx.simulator().model().dofTotalVolume(globalDofIdx)
973 * getValue(intQuants.porosity());
974 }
975 else if ((key.first == "BPPO") ||
976 (key.first == "BPPG") ||
977 (key.first == "BPPW"))
978 {
979 auto phase = RegionPhasePoreVolAverage::Phase{};
980
981 if (key.first == "BPPO") {
982 phase.ix = oilPhaseIdx;
983 }
984 else if (key.first == "BPPG") {
985 phase.ix = gasPhaseIdx;
986 }
987 else { // BPPW
988 phase.ix = waterPhaseIdx;
989 }
990
991 // Note different region handling here. FIPNUM is
992 // one-based, but we need zero-based lookup in
993 // DatumDepth. On the other hand, pvtRegionIndex is
994 // zero-based but we need one-based lookup in
995 // RegionPhasePoreVolAverage.
996
997 // Subtract one to convert FIPNUM to region index.
998 const auto datum = this->eclState_.getSimulationConfig()
999 .datumDepths()(this->regions_["FIPNUM"][dofIdx] - 1);
1000
1001 // Add one to convert region index to region ID.
1002 const auto region = RegionPhasePoreVolAverage::Region {
1003 elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex() + 1
1004 };
1005
1006 const auto density = this->regionAvgDensity_
1007 ->value("PVTNUM", phase, region);
1008
1009 const auto press = getValue(fs.pressure(phase.ix));
1010 const auto grav =
1011 elemCtx.problem().gravity()[GridView::dimensionworld - 1];
1012 const auto dz = problem.dofCenterDepth(globalDofIdx) - datum;
1013
1014 val.second = press - density*dz*grav;
1015 }
1016 else if ((key.first == "BFLOWI") ||
1017 (key.first == "BFLOWJ") ||
1018 (key.first == "BFLOWK"))
1019 {
1020 auto dir = FaceDir::ToIntersectionIndex(Dir::XPlus);
1021
1022 if (key.first == "BFLOWJ") {
1023 dir = FaceDir::ToIntersectionIndex(Dir::YPlus);
1024 }
1025 else if (key.first == "BFLOWK") {
1026 dir = FaceDir::ToIntersectionIndex(Dir::ZPlus);
1027 }
1028
1029 val.second = this->flows_[dir][waterCompIdx][globalDofIdx];
1030 }
1031 else {
1032 std::string logstring = "Keyword '";
1033 logstring.append(key.first);
1034 logstring.append("' is unhandled for output to summary file.");
1035 OpmLog::warning("Unhandled output keyword", logstring);
1036 }
1037 }
1038 }
1039 }
1040 }
1041
1070 template <class ActiveIndex, class CartesianIndex>
1071 void processFluxes(const ElementContext& elemCtx,
1072 ActiveIndex&& activeIndex,
1073 CartesianIndex&& cartesianIndex)
1074 {
1076 const auto identifyCell = [&activeIndex, &cartesianIndex](const Element& elem)
1078 {
1079 const auto cellIndex = activeIndex(elem);
1080
1081 return {
1082 static_cast<int>(cellIndex),
1083 cartesianIndex(cellIndex),
1084 elem.partitionType() == Dune::InteriorEntity
1085 };
1086 };
1087
1088 const auto timeIdx = 0u;
1089 const auto& stencil = elemCtx.stencil(timeIdx);
1090 const auto numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
1091
1092 for (auto scvfIdx = 0 * numInteriorFaces; scvfIdx < numInteriorFaces; ++scvfIdx) {
1093 const auto& face = stencil.interiorFace(scvfIdx);
1094 const auto left = identifyCell(stencil.element(face.interiorIndex()));
1095 const auto right = identifyCell(stencil.element(face.exteriorIndex()));
1096
1097 const auto rates = this->
1098 getComponentSurfaceRates(elemCtx, face.area(), scvfIdx, timeIdx);
1099
1100 this->interRegionFlows_.addConnection(left, right, rates);
1101 }
1102 }
1103
1109 {
1110 // Inter-region flow rates. Note: ".clear()" prepares to accumulate
1111 // contributions per bulk connection between FIP regions.
1112 this->interRegionFlows_.clear();
1113 }
1114
1119 {
1120 this->interRegionFlows_.compress();
1121 }
1122
1127 {
1128 return this->interRegionFlows_;
1129 }
1130
1131 template <class FluidState>
1132 void assignToFluidState(FluidState& fs, unsigned elemIdx) const
1133 {
1134 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1135 if (this->saturation_[phaseIdx].empty())
1136 continue;
1137
1138 fs.setSaturation(phaseIdx, this->saturation_[phaseIdx][elemIdx]);
1139 }
1140
1141 if (!this->fluidPressure_.empty()) {
1142 // this assumes that capillary pressures only depend on the phase saturations
1143 // and possibly on temperature. (this is always the case for ECL problems.)
1144 std::array<Scalar, numPhases> pc = {0};
1145 const MaterialLawParams& matParams = simulator_.problem().materialLawParams(elemIdx);
1146 MaterialLaw::capillaryPressures(pc, matParams, fs);
1147 Valgrind::CheckDefined(this->fluidPressure_[elemIdx]);
1148 Valgrind::CheckDefined(pc);
1149 const auto& pressure = this->fluidPressure_[elemIdx];
1150 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
1151 if (!FluidSystem::phaseIsActive(phaseIdx))
1152 continue;
1153
1154 if (Indices::oilEnabled)
1155 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
1156 else if (Indices::gasEnabled)
1157 fs.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
1158 else if (Indices::waterEnabled)
1159 //single (water) phase
1160 fs.setPressure(phaseIdx, pressure);
1161 }
1162 }
1163
1164 if (!this->temperature_.empty())
1165 fs.setTemperature(this->temperature_[elemIdx]);
1166 if (!this->rs_.empty())
1167 fs.setRs(this->rs_[elemIdx]);
1168 if (!this->rsw_.empty())
1169 fs.setRsw(this->rsw_[elemIdx]);
1170 if (!this->rv_.empty())
1171 fs.setRv(this->rv_[elemIdx]);
1172 if (!this->rvw_.empty())
1173 fs.setRvw(this->rvw_[elemIdx]);
1174 }
1175
1176 void initHysteresisParams(Simulator& simulator, unsigned elemIdx) const
1177 {
1178 if (!this->soMax_.empty())
1179 simulator.problem().setMaxOilSaturation(elemIdx, this->soMax_[elemIdx]);
1180
1181 if (simulator.problem().materialLawManager()->enableHysteresis()) {
1182 auto matLawManager = simulator.problem().materialLawManager();
1183
1184 if (FluidSystem::phaseIsActive(oilPhaseIdx)
1185 && FluidSystem::phaseIsActive(waterPhaseIdx)) {
1186 Scalar somax = 2.0;
1187 Scalar swmax = -2.0;
1188 Scalar swmin = 2.0;
1189
1190 if (matLawManager->enableNonWettingHysteresis()) {
1191 if (!this->soMax_.empty()) {
1192 somax = this->soMax_[elemIdx];
1193 }
1194 }
1195 if (matLawManager->enableWettingHysteresis()) {
1196 if (!this->swMax_.empty()) {
1197 swmax = this->swMax_[elemIdx];
1198 }
1199 }
1200 if (matLawManager->enablePCHysteresis()) {
1201 if (!this->swmin_.empty()) {
1202 swmin = this->swmin_[elemIdx];
1203 }
1204 }
1205 matLawManager->setOilWaterHysteresisParams(
1207 }
1208 if (FluidSystem::phaseIsActive(oilPhaseIdx)
1209 && FluidSystem::phaseIsActive(gasPhaseIdx)) {
1210 Scalar sgmax = 2.0;
1211 Scalar shmax = -2.0;
1212 Scalar somin = 2.0;
1213
1214 if (matLawManager->enableNonWettingHysteresis()) {
1215 if (!this->sgmax_.empty()) {
1216 sgmax = this->sgmax_[elemIdx];
1217 }
1218 }
1219 if (matLawManager->enableWettingHysteresis()) {
1220 if (!this->shmax_.empty()) {
1221 shmax = this->shmax_[elemIdx];
1222 }
1223 }
1224 if (matLawManager->enablePCHysteresis()) {
1225 if (!this->somin_.empty()) {
1226 somin = this->somin_[elemIdx];
1227 }
1228 }
1229 matLawManager->setGasOilHysteresisParams(
1231 }
1232
1233 }
1234
1235 if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) {
1236 simulator.problem().materialLawManager()
1237 ->applyRestartSwatInit(elemIdx, this->ppcw_[elemIdx]);
1238 }
1239 }
1240
1241 void updateFluidInPlace(const ElementContext& elemCtx)
1242 {
1243 for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
1244 updateFluidInPlace_(elemCtx, dofIdx);
1245 }
1246 }
1247
1248 void updateFluidInPlace(const unsigned globalDofIdx,
1249 const IntensiveQuantities& intQuants,
1250 const double totVolume)
1251 {
1252 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
1253 }
1254
1255private:
1256 bool isDefunctParallelWell(std::string wname) const override
1257 {
1258 if (simulator_.gridView().comm().size() == 1)
1259 return false;
1260 const auto& parallelWells = simulator_.vanguard().parallelWells();
1261 std::pair<std::string, bool> value {wname, true};
1262 auto candidate = std::lower_bound(parallelWells.begin(), parallelWells.end(), value);
1263 return candidate == parallelWells.end() || *candidate != value;
1264 }
1265
1266 void updateFluidInPlace_(const ElementContext& elemCtx, const unsigned dofIdx)
1267 {
1268 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
1269 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
1270 const auto totVolume = elemCtx.simulator().model().dofTotalVolume(globalDofIdx);
1271
1272 this->updateFluidInPlace_(globalDofIdx, intQuants, totVolume);
1273 }
1274
1275 void updateFluidInPlace_(const unsigned globalDofIdx,
1276 const IntensiveQuantities& intQuants,
1277 const double totVolume)
1278 {
1279 OPM_TIMEBLOCK_LOCAL(updateFluidInPlace);
1280
1281 this->updateTotalVolumesAndPressures_(globalDofIdx, intQuants, totVolume);
1282
1283 if (this->computeFip_) {
1284 this->updatePhaseInplaceVolumes_(globalDofIdx, intQuants, totVolume);
1285 }
1286 }
1287
1288 void createLocalRegion_(std::vector<int>& region)
1289 {
1290 std::size_t elemIdx = 0;
1291 for (const auto& elem : elements(simulator_.gridView())) {
1292 if (elem.partitionType() != Dune::InteriorEntity) {
1293 region[elemIdx] = 0;
1294 }
1295
1296 ++elemIdx;
1297 }
1298 }
1299
1300 template <typename FluidState>
1301 void aggregateAverageDensityContributions_(const FluidState& fs,
1302 const unsigned int globalDofIdx,
1303 const double porv)
1304 {
1305 auto pvCellValue = RegionPhasePoreVolAverage::CellValue{};
1306 pvCellValue.porv = porv;
1307
1308 for (auto phaseIdx = 0*FluidSystem::numPhases;
1309 phaseIdx < FluidSystem::numPhases; ++phaseIdx)
1310 {
1311 if (! FluidSystem::phaseIsActive(phaseIdx)) {
1312 continue;
1313 }
1314
1315 pvCellValue.value = getValue(fs.density(phaseIdx));
1316 pvCellValue.sat = getValue(fs.saturation(phaseIdx));
1317
1318 this->regionAvgDensity_
1319 ->addCell(globalDofIdx,
1320 RegionPhasePoreVolAverage::Phase { phaseIdx },
1321 pvCellValue);
1322 }
1323 }
1324
1340 data::InterRegFlowMap::FlowRates
1341 getComponentSurfaceRates(const ElementContext& elemCtx,
1342 const Scalar faceArea,
1343 const std::size_t scvfIdx,
1344 const std::size_t timeIdx) const
1345 {
1346 using Component = data::InterRegFlowMap::Component;
1347
1348 auto rates = data::InterRegFlowMap::FlowRates {};
1349
1350 const auto& extQuant = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
1351
1352 const auto alpha = getValue(extQuant.extrusionFactor()) * faceArea;
1353
1354 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1355 const auto& up = elemCtx
1356 .intensiveQuantities(extQuant.upstreamIndex(oilPhaseIdx), timeIdx);
1357
1358 using FluidState = std::remove_cv_t<std::remove_reference_t<
1359 decltype(up.fluidState())>>;
1360
1361 const auto pvtReg = up.pvtRegionIndex();
1362
1364 (up.fluidState(), oilPhaseIdx, pvtReg));
1365
1366 const auto qO = alpha * bO * getValue(extQuant.volumeFlux(oilPhaseIdx));
1367
1368 rates[Component::Oil] += qO;
1369
1370 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1371 const auto Rs = getValue(
1372 BlackOil::getRs_<FluidSystem, FluidState, Scalar>
1373 (up.fluidState(), pvtReg));
1374
1375 rates[Component::Gas] += qO * Rs;
1376 rates[Component::Disgas] += qO * Rs;
1377 }
1378 }
1379
1380 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1381 const auto& up = elemCtx
1382 .intensiveQuantities(extQuant.upstreamIndex(gasPhaseIdx), timeIdx);
1383
1384 using FluidState = std::remove_cv_t<std::remove_reference_t<
1385 decltype(up.fluidState())>>;
1386
1387 const auto pvtReg = up.pvtRegionIndex();
1388
1390 (up.fluidState(), gasPhaseIdx, pvtReg));
1391
1392 const auto qG = alpha * bG * getValue(extQuant.volumeFlux(gasPhaseIdx));
1393
1394 rates[Component::Gas] += qG;
1395
1396 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1397 const auto Rv = getValue(
1398 BlackOil::getRv_<FluidSystem, FluidState, Scalar>
1399 (up.fluidState(), pvtReg));
1400
1401 rates[Component::Oil] += qG * Rv;
1402 rates[Component::Vapoil] += qG * Rv;
1403 }
1404 }
1405
1406 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
1407 const auto& up = elemCtx
1408 .intensiveQuantities(extQuant.upstreamIndex(waterPhaseIdx), timeIdx);
1409
1410 using FluidState = std::remove_cv_t<std::remove_reference_t<
1411 decltype(up.fluidState())>>;
1412
1413 const auto pvtReg = up.pvtRegionIndex();
1414
1416 (up.fluidState(), waterPhaseIdx, pvtReg));
1417
1418 rates[Component::Water] +=
1419 alpha * bW * getValue(extQuant.volumeFlux(waterPhaseIdx));
1420 }
1421
1422 return rates;
1423 }
1424
1425 template <typename FluidState>
1426 Scalar hydroCarbonFraction(const FluidState& fs) const
1427 {
1428 if (this->eclState_.runspec().co2Storage()) {
1429 // CO2 storage: Hydrocarbon volume is full pore-volume.
1430 return 1.0;
1431 }
1432
1433 // Common case. Hydrocarbon volume is fraction occupied by actual
1434 // hydrocarbons.
1435 auto hydrocarbon = Scalar {0};
1436 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1437 hydrocarbon += getValue(fs.saturation(oilPhaseIdx));
1438 }
1439
1440 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1441 hydrocarbon += getValue(fs.saturation(gasPhaseIdx));
1442 }
1443
1444 return hydrocarbon;
1445 }
1446
1447 void updateTotalVolumesAndPressures_(const unsigned globalDofIdx,
1448 const IntensiveQuantities& intQuants,
1449 const double totVolume)
1450 {
1451 const auto& fs = intQuants.fluidState();
1452
1453 const double pv = totVolume * intQuants.porosity().value();
1454 const auto hydrocarbon = this->hydroCarbonFraction(fs);
1455
1456 if (! this->hydrocarbonPoreVolume_.empty()) {
1457 this->fip_[Inplace::Phase::PoreVolume][globalDofIdx] =
1458 totVolume * intQuants.referencePorosity();
1459
1460 this->dynamicPoreVolume_[globalDofIdx] = pv;
1461 this->hydrocarbonPoreVolume_[globalDofIdx] = pv * hydrocarbon;
1462 }
1463
1464 if (!this->pressureTimesHydrocarbonVolume_.empty() &&
1465 !this->pressureTimesPoreVolume_.empty())
1466 {
1467 assert(this->hydrocarbonPoreVolume_.size() == this->pressureTimesHydrocarbonVolume_.size());
1468 assert(this->fip_[Inplace::Phase::PoreVolume].size() == this->pressureTimesPoreVolume_.size());
1469
1470 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
1471 this->pressureTimesPoreVolume_[globalDofIdx] =
1472 getValue(fs.pressure(oilPhaseIdx)) * pv;
1473
1474 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
1475 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
1476 }
1477 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
1478 this->pressureTimesPoreVolume_[globalDofIdx] =
1479 getValue(fs.pressure(gasPhaseIdx)) * pv;
1480
1481 this->pressureTimesHydrocarbonVolume_[globalDofIdx] =
1482 this->pressureTimesPoreVolume_[globalDofIdx] * hydrocarbon;
1483 }
1484 else if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
1485 this->pressureTimesPoreVolume_[globalDofIdx] =
1486 getValue(fs.pressure(waterPhaseIdx)) * pv;
1487 }
1488 }
1489 }
1490
1491 void updatePhaseInplaceVolumes_(const unsigned globalDofIdx,
1492 const IntensiveQuantities& intQuants,
1493 const double totVolume)
1494 {
1495 std::array<Scalar, FluidSystem::numPhases> fip {};
1496 std::array<Scalar, FluidSystem::numPhases> fipr{}; // at reservoir condition
1497
1498 const auto& fs = intQuants.fluidState();
1499 const auto pv = totVolume * intQuants.porosity().value();
1500
1501 for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) {
1502 if (!FluidSystem::phaseIsActive(phaseIdx)) {
1503 continue;
1504 }
1505
1506 const auto b = getValue(fs.invB(phaseIdx));
1507 const auto s = getValue(fs.saturation(phaseIdx));
1508
1509 fipr[phaseIdx] = s * pv;
1510 fip [phaseIdx] = b * fipr[phaseIdx];
1511 }
1512
1513 this->updateInplaceVolumesSurface(globalDofIdx, fip);
1514 this->updateInplaceVolumesReservoir(globalDofIdx, fs, fipr);
1515
1516 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1517 FluidSystem::phaseIsActive(gasPhaseIdx))
1518 {
1519 this->updateOilGasDistribution(globalDofIdx, fs, fip);
1520 }
1521
1522 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1523 FluidSystem::phaseIsActive(gasPhaseIdx))
1524 {
1525 this->updateGasWaterDistribution(globalDofIdx, fs, fip);
1526 }
1527
1528 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1529 (!this->fip_[Inplace::Phase::CO2InGasPhaseInMob].empty() ||
1530 !this->fip_[Inplace::Phase::CO2InGasPhaseMob].empty() ||
1531 !this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob].empty() ||
1532 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMob].empty() ||
1533 !this->fip_[Inplace::Phase::CO2Mass].empty() ||
1534 !this->fip_[Inplace::Phase::CO2MassInGasPhase].empty() ||
1535 !this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg].empty() ||
1536 !this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg].empty() ||
1537 !this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg].empty() ||
1538 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg].empty() ||
1539 !this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped].empty() ||
1540 !this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped].empty() ||
1541 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped].empty() ||
1542 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped].empty()))
1543 {
1544 this->updateCO2InGas(globalDofIdx, pv, intQuants);
1545 }
1546
1547 if ((!this->fip_[Inplace::Phase::CO2InWaterPhase].empty() ||
1548 !this->fip_[Inplace::Phase::CO2MassInWaterPhase].empty() ||
1549 !this->fip_[Inplace::Phase::CO2Mass].empty()) &&
1550 (FluidSystem::phaseIsActive(waterPhaseIdx) ||
1551 FluidSystem::phaseIsActive(oilPhaseIdx)))
1552 {
1553 this->updateCO2InWater(globalDofIdx, pv, fs);
1554 }
1555 }
1556
1557 template <typename FIPArray>
1558 void updateInplaceVolumesSurface(const unsigned globalDofIdx,
1559 const FIPArray& fip)
1560 {
1561 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1562 !this->fip_[Inplace::Phase::OIL].empty())
1563 {
1564 this->fip_[Inplace::Phase::OIL][globalDofIdx] = fip[oilPhaseIdx];
1565 }
1566
1567 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1568 !this->fip_[Inplace::Phase::OilInLiquidPhase].empty())
1569 {
1570 this->fip_[Inplace::Phase::OilInLiquidPhase][globalDofIdx] = fip[oilPhaseIdx];
1571 }
1572
1573 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1574 !this->fip_[Inplace::Phase::GAS].empty())
1575 {
1576 this->fip_[Inplace::Phase::GAS][globalDofIdx] = fip[gasPhaseIdx];
1577 }
1578
1579 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1580 !this->fip_[Inplace::Phase::GasInGasPhase].empty())
1581 {
1582 this->fip_[Inplace::Phase::GasInGasPhase][globalDofIdx] = fip[gasPhaseIdx];
1583 }
1584
1585 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1586 !this->fip_[Inplace::Phase::WATER].empty())
1587 {
1588 this->fip_[Inplace::Phase::WATER][globalDofIdx] = fip[waterPhaseIdx];
1589 }
1590 }
1591
1592 template <typename FluidState, typename FIPArray>
1593 void updateInplaceVolumesReservoir(const unsigned globalDofIdx,
1594 const FluidState& fs,
1595 const FIPArray& fipr)
1596 {
1597 if (FluidSystem::phaseIsActive(oilPhaseIdx) &&
1598 !this->fip_[Inplace::Phase::OilResVolume].empty())
1599 {
1600 this->fip_[Inplace::Phase::OilResVolume][globalDofIdx] = fipr[oilPhaseIdx];
1601 }
1602
1603 if (FluidSystem::phaseIsActive(gasPhaseIdx) &&
1604 !this->fip_[Inplace::Phase::GasResVolume].empty())
1605 {
1606 this->fip_[Inplace::Phase::GasResVolume][globalDofIdx] = fipr[gasPhaseIdx];
1607 }
1608
1609 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1610 !this->fip_[Inplace::Phase::WaterResVolume].empty())
1611 {
1612 this->fip_[Inplace::Phase::WaterResVolume][globalDofIdx] = fipr[waterPhaseIdx];
1613 }
1614
1615 if (FluidSystem::phaseIsActive(waterPhaseIdx) &&
1616 !this->fip_[Inplace::Phase::SALT].empty())
1617 {
1618 this->fip_[Inplace::Phase::SALT][globalDofIdx] =
1619 fipr[waterPhaseIdx] * fs.saltConcentration().value();
1620 }
1621 }
1622
1623 template <typename FluidState, typename FIPArray>
1624 void updateOilGasDistribution(const unsigned globalDofIdx,
1625 const FluidState& fs,
1626 const FIPArray& fip)
1627 {
1628 // Gas dissolved in oil and vaporized oil
1629 const auto gasInPlaceLiquid = getValue(fs.Rs()) * fip[oilPhaseIdx];
1630 const auto oilInPlaceGas = getValue(fs.Rv()) * fip[gasPhaseIdx];
1631
1632 if (!this->fip_[Inplace::Phase::GasInLiquidPhase].empty()) {
1633 this->fip_[Inplace::Phase::GasInLiquidPhase][globalDofIdx] = gasInPlaceLiquid;
1634 }
1635
1636 if (!this->fip_[Inplace::Phase::OilInGasPhase].empty()) {
1637 this->fip_[Inplace::Phase::OilInGasPhase][globalDofIdx] = oilInPlaceGas;
1638 }
1639
1640 // Add dissolved gas and vaporized oil to total Fip
1641 if (!this->fip_[Inplace::Phase::OIL].empty()) {
1642 this->fip_[Inplace::Phase::OIL][globalDofIdx] += oilInPlaceGas;
1643 }
1644
1645 if (!this->fip_[Inplace::Phase::GAS].empty()) {
1646 this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceLiquid;
1647 }
1648 }
1649
1650 template <typename FluidState, typename FIPArray>
1651 void updateGasWaterDistribution(const unsigned globalDofIdx,
1652 const FluidState& fs,
1653 const FIPArray& fip)
1654 {
1655 // Gas dissolved in water and vaporized water
1656 const auto gasInPlaceWater = getValue(fs.Rsw()) * fip[waterPhaseIdx];
1657 const auto waterInPlaceGas = getValue(fs.Rvw()) * fip[gasPhaseIdx];
1658
1659 if (!this->fip_[Inplace::Phase::WaterInGasPhase].empty()) {
1660 this->fip_[Inplace::Phase::WaterInGasPhase][globalDofIdx] = waterInPlaceGas;
1661 }
1662
1663 if (!this->fip_[Inplace::Phase::WaterInWaterPhase].empty()) {
1664 this->fip_[Inplace::Phase::WaterInWaterPhase][globalDofIdx] = fip[waterPhaseIdx];
1665 }
1666
1667 // For water+gas cases the gas in water is added to the GIPL value
1668 if (!this->fip_[Inplace::Phase::GasInLiquidPhase].empty() &&
1669 !FluidSystem::phaseIsActive(oilPhaseIdx))
1670 {
1671 this->fip_[Inplace::Phase::GasInLiquidPhase][globalDofIdx] = gasInPlaceWater;
1672 }
1673
1674 // Add dissolved gas and vaporized water to total Fip
1675 if (!this->fip_[Inplace::Phase::WATER].empty()) {
1676 this->fip_[Inplace::Phase::WATER][globalDofIdx] += waterInPlaceGas;
1677 }
1678
1679 if (!this->fip_[Inplace::Phase::GAS].empty()) {
1680 this->fip_[Inplace::Phase::GAS][globalDofIdx] += gasInPlaceWater;
1681 }
1682 }
1683
1684 template <typename IntensiveQuantities>
1685 void updateCO2InGas(const unsigned globalDofIdx,
1686 const double pv,
1687 const IntensiveQuantities& intQuants)
1688 {
1689 const auto& scaledDrainageInfo = this->simulator_.problem().materialLawManager()
1690 ->oilWaterScaledEpsInfoDrainage(globalDofIdx);
1691
1692 const auto& fs = intQuants.fluidState();
1693 Scalar sgcr = scaledDrainageInfo.Sgcr;
1694 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1695 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1696 sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/false);
1697 }
1698
1699 const Scalar sg = getValue(fs.saturation(gasPhaseIdx));
1700 const Scalar rhog = getValue(fs.density(gasPhaseIdx));
1701 const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx)
1702 ? FluidSystem::convertRvwToXgW(getValue(fs.Rvw()), fs.pvtRegionIndex())
1703 : FluidSystem::convertRvToXgO(getValue(fs.Rv()), fs.pvtRegionIndex());
1704
1705 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1706 const Scalar massGas = (1 - xgW) * pv * rhog;
1707 if (!this->fip_[Inplace::Phase::CO2Mass].empty()) {
1708 this->fip_[Inplace::Phase::CO2Mass][globalDofIdx] = massGas * sg;
1709 }
1710
1711 if (!this->fip_[Inplace::Phase::CO2MassInGasPhase].empty()) {
1712 this->fip_[Inplace::Phase::CO2MassInGasPhase][globalDofIdx] = massGas * sg;
1713 }
1714
1715 if (!this->fip_[Inplace::Phase::CO2InGasPhaseInMob].empty()) {
1716 const Scalar imMobileGas = massGas / mM * std::min(sgcr , sg);
1717 this->fip_[Inplace::Phase::CO2InGasPhaseInMob][globalDofIdx] = imMobileGas;
1718 }
1719
1720 if (!this->fip_[Inplace::Phase::CO2InGasPhaseMob].empty()) {
1721 const Scalar mobileGas = massGas / mM * std::max(Scalar{0.0}, sg - sgcr);
1722 this->fip_[Inplace::Phase::CO2InGasPhaseMob][globalDofIdx] = mobileGas;
1723 }
1724
1725 if (!this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg].empty()) {
1726 if (sgcr >= sg) {
1727 const Scalar imMobileGasKrg = massGas / mM * sg;
1728 this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg][globalDofIdx] = imMobileGasKrg;
1729 } else {
1730 this->fip_[Inplace::Phase::CO2InGasPhaseInMobKrg][globalDofIdx] = 0;
1731 }
1732 }
1733
1734 if (!this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg].empty()) {
1735 if (sg > sgcr) {
1736 const Scalar mobileGasKrg = massGas / mM * sg;
1737 this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg][globalDofIdx] = mobileGasKrg;
1738 } else {
1739 this->fip_[Inplace::Phase::CO2InGasPhaseMobKrg][globalDofIdx] = 0;
1740 }
1741 }
1742
1743 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob].empty()) {
1744 const Scalar imMobileMassGas = massGas * std::min(sgcr , sg);
1745 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMob][globalDofIdx] = imMobileMassGas;
1746 }
1747
1748 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMob].empty()) {
1749 const Scalar mobileMassGas = massGas * std::max(Scalar{0.0}, sg - sgcr);
1750 this->fip_[Inplace::Phase::CO2MassInGasPhaseMob][globalDofIdx] = mobileMassGas;
1751 }
1752
1753 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg].empty()) {
1754 if (sgcr >= sg) {
1755 const Scalar imMobileMassGasKrg = massGas * sg;
1756 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg][globalDofIdx] = imMobileMassGasKrg;
1757 } else {
1758 this->fip_[Inplace::Phase::CO2MassInGasPhaseInMobKrg][globalDofIdx] = 0;
1759 }
1760 }
1761
1762 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg].empty()) {
1763 if (sg > sgcr) {
1764 const Scalar mobileMassGasKrg = massGas * sg;
1765 this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg][globalDofIdx] = mobileMassGasKrg;
1766 } else {
1767 this->fip_[Inplace::Phase::CO2MassInGasPhaseMobKrg][globalDofIdx] = 0;
1768 }
1769 }
1770
1771 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped].empty() ||
1772 !this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped].empty() ) {
1774 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1775 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1776 // Get the maximum trapped gas saturation
1777 trappedGasSaturation = MaterialLaw::trappedGasSaturation(matParams, /*maximumTrapping*/true);
1778 }
1779 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped].empty()) {
1780 const Scalar imMobileMassGas = massGas * std::min(trappedGasSaturation , sg);
1781 this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumTrapped][globalDofIdx] = imMobileMassGas;
1782 }
1783 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped].empty()) {
1784 const Scalar mobileMassGas = massGas * std::max(Scalar{0.0}, sg - trappedGasSaturation);
1785 this->fip_[Inplace::Phase::CO2MassInGasPhaseMaximumUnTrapped][globalDofIdx] = mobileMassGas;
1786 }
1787 }
1788
1789 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped].empty() ||
1790 !this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped].empty()) {
1792 if (this->simulator_.problem().materialLawManager()->enableHysteresis()) {
1793 const auto& matParams = simulator_.problem().materialLawParams(globalDofIdx);
1794 const double krg = getValue(intQuants.relativePermeability(gasPhaseIdx));
1795 trappedGasSaturation = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
1796 }
1797 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped].empty()) {
1798 const Scalar imMobileMassGas = massGas * std::min(trappedGasSaturation , sg);
1799 this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveTrapped][globalDofIdx] = imMobileMassGas;
1800 }
1801 if (!this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped].empty()) {
1802 const Scalar mobileMassGas = massGas * std::max(Scalar{0.0}, sg - trappedGasSaturation);
1803 this->fip_[Inplace::Phase::CO2MassInGasPhaseEffectiveUnTrapped][globalDofIdx] = mobileMassGas;
1804 }
1805 }
1806 }
1807
1808 template <typename FluidState>
1809 void updateCO2InWater(const unsigned globalDofIdx,
1810 const double pv,
1811 const FluidState& fs)
1812 {
1813 const auto co2InWater = FluidSystem::phaseIsActive(oilPhaseIdx)
1814 ? this->co2InWaterFromOil(fs, pv)
1815 : this->co2InWaterFromWater(fs, pv);
1816
1817 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1818 if (!this->fip_[Inplace::Phase::CO2Mass].empty()) {
1819 this->fip_[Inplace::Phase::CO2Mass][globalDofIdx] += co2InWater * mM;
1820 }
1821 if (!this->fip_[Inplace::Phase::CO2MassInWaterPhase].empty()) {
1822 this->fip_[Inplace::Phase::CO2MassInWaterPhase][globalDofIdx] = co2InWater * mM;
1823 }
1824 if (!this->fip_[Inplace::Phase::CO2InWaterPhase].empty()) {
1825 this->fip_[Inplace::Phase::CO2InWaterPhase][globalDofIdx] = co2InWater;
1826 }
1827 }
1828
1829 template <typename FluidState>
1830 Scalar co2InWaterFromWater(const FluidState& fs, const double pv) const
1831 {
1832 const double rhow = getValue(fs.density(waterPhaseIdx));
1833 const double sw = getValue(fs.saturation(waterPhaseIdx));
1834 const double xwG = FluidSystem::convertRswToXwG(getValue(fs.Rsw()), fs.pvtRegionIndex());
1835
1836 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1837
1838 return xwG * pv * rhow * sw / mM;
1839 }
1840
1841 template <typename FluidState>
1842 Scalar co2InWaterFromOil(const FluidState& fs, const double pv) const
1843 {
1844 const double rhoo = getValue(fs.density(oilPhaseIdx));
1845 const double so = getValue(fs.saturation(oilPhaseIdx));
1846 const double xoG = FluidSystem::convertRsToXoG(getValue(fs.Rs()), fs.pvtRegionIndex());
1847
1848 const Scalar mM = FluidSystem::molarMass(gasCompIdx, fs.pvtRegionIndex());
1849
1850 return xoG * pv * rhoo * so / mM;
1851 }
1852
1853 const Simulator& simulator_;
1854};
1855
1856} // namespace Opm
1857
1858#endif // OPM_OUTPUT_BLACK_OIL_MODULE_HPP
Helper class for grid instantiation of ECL file-format using problems.
Output module for the results black oil model writing in ECL binary format.
Declares the properties required by the black oil model.
The base class for the element-centered finite-volume discretization scheme.
Definition ecfvdiscretization.hh:147
Definition GenericOutputBlackoilModule.hpp:61
Inter-region flow accumulation maps for all region definition arrays.
Definition InterRegFlows.hpp:179
void compress()
Form CSR adjacency matrix representation of input graph from connections established in previous call...
Definition InterRegFlows.cpp:165
void addConnection(const Cell &source, const Cell &destination, const data::InterRegFlowMap::FlowRates &rates)
Add flow rate connection between regions for all region definitions.
Definition InterRegFlows.cpp:156
void clear()
Clear all internal buffers, but preserve allocated capacity.
Definition InterRegFlows.cpp:172
Output module for the results black oil model writing in ECL binary format.
Definition OutputBlackoilModule.hpp:90
static void registerParameters()
Register all run-time parameters for the Vtk output module.
Definition OutputBlackoilModule.hpp:182
void processElement(const ElementContext &elemCtx)
Modify the internal buffers according to the intensive quanties relevant for an element.
Definition OutputBlackoilModule.hpp:272
void initializeFluxData()
Prepare for capturing connection fluxes, particularly to account for inter-region flows.
Definition OutputBlackoilModule.hpp:1108
void allocBuffers(const unsigned bufferSize, const unsigned reportStepNum, const bool substep, const bool log, const bool isRestart)
Allocate memory for the scalar fields we would like to write to ECL output files.
Definition OutputBlackoilModule.hpp:197
void processFluxes(const ElementContext &elemCtx, ActiveIndex &&activeIndex, CartesianIndex &&cartesianIndex)
Capture connection fluxes, particularly to account for inter-region flows.
Definition OutputBlackoilModule.hpp:1071
const InterRegFlowMap & getInterRegFlows() const
Get read-only access to collection of inter-region flows.
Definition OutputBlackoilModule.hpp:1126
void finalizeFluxData()
Finalize capturing connection fluxes.
Definition OutputBlackoilModule.hpp:1118
Declare the properties used by the infrastructure code of the finite volume discretizations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
std::string moduleVersionName()
Return the version name of the module, for example "2015.10" (for a release branch) or "2016....
Definition moduleVersion.cpp:34
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
This file provides the infrastructure to retrieve run-time parameters.
The Opm property system, traits with inheritance.
Minimal characteristics of a cell from a simulation grid.
Definition InterRegFlows.hpp:50
Definition OutputBlackoilModule.hpp:70
Definition OutputBlackoilModule.hpp:71