From 85513754bce10c4455a64f3b7cff3a8dba78a8c6 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 11 Sep 2024 14:25:34 +0200 Subject: [PATCH] splitting Blackoil related to FlowProblemBlackoil so FlowProblem can be used for compositional or other setting --- CMakeLists_files.cmake | 4 +- flowexperimental/flowexp.hpp | 12 +- flowexperimental/flowexp_blackoil.cpp | 2 +- opm/simulators/flow/BaseAquiferModel.hpp | 7 + opm/simulators/flow/BlackoilModel.hpp | 4 +- ...ties.hpp => FlowBaseProblemProperties.hpp} | 63 +- opm/simulators/flow/FlowProblem.hpp | 1332 +-------------- opm/simulators/flow/FlowProblemBlackoil.hpp | 1465 +++++++++++++++++ .../flow/FlowProblemBlackoilProperties.hpp | 99 ++ opm/simulators/flow/MixingRateControls.cpp | 2 +- opm/simulators/flow/MixingRateControls.hpp | 4 +- opm/simulators/linalg/ISTLSolver.hpp | 2 +- tests/TestTypeTag.hpp | 6 +- tests/test_RestartSerialization.cpp | 2 +- tests/test_equil.cpp | 7 +- tests/test_glift1.cpp | 4 +- tests/test_wellmodel.cpp | 2 +- 17 files changed, 1666 insertions(+), 1351 deletions(-) rename opm/simulators/flow/{FlowProblemProperties.hpp => FlowBaseProblemProperties.hpp} (80%) create mode 100644 opm/simulators/flow/FlowProblemBlackoil.hpp create mode 100644 opm/simulators/flow/FlowProblemBlackoilProperties.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 34d29beb1..052e397b8 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -750,14 +750,16 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/flow/ExtraConvergenceOutputThread.hpp opm/simulators/flow/FemCpGridCompat.hpp opm/simulators/flow/FIBlackoilModel.hpp + opm/simulators/flow/FlowBaseProblemProperties.hpp opm/simulators/flow/FlowBaseVanguard.hpp opm/simulators/flow/FlowGenericProblem.hpp opm/simulators/flow/FlowGenericProblem_impl.hpp opm/simulators/flow/FlowGenericVanguard.hpp opm/simulators/flow/FlowMain.hpp opm/simulators/flow/FlowProblem.hpp + opm/simulators/flow/FlowProblemBlackoil.hpp opm/simulators/flow/FlowProblemParameters.hpp - opm/simulators/flow/FlowProblemProperties.hpp + opm/simulators/flow/FlowProblemBlackoilProperties.hpp opm/simulators/flow/FlowUtils.hpp opm/simulators/flow/FlowsData.hpp opm/simulators/flow/FlowThresholdPressure.hpp diff --git a/flowexperimental/flowexp.hpp b/flowexperimental/flowexp.hpp index 432af650b..c12591b8f 100644 --- a/flowexperimental/flowexp.hpp +++ b/flowexperimental/flowexp.hpp @@ -37,8 +37,8 @@ #include -#include -#include +#include +#include #include @@ -57,7 +57,7 @@ namespace TTag { struct FlowExpTypeTag { - using InheritsFrom = std::tuple; + using InheritsFrom = std::tuple; }; } @@ -126,9 +126,9 @@ struct LinearSolverBackend { namespace Opm { template -class FlowExpProblem : public FlowProblem //, public FvBaseProblem +class FlowExpProblem : public FlowProblemBlackoil //, public FvBaseProblem { - typedef FlowProblem ParentType; + typedef FlowProblemBlackoil ParentType; using BaseType = ParentType; // GetPropType; using Scalar = GetPropType; @@ -203,7 +203,7 @@ public: } // inherit the constructors - using ParentType::FlowProblem; + using ParentType::FlowProblemBlackoil; }; } diff --git a/flowexperimental/flowexp_blackoil.cpp b/flowexperimental/flowexp_blackoil.cpp index 36bbc09b6..2ec5c4cd3 100644 --- a/flowexperimental/flowexp_blackoil.cpp +++ b/flowexperimental/flowexp_blackoil.cpp @@ -23,7 +23,7 @@ #include #include #include -#include +#include #include "BlackOilIntensiveQuantitiesGlobalIndex.hpp" #include "FIBlackOilModelNoCache.hpp" diff --git a/opm/simulators/flow/BaseAquiferModel.hpp b/opm/simulators/flow/BaseAquiferModel.hpp index 0a52d6bd6..ba32fe5d5 100644 --- a/opm/simulators/flow/BaseAquiferModel.hpp +++ b/opm/simulators/flow/BaseAquiferModel.hpp @@ -108,6 +108,13 @@ public: unsigned) const { } + + void addToSource(RateVector&, + unsigned, + unsigned) const + { } + + /*! * \brief This method is called after each Newton-Raphson successful iteration. * diff --git a/opm/simulators/flow/BlackoilModel.hpp b/opm/simulators/flow/BlackoilModel.hpp index 3d602a7fe..4182958b0 100644 --- a/opm/simulators/flow/BlackoilModel.hpp +++ b/opm/simulators/flow/BlackoilModel.hpp @@ -38,7 +38,7 @@ #include #include #include -#include +#include #include #include #include @@ -77,7 +77,7 @@ namespace Opm::Properties { namespace TTag { -struct FlowProblem { using InheritsFrom = std::tuple; }; +struct FlowProblem { using InheritsFrom = std::tuple; }; } diff --git a/opm/simulators/flow/FlowProblemProperties.hpp b/opm/simulators/flow/FlowBaseProblemProperties.hpp similarity index 80% rename from opm/simulators/flow/FlowProblemProperties.hpp rename to opm/simulators/flow/FlowBaseProblemProperties.hpp index 391c787e5..a1982cc07 100644 --- a/opm/simulators/flow/FlowProblemProperties.hpp +++ b/opm/simulators/flow/FlowBaseProblemProperties.hpp @@ -23,12 +23,12 @@ /*! * \file * - * \copydoc Opm::FlowProblem + * \copydoc Opm::FlowBaseProblem */ -#ifndef OPM_FLOW_PROBLEM_PROPERTIES_HPP -#define OPM_FLOW_PROBLEM_PROPERTIES_HPP +#ifndef OPM_FLOW_BASE_PROBLEM_PROPERTIES_HPP +#define OPM_FLOW_BASE_PROBLEM_PROPERTIES_HPP + -#include #include #include @@ -39,10 +39,7 @@ #include #include #include -#include -#include -#include -#include +#include #if HAVE_DAMARIS #include @@ -50,11 +47,6 @@ #include -namespace Opm { -template -class FlowProblem; -} - namespace Opm::Properties { namespace TTag { @@ -89,14 +81,14 @@ struct EnableThermalFluxBoundaries { using type = UndefinedProperty; }; template struct WellModel { using type = UndefinedProperty; }; -// Set the problem property -template -struct Problem -{ using type = FlowProblem; }; +// Tracer might be moved to the blackoil side +// The class that deals with the tracer +template +struct TracerModel { using type = UndefinedProperty; }; -template -struct Model -{ using type = FIBlackOilModel; }; +template +struct TracerModel +{ using type = ::Opm::TracerModel; }; // Select the element centered finite volume method as spatial discretization template @@ -123,25 +115,6 @@ template struct GridView { using type = typename GetPropType::LeafGridView; }; -// Set the material law for fluid fluxes -template -struct MaterialLaw -{ -private: - using Scalar = GetPropType; - using FluidSystem = GetPropType; - - using Traits = ThreePhaseMaterialTraits; - -public: - using EclMaterialLawManager = ::Opm::EclMaterialLawManager; - - using type = typename EclMaterialLawManager::MaterialLaw; -}; - // Set the material law for energy storage in rock template struct SolidEnergyLaw @@ -212,16 +185,6 @@ template struct EnableApiTracking { static constexpr bool value = false; }; -// Use the "velocity module" which uses the Eclipse "NEWTRAN" transmissibilities -template -struct FluxModule -{ using type = NewTranFluxModule; }; - -// Use the dummy gradient calculator in order not to do unnecessary work. -template -struct GradientCalculator -{ using type = DummyGradientCalculator; }; - // store temperature (but do not conserve energy, as long as EnableEnergy is false) template struct EnableTemperature @@ -275,4 +238,4 @@ struct EnableDebuggingChecks } // namespace Opm::Properties -#endif // OPM_FLOW_PROBLEM_PROPERTIES_HPP +#endif // OPM_BASE_FLOW_PROBLEM_PROPERTIES_HPP diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index a7125da59..f49e3be94 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -44,41 +44,23 @@ #include #include #include -#include -#include -#include -#include -#include -#include -#include -#include #include #include #include -#include -#include #include -#include #include -#include -#include #include #include #include #include -#include #include -#include -#include -#include -#include -#include +// TODO: maybe we can name it FlowProblemProperties.hpp +#include #include #include -#include #include #include @@ -88,12 +70,6 @@ #include -#include - -#if HAVE_DAMARIS -#include -#endif - #include #include #include @@ -113,6 +89,7 @@ class FlowProblem : public GetPropType , public FlowGenericProblem, GetPropType> { +protected: using BaseType = FlowGenericProblem, GetPropType>; using ParentType = GetPropType; @@ -134,32 +111,35 @@ class FlowProblem : public GetPropType enum { numEq = getPropValue() }; enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; - enum { enableExperiments = getPropValue() }; - enum { enableSolvent = getPropValue() }; - enum { enablePolymer = getPropValue() }; + + enum { enableConvectiveMixing = getPropValue() }; enum { enableBrine = getPropValue() }; - enum { enableSaltPrecipitation = getPropValue() }; - enum { enablePolymerMolarWeight = getPropValue() }; - enum { enableFoam = getPropValue() }; - enum { enableExtbo = getPropValue() }; - enum { enableTemperature = getPropValue() }; - enum { enableEnergy = getPropValue() }; enum { enableDiffusion = getPropValue() }; enum { enableDispersion = getPropValue() }; - enum { enableConvectiveMixing = getPropValue() }; - enum { enableThermalFluxBoundaries = getPropValue() }; - enum { enableApiTracking = getPropValue() }; + enum { enableEnergy = getPropValue() }; + enum { enableExperiments = getPropValue() }; + enum { enableExtbo = getPropValue() }; + enum { enableFoam = getPropValue() }; enum { enableMICP = getPropValue() }; + enum { enablePolymer = getPropValue() }; + enum { enablePolymerMolarWeight = getPropValue() }; + enum { enableSaltPrecipitation = getPropValue() }; + enum { enableSolvent = getPropValue() }; + enum { enableTemperature = getPropValue() }; + enum { enableThermalFluxBoundaries = getPropValue() }; + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + + // TODO: later, gasCompIdx, oilCompIdx and waterCompIdx should go to the FlowProblemBlackoil in the future + // we do not want them in the compositional setting enum { gasCompIdx = FluidSystem::gasCompIdx }; enum { oilCompIdx = FluidSystem::oilCompIdx }; enum { waterCompIdx = FluidSystem::waterCompIdx }; using PrimaryVariables = GetPropType; using RateVector = GetPropType; - using BoundaryRateVector = GetPropType; using Simulator = GetPropType; using Element = typename GridView::template Codim<0>::Entity; using ElementContext = GetPropType; @@ -176,28 +156,11 @@ class FlowProblem : public GetPropType using WellModel = GetPropType; using AquiferModel = GetPropType; - using SolventModule = BlackOilSolventModule; - using PolymerModule = BlackOilPolymerModule; - using FoamModule = BlackOilFoamModule; - using BrineModule = BlackOilBrineModule; - using ExtboModule = BlackOilExtboModule; - using MICPModule = BlackOilMICPModule; - using DispersionModule = BlackOilDispersionModule; - using DiffusionModule = BlackOilDiffusionModule; - using ConvectiveMixingModule = BlackOilConvectiveMixingModule; - using ModuleParams = typename BlackOilLocalResidualTPFA::ModuleParams; - - using InitialFluidState = typename EquilInitializer::ScalarFluidState; - using Toolbox = MathToolbox; using DimMatrix = Dune::FieldMatrix; - using EclWriterType = EclWriter; -#if HAVE_DAMARIS - using DamarisWriterType = DamarisWriter; -#endif - using TracerModel = ::Opm::TracerModel; + using TracerModel = GetPropType; using DirectionalMobilityPtr = Utility::CopyablePtr>; public: @@ -215,12 +178,7 @@ public: static void registerParameters() { ParentType::registerParameters(); - EclWriterType::registerParameters(); -#if HAVE_DAMARIS - DamarisWriterType::registerParameters(); -#endif - VtkTracerModule::registerParameters(); registerFlowProblemParameters(); } @@ -247,7 +205,7 @@ public: /*! * \copydoc Doxygen::defaultProblemConstructor */ - FlowProblem(Simulator& simulator) + explicit FlowProblem(Simulator& simulator) : ParentType(simulator) , BaseType(simulator.vanguard().eclState(), simulator.vanguard().schedule(), @@ -260,61 +218,15 @@ public: enableEnergy, enableDiffusion, enableDispersion) - , thresholdPressures_(simulator) , wellModel_(simulator) , aquiferModel_(simulator) , pffDofData_(simulator.gridView(), this->elementMapper()) , tracerModel_(simulator) - , mixControls_(simulator.vanguard().schedule()) - , actionHandler_(simulator.vanguard().eclState(), - simulator.vanguard().schedule(), - simulator.vanguard().actionState(), - simulator.vanguard().summaryState(), - wellModel_, - simulator.vanguard().grid().comm()) { - this->model().addOutputModule(new VtkTracerModule(simulator)); - // Tell the black-oil extensions to initialize their internal data structures const auto& vanguard = simulator.vanguard(); - BlackOilBrineParams brineParams; - brineParams.template initFromState(vanguard.eclState()); - BrineModule::setParams(std::move(brineParams)); - - DiffusionModule::initFromState(vanguard.eclState()); - DispersionModule::initFromState(vanguard.eclState()); - - BlackOilExtboParams extboParams; - extboParams.template initFromState(vanguard.eclState()); - ExtboModule::setParams(std::move(extboParams)); - - BlackOilFoamParams foamParams; - foamParams.template initFromState(vanguard.eclState()); - FoamModule::setParams(std::move(foamParams)); - - BlackOilMICPParams micpParams; - micpParams.template initFromState(vanguard.eclState()); - MICPModule::setParams(std::move(micpParams)); - - BlackOilPolymerParams polymerParams; - polymerParams.template initFromState(vanguard.eclState()); - PolymerModule::setParams(std::move(polymerParams)); - - BlackOilSolventParams solventParams; - solventParams.template initFromState(vanguard.eclState(), vanguard.schedule()); - SolventModule::setParams(std::move(solventParams)); - - // create the ECL writer - eclWriter_ = std::make_unique(simulator); -#if HAVE_DAMARIS - // create Damaris writer - damarisWriter_ = std::make_unique(simulator); - enableDamarisOutput_ = Parameters::Get(); -#endif enableDriftCompensation_ = Parameters::Get(); - enableEclOutput_ = Parameters::Get(); enableVtkOutput_ = Parameters::Get(); this->enableTuning_ = Parameters::Get(); @@ -339,149 +251,6 @@ public: relpermDiagnostics.diagnosis(vanguard.eclState(), vanguard.cartesianIndexMapper()); } - /*! - * \copydoc FvBaseProblem::finishInit - */ - void finishInit() - { - ParentType::finishInit(); - auto& simulator = this->simulator(); - - auto finishTransmissibilities = [updated = false, this]() mutable - { - if (updated) { return; } - this->transmissibilities_.finishInit([&vg = this->simulator().vanguard()](const unsigned int it) { - return vg.gridIdxToEquilGridIdx(it); - }); - updated = true; - }; - - // calculating the TRANX, TRANY, TRANZ and NNC for output purpose - // for parallel running, it is based on global trans_ - // for serial running, it is based on the transmissibilities_ - // we try to avoid for the parallel running, has both global trans_ and transmissibilities_ allocated at the same time - if (enableEclOutput_) { - if (simulator.vanguard().grid().comm().size() > 1) { - if (simulator.vanguard().grid().comm().rank() == 0) - eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility()); - } else { - finishTransmissibilities(); - eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities()); - } - - std::function equilGridToGrid = [&simulator](unsigned int i) { - return simulator.vanguard().gridEquilIdxToGridIdx(i); - }; - eclWriter_->extractOutputTransAndNNC(equilGridToGrid); - } - simulator.vanguard().releaseGlobalTransmissibilities(); - - const auto& eclState = simulator.vanguard().eclState(); - const auto& schedule = simulator.vanguard().schedule(); - - // Set the start time of the simulation - simulator.setStartTime(schedule.getStartTime()); - simulator.setEndTime(schedule.simTime(schedule.size() - 1)); - - // We want the episode index to be the same as the report step index to make - // things simpler, so we have to set the episode index to -1 because it is - // incremented by endEpisode(). The size of the initial time step and - // length of the initial episode is set to zero for the same reason. - simulator.setEpisodeIndex(-1); - simulator.setEpisodeLength(0.0); - - // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false - // disables gravity, else the standard value of the gravity constant at sea level - // on earth is used - this->gravity_ = 0.0; - if (Parameters::Get()) - this->gravity_[dim - 1] = 9.80665; - if (!eclState.getInitConfig().hasGravity()) - this->gravity_[dim - 1] = 0.0; - - if (this->enableTuning_) { - // if support for the TUNING keyword is enabled, we get the initial time - // steping parameters from it instead of from command line parameters - const auto& tuning = schedule[0].tuning(); - this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0; - this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC; - } - - this->initFluidSystem_(); - - // deal with DRSDT - this->mixControls_.init(this->model().numGridDof(), - this->episodeIndex(), - eclState.runspec().tabdims().getNumPVTTables()); - - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && - FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { - this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0); - } - - this->readRockParameters_(simulator.vanguard().cellCenterDepths(), - [&simulator](const unsigned idx) - { - std::array coords; - simulator.vanguard().cartesianCoordinate(idx, coords); - for (auto& c : coords) { - ++c; - } - return coords; - }); - readMaterialParameters_(); - readThermalParameters_(); - - // write the static output files (EGRID, INIT) - if (enableEclOutput_) { - eclWriter_->writeInit(); - } - - finishTransmissibilities(); - - const auto& initconfig = eclState.getInitConfig(); - tracerModel_.init(initconfig.restartRequested()); - if (initconfig.restartRequested()) - readEclRestartSolution_(); - else - readInitialCondition_(); - - tracerModel_.prepareTracerBatches(); - - updatePffDofData_(); - - if constexpr (getPropValue()) { - const auto& vanguard = this->simulator().vanguard(); - const auto& gridView = vanguard.gridView(); - int numElements = gridView.size(/*codim=*/0); - this->polymer_.maxAdsorption.resize(numElements, 0.0); - } - - readBoundaryConditions_(); - - // compute and set eq weights based on initial b values - computeAndSetEqWeights_(); - - if (enableDriftCompensation_) { - drift_.resize(this->model().numGridDof()); - drift_ = 0.0; - } - - if (enableVtkOutput_ && eclState.getIOConfig().initOnly()) { - simulator.setTimeStepSize(0.0); - ParentType::writeOutput(true); - } - - // after finishing the initialization and writing the initial solution, we move - // to the first "real" episode/report step - // for restart the episode index and start is already set - if (!initconfig.restartRequested()) { - simulator.startNextEpisode(schedule.seconds(1)); - simulator.setEpisodeIndex(0); - simulator.setTimeStepIndex(0); - } - } - void prefetch(const Element& elem) const { pffDofData_.prefetch(elem); } @@ -500,7 +269,7 @@ public: void deserialize(Restarter& res) { // reload the current episode/report step from the deck - beginEpisode(); + this->beginEpisode(); // deserialize the wells wellModel_.deserialize(res); @@ -531,7 +300,7 @@ public: /*! * \brief Called by the simulator before an episode begins. */ - void beginEpisode() + virtual void beginEpisode() { OPM_TIMEBLOCK(beginEpisode); // Proceed to the next report step @@ -568,10 +337,6 @@ public: } bool tuningEvent = this->beginEpisode_(enableExperiments, this->episodeIndex()); - this->mixControls_.init(this->model().numGridDof(), - episodeIdx, - eclState.runspec().tabdims().getNumPVTTables()); - // set up the wells for the next episode. wellModel_.beginEpisode(); @@ -587,21 +352,6 @@ public: // if TUNING is enabled, also limit the time step size after a tuning event to TSINIT dt = std::min(dt, this->initialTimeStepSize_); simulator.setTimeStepSize(dt); - - // Evaluate UDQ assign statements to make sure the settings are - // available as UDA controls for the current report step. - actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState()); - - if (episodeIdx >= 0) { - const auto& oilVap = schedule[episodeIdx].oilvap(); - if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) { - FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2()); - } else { - FluidSystem::setVapPars(0.0, 0.0); - } - } - - ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), simulator.vanguard().schedule(), episodeIdx, moduleParams_.convectiveMixingModuleParam); } /*! @@ -624,7 +374,7 @@ public: // update maximum water saturation and minimum pressure // used when ROCKCOMP is activated // Do not update max RS first step after a restart - asImp_().updateExplicitQuantities_(episodeIdx, timeStepSize, first_step_ && (episodeIdx > 0)); + this->updateExplicitQuantities_(episodeIdx, timeStepSize, first_step_ && (episodeIdx > 0)); first_step_ = false; if (nonTrivialBoundaryConditions()) { @@ -660,7 +410,7 @@ public: /*! * \brief Called by the simulator after each time integration. */ - void endTimeStep() + virtual void endTimeStep() { OPM_TIMEBLOCK(endTimeStep); @@ -706,72 +456,15 @@ public: } } } - - const bool isSubStep = !this->simulator().episodeWillBeOver(); - - // For CpGrid with LGRs, ecl/vtk output is not supported yet. - const auto& grid = this->simulator().vanguard().gridView().grid(); - - using GridType = std::remove_cv_t>; - constexpr bool isCpGrid = std::is_same_v; - if (!isCpGrid || (grid.maxLevel() == 0)) { - this->eclWriter_->evalSummaryState(isSubStep); - } - - { - OPM_TIMEBLOCK(applyActions); - - const int episodeIdx = this->episodeIndex(); - - // Re-ordering in case of Alugrid - this->actionHandler_ - .applyActions(episodeIdx, simulator.time() + simulator.timeStepSize(), - [this](const bool global) - { - using TransUpdateQuantities = typename Vanguard::TransmissibilityType::TransUpdateQuantities; - this->transmissibilities_ - .update(global, TransUpdateQuantities::All, [&vg = this->simulator().vanguard()] - (const unsigned int i) - { - return vg.gridIdxToEquilGridIdx(i); - }); - }); - } - - // Deal with "clogging" for the MICP model - if constexpr (enableMICP) { - auto& model = this->model(); - const auto& residual = model.linearizer().residual(); - - for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); ++globalDofIdx) { - auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx]; - MICPModule::checkCloggingMICP(model, phi, globalDofIdx); - } - } } /*! * \brief Called by the simulator after the end of an episode. */ - void endEpisode() + virtual void endEpisode() { - OPM_TIMEBLOCK(endEpisode); - const int episodeIdx = this->episodeIndex(); - // Rerun UDQ assignents following action processing on the final - // time step of this episode to make sure that any UDQ ASSIGN - // operations triggered in action blocks take effect. This is - // mainly to work around a shortcoming of the ScheduleState copy - // constructor which clears pending UDQ assignments under the - // assumption that all such assignments have been processed. If an - // action block happens to trigger on the final time step of an - // episode and that action block runs a UDQ assignment, then that - // assignment would be dropped and the rest of the simulator will - // never see its effect without this hack. - this->actionHandler_ - .evalUDQAssignments(episodeIdx, this->simulator().vanguard().udqState()); - this->wellModel_.endEpisode(); this->aquiferModel_.endEpisode(); @@ -791,7 +484,7 @@ public: * \brief Write the requested quantities of the current solution into the output * files. */ - void writeOutput(const SimulatorTimer& timer, bool verbose = true) + void writeOutput(bool verbose = true) { OPM_TIMEBLOCK(problemWriteOutput); // use the generic code to prepare the output fields and to @@ -800,30 +493,8 @@ public: this->simulator().episodeWillBeOver()) { ParentType::writeOutput(verbose); } - - bool isSubStep = !this->simulator().episodeWillBeOver(); - - data::Solution localCellData = {}; -#if HAVE_DAMARIS - // N.B. the Damaris output has to be done before the ECL output as the ECL one - // does all kinds of std::move() relocation of data - if (enableDamarisOutput_) { - damarisWriter_->writeOutput(localCellData, isSubStep) ; - } -#endif - if (enableEclOutput_){ - eclWriter_->writeOutput(std::move(localCellData), timer, isSubStep); - } } - void finalizeOutput() { - OPM_TIMEBLOCK(finalizeOutput); - // this will write all pending output to disk - // to avoid corruption of output files - eclWriter_.reset(); - } - - /*! * \copydoc FvBaseMultiPhaseProblem::intrinsicPermeability */ @@ -976,17 +647,6 @@ public: const typename Vanguard::TransmissibilityType& eclTransmissibilities() const { return transmissibilities_; } - /*! - * \copydoc BlackOilBaseProblem::thresholdPressure - */ - Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const - { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); } - - const FlowThresholdPressure& thresholdPressure() const - { return thresholdPressures_; } - - FlowThresholdPressure& thresholdPressure() - { return thresholdPressures_; } const TracerModel& tracerModel() const { return tracerModel_; } @@ -1175,6 +835,7 @@ public: unsigned plmixnumRegionIndex(const Context& context, unsigned spaceIdx, unsigned timeIdx) const { return this->plmixnumRegionIndex(context.globalSpaceIndex(spaceIdx, timeIdx)); } + // TODO: polymer related might need to go to the blackoil side using BaseType::maxPolymerAdsorption; /*! * \brief Returns the max polymer adsorption value @@ -1198,7 +859,7 @@ public: // use the initial temperature of the DOF if temperature is not a primary // variable unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - return initialFluidStates_[globalDofIdx].temperature(/*phaseIdx=*/0); + return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0); } @@ -1206,7 +867,7 @@ public: { // use the initial temperature of the DOF if temperature is not a primary // variable - return initialFluidStates_[globalDofIdx].temperature(/*phaseIdx=*/0); + return asImp_().initialFluidState(globalDofIdx).temperature(/*phaseIdx=*/0); } const SolidEnergyLawParams& @@ -1222,49 +883,6 @@ public: return this->thermalLawManager_->thermalConductionLawParams(globalSpaceIdx); } - /*! - * \copydoc FvBaseProblem::boundary - * - * Reservoir simulation uses no-flow conditions as default for all boundaries. - */ - template - void boundary(BoundaryRateVector& values, - const Context& context, - unsigned spaceIdx, - unsigned timeIdx) const - { - OPM_TIMEBLOCK_LOCAL(eclProblemBoundary); - if (!context.intersection(spaceIdx).boundary()) - return; - - if constexpr (!enableEnergy || !enableThermalFluxBoundaries) - values.setNoFlow(); - else { - // in the energy case we need to specify a non-trivial boundary condition - // because the geothermal gradient needs to be maintained. for this, we - // simply assume the initial temperature at the boundary and specify the - // thermal flow accordingly. in this context, "thermal flow" means energy - // flow due to a temerature gradient while assuming no-flow for mass - unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx); - unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx); - values.setThermalFlow(context, spaceIdx, timeIdx, initialFluidStates_[globalDofIdx]); - } - - if (nonTrivialBoundaryConditions()) { - unsigned indexInInside = context.intersection(spaceIdx).indexInInside(); - unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx); - unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx); - unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx); - const auto [type, massrate] = boundaryCondition(globalDofIdx, indexInInside); - if (type == BCType::THERMAL) - values.setThermalFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside)); - else if (type == BCType::FREE || type == BCType::DIRICHLET) - values.setFreeFlow(context, spaceIdx, timeIdx, boundaryFluidState(globalDofIdx, indexInInside)); - else if (type == BCType::RATE) - values.setMassRate(massrate, pvtRegionIdx); - } - } - /*! * \brief Returns an element's historic maximum oil phase saturation that was * observed during the simulation. @@ -1299,93 +917,10 @@ public: this->maxOilSaturation_[globalDofIdx] = value; } - /*! - * \brief Returns the maximum value of the gas dissolution factor at the current time - * for a given degree of freedom. - */ - Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const - { - return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx, - this->episodeIndex(), - this->pvtRegionIndex(globalDofIdx)); - } - - /*! - * \brief Returns the maximum value of the oil vaporization factor at the current - * time for a given degree of freedom. - */ - Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const - { - return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx, - this->episodeIndex(), - this->pvtRegionIndex(globalDofIdx)); - } - - /*! - * \brief Return if the storage term of the first iteration is identical to the storage - * term for the solution of the previous time step. - * - * For quite technical reasons, the storage term cannot be recycled if either DRSDT - * or DRVDT are active. Nor if the porosity is changes between timesteps - * using a pore volume multiplier (i.e., poreVolumeMultiplier() != 1.0) - */ - bool recycleFirstIterationStorage() const - { - int episodeIdx = this->episodeIndex(); - return !this->mixControls_.drsdtActive(episodeIdx) && - !this->mixControls_.drvdtActive(episodeIdx) && - this->rockCompPoroMultWc_.empty() && - this->rockCompPoroMult_.empty(); - } - - /*! - * \copydoc FvBaseProblem::initial - * - * The reservoir problem uses a constant boundary condition for - * the whole domain. - */ - template - void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const - { - unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx); - - values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx)); - values.assignNaive(initialFluidStates_[globalDofIdx]); - - SolventModule::assignPrimaryVars(values, - enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0, - enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0); - - if constexpr (enablePolymer) - values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx]; - - if constexpr (enablePolymerMolarWeight) - values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx]; - - if constexpr (enableBrine) { - if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) { - values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation(); - } - else { - values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration(); - } - } - - if constexpr (enableMICP){ - values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx]; - values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx]; - values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx]; - values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx]; - values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx]; - } - - values.checkDefined(); - } - /*! * \copydoc FvBaseProblem::initialSolutionApplied() */ - void initialSolutionApplied() + virtual void initialSolutionApplied() { // Calculate all intensive quantities. this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx*/0); @@ -1403,17 +938,8 @@ public: // the well model uses these... wellModel_.init(); - // let the object for threshold pressures initialize itself. this is done only at - // this point, because determining the threshold pressures may require to access - // the initial solution. - thresholdPressures_.finishInit(); - aquiferModel_.initialSolutionApplied(); - if (this->simulator().episodeIndex() == 0) { - eclWriter_->writeInitialFIPReport(); - } - const bool invalidateFromHyst = updateHysteresis_(); if (invalidateFromHyst) { OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities); @@ -1459,101 +985,9 @@ public: addToSourceDense(rate, globalDofIdx, timeIdx); } - void addToSourceDense(RateVector& rate, - unsigned globalDofIdx, - unsigned timeIdx) const - { - aquiferModel_.addToSource(rate, globalDofIdx, timeIdx); - - // Add source term from deck - const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source(); - std::array ijk; - this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk); - - if (source.hasSource(ijk)) { - const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx); - static std::array sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS}; - static std::array phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}; - static std::array cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx}; - - for (unsigned i = 0; i < phidx_map.size(); ++i) { - const auto phaseIdx = phidx_map[i]; - const auto sourceComp = sc_map[i]; - const auto compIdx = cidx_map[i]; - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx); - if constexpr (getPropValue()) { - mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx); - } - rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate; - } - - if constexpr (enableSolvent) { - Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx); - if constexpr (getPropValue()) { - const auto& solventPvt = SolventModule::solventPvt(); - mass_rate /= solventPvt.referenceDensity(pvtRegionIdx); - } - rate[Indices::contiSolventEqIdx] += mass_rate; - } - if constexpr (enablePolymer) { - rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx); - } - if constexpr (enableEnergy) { - for (unsigned i = 0; i < phidx_map.size(); ++i) { - const auto phaseIdx = phidx_map[i]; - if (!FluidSystem::phaseIsActive(phaseIdx)) { - continue; - } - const auto sourceComp = sc_map[i]; - if (source.hasHrate({ijk, sourceComp})) { - rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx); - } else { - const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, /*timeIdx*/ 0); - auto fs = intQuants.fluidState(); - // if temperature is not set, use cell temperature as default - if (source.hasTemperature({ijk, sourceComp})) { - Scalar temperature = source.temperature({ijk, sourceComp}); - fs.setTemperature(temperature); - } - const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx); - Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx); - Scalar energy_rate = getValue(h)*mass_rate; - rate[Indices::contiEnergyEqIdx] += energy_rate; - } - } - } - } - - // if requested, compensate systematic mass loss for cells which were "well - // behaved" in the last time step - if (enableDriftCompensation_) { - const auto& simulator = this->simulator(); - const auto& model = this->model(); - - // we use a lower tolerance for the compensation too - // assure the added drift from the last step does not - // cause convergence issues on the current step - Scalar maxCompensation = model.newtonMethod().tolerance()/10; - Scalar poro = this->porosity(globalDofIdx, timeIdx); - Scalar dt = simulator.timeStepSize(); - EqVector dofDriftRate = drift_[globalDofIdx]; - dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx); - - // restrict drift compensation to the CNV tolerance - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) { - Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro; - if (cnv > maxCompensation) { - dofDriftRate[eqIdx] *= maxCompensation/cnv; - } - } - - for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) - rate[eqIdx] -= dofDriftRate[eqIdx]; - } - } + virtual void addToSourceDense(RateVector& rate, + unsigned globalDofIdx, + unsigned timeIdx) const = 0; /*! * \brief Returns a reference to the ECL well manager used by the problem. @@ -1572,127 +1006,9 @@ public: AquiferModel& mutableAquiferModel() { return aquiferModel_; } - // temporary solution to facilitate output of initial state from flow - const InitialFluidState& initialFluidState(unsigned globalDofIdx) const - { return initialFluidStates_[globalDofIdx]; } - - const EclipseIO& eclIO() const - { return eclWriter_->eclIO(); } - - void setSubStepReport(const SimulatorReportSingle& report) - { return eclWriter_->setSubStepReport(report); } - - void setSimulationReport(const SimulatorReport& report) - { return eclWriter_->setSimulationReport(report); } - bool nonTrivialBoundaryConditions() const { return nonTrivialBoundaryConditions_; } - const InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const - { - OPM_TIMEBLOCK_LOCAL(boundaryFluidState); - const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop; - if (bcprop.size() > 0) { - FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId); - - // index == 0: no boundary conditions for this - // global cell and direction - if (bcindex_(dir)[globalDofIdx] == 0) - return initialFluidStates_[globalDofIdx]; - - const auto& bc = bcprop[bcindex_(dir)[globalDofIdx]]; - if (bc.bctype == BCType::DIRICHLET ) - { - InitialFluidState fluidState; - const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx); - fluidState.setPvtRegionIndex(pvtRegionIdx); - - switch (bc.component) { - case BCComponent::OIL: - if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) - throw std::logic_error("oil is not active and you're trying to add oil BC"); - - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0); - break; - case BCComponent::GAS: - if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) - throw std::logic_error("gas is not active and you're trying to add gas BC"); - - fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0); - break; - case BCComponent::WATER: - if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) - throw std::logic_error("water is not active and you're trying to add water BC"); - - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0); - break; - case BCComponent::SOLVENT: - case BCComponent::POLYMER: - case BCComponent::NONE: - throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC"); - break; - } - fluidState.setTotalSaturation(1.0); - double pressure = initialFluidStates_[globalDofIdx].pressure(refPressurePhaseIdx_()); - const auto pressure_input = bc.pressure; - if (pressure_input) { - pressure = *pressure_input; - } - - std::array pc = {0}; - const auto& matParams = materialLawParams(globalDofIdx); - MaterialLaw::capillaryPressures(pc, matParams, fluidState); - Valgrind::CheckDefined(pressure); - Valgrind::CheckDefined(pc); - for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) { - const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx); - - fluidState.setPc(phaseIdx, pc[phaseIdx]); - if (Indices::oilEnabled) - fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx])); - else if (Indices::gasEnabled) - fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx])); - else if (Indices::waterEnabled) - //single (water) phase - fluidState.setPressure(phaseIdx, pressure); - } - - double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature - const auto temperature_input = bc.temperature; - if(temperature_input) - temperature = *temperature_input; - fluidState.setTemperature(temperature); - - if (FluidSystem::enableDissolvedGas()) { - fluidState.setRs(0.0); - fluidState.setRv(0.0); - } - if (FluidSystem::enableDissolvedGasInWater()) { - fluidState.setRsw(0.0); - } - if (FluidSystem::enableVaporizedWater()) - fluidState.setRvw(0.0); - - for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) { - const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx); - - const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx); - fluidState.setInvB(phaseIdx, b); - - const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx); - fluidState.setDensity(phaseIdx, rho); - if (enableEnergy) { - const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx); - fluidState.setEnthalpy(phaseIdx, h); - } - } - fluidState.checkDefined(); - return fluidState; - } - } - return initialFluidStates_[globalDofIdx]; - } - /*! * \brief Propose the size of the next time step to the simulator. * @@ -1755,7 +1071,7 @@ public: // water compaction assert(!this->rockCompPoroMultWc_.empty()); LhsEval SwMax = max(decay(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]); - LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx); + LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx); return this->rockCompPoroMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true); } @@ -1773,25 +1089,6 @@ public: : this->simulator().problem().getRockCompTransMultVal(elementIdx); } - /*! - * \brief Calculate the transmissibility multiplier due to porosity reduction. - * - * TODO: The API of this is a bit ad-hoc, it would be better to use context objects. - */ - template - LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants) const - { - OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier); - if (!enableSaltPrecipitation) - return 1.0; - - const auto& fs = intQuants.fluidState(); - unsigned tableIdx = fs.pvtRegionIndex(); - LhsEval porosityFactor = decay(1. - fs.saltSaturation()); - porosityFactor = min(porosityFactor, 1.0); - const auto& permfactTable = BrineModule::permfactTable(tableIdx); - return permfactTable.eval(porosityFactor, /*extrapolation=*/true); - } /*! * \brief Return the well transmissibility multiplier due to rock changues. @@ -1840,16 +1137,10 @@ public: rate[Indices::canonicalToActiveComponentIndex(waterCompIdx)] = bc.rate; break; case BCComponent::SOLVENT: - if constexpr (!enableSolvent) - throw std::logic_error("solvent is disabled and you're trying to add solvent to BC"); - - rate[Indices::solventSaturationIdx] = bc.rate; + this->handleSolventBC(bc, rate); break; case BCComponent::POLYMER: - if constexpr (!enablePolymer) - throw std::logic_error("polymer is disabled and you're trying to add polymer to BC"); - - rate[Indices::polymerConcentrationIdx] = bc.rate; + this->handlePolymerBC(bc, rate); break; case BCComponent::NONE: throw std::logic_error("you need to specify the component when RATE type is set in BC"); @@ -1859,19 +1150,6 @@ public: return {bc.bctype, rate}; } - const std::unique_ptr& eclWriter() const - { - return eclWriter_; - } - - void setConvData(const std::vector>& data) - { - eclWriter_->mutableOutputModule().setCnvData(data); - } - - const ModuleParams& moduleParams() const { - return moduleParams_; - } template void serializeOp(Serializer& serializer) @@ -1881,49 +1159,17 @@ public: serializer(wellModel_); serializer(aquiferModel_); serializer(tracerModel_); - serializer(mixControls_); serializer(*materialLawManager_); - serializer(*eclWriter_); } - Scalar drsdtcon(unsigned elemIdx, int episodeIdx) const - { - return this->mixControls_.drsdtcon(elemIdx, episodeIdx, - this->pvtRegionIndex(elemIdx)); - } private: Implementation& asImp_() { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + protected: - void updateExplicitQuantities_(int episodeIdx, int timeStepSize, const bool first_step_after_restart = false) - { - OPM_TIMEBLOCK(updateExplicitQuantities); - const bool invalidateFromMaxWaterSat = updateMaxWaterSaturation_(); - const bool invalidateFromMinPressure = updateMinPressure_(); - - // update hysteresis and max oil saturation used in vappars - const bool invalidateFromHyst = updateHysteresis_(); - const bool invalidateFromMaxOilSat = updateMaxOilSaturation_(); - - - // deal with DRSDT and DRVDT - const bool invalidateDRDT = !first_step_after_restart && this->asImp_().updateCompositionChangeLimits_(); - - // the derivatives may have change - bool invalidateIntensiveQuantities - = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat || invalidateDRDT; - if (invalidateIntensiveQuantities) { - OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities); - this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); - } - - if constexpr (getPropValue()) - updateMaxPolymerAdsorption_(); - - updateRockCompTransMultVal_(); - mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize); - } - template void updateProperty_(const std::string& failureMsg, UpdateFunc func) @@ -1944,41 +1190,6 @@ protected: OPM_END_PARALLEL_TRY_CATCH(failureMsg, vanguard.grid().comm()); } - // update the parameters needed for DRSDT and DRVDT - bool updateCompositionChangeLimits_() - { - OPM_TIMEBLOCK(updateCompositionChangeLimits); - // update the "last Rs" values for all elements, including the ones in the ghost - // and overlap regions - int episodeIdx = this->episodeIndex(); - std::array active{this->mixControls_.drsdtConvective(episodeIdx), - this->mixControls_.drsdtActive(episodeIdx), - this->mixControls_.drvdtActive(episodeIdx)}; - if (!active[0] && !active[1] && !active[2]) { - return false; - } - - this->updateProperty_("FlowProblem::updateCompositionChangeLimits_()) failed:", - [this,episodeIdx,active](unsigned compressedDofIdx, - const IntensiveQuantities& iq) - { - const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx); - const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0; - const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx); - this->mixControls_.update(compressedDofIdx, - iq, - episodeIdx, - this->gravity_[dim - 1], - perm[dim - 1][dim - 1], - distZ, - pvtRegionIdx); - } - ); - - return true; - } - - bool updateMaxOilSaturation_() { OPM_TIMEBLOCK(updateMaxOilSaturation); @@ -2173,8 +1384,9 @@ protected: } } - void readInitialCondition_() + virtual void readInitialCondition_() { + // TODO: whether we should move this to FlowProblemBlackoil const auto& simulator = this->simulator(); const auto& vanguard = simulator.vanguard(); const auto& eclState = vanguard.eclState(); @@ -2184,17 +1396,10 @@ protected: else readExplicitInitialCondition_(); - if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP) - this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(), - enableSolvent, - enablePolymer, - enablePolymerMolarWeight, - enableMICP); - //initialize min/max values std::size_t numElems = this->model().numGridDof(); for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { - const auto& fs = initialFluidStates_[elemIdx]; + const auto& fs = asImp_().initialFluidStates()[elemIdx]; if (!this->maxWaterSaturation_.empty()) this->maxWaterSaturation_[elemIdx] = std::max(this->maxWaterSaturation_[elemIdx], fs.saturation(waterPhaseIdx)); if (!this->maxOilSaturation_.empty()) @@ -2202,373 +1407,10 @@ protected: if (!this->minRefPressure_.empty()) this->minRefPressure_[elemIdx] = std::min(this->minRefPressure_[elemIdx], fs.pressure(refPressurePhaseIdx_())); } - - } - void readEquilInitialCondition_() - { - const auto& simulator = this->simulator(); - - // initial condition corresponds to hydrostatic conditions. - EquilInitializer equilInitializer(simulator, *materialLawManager_); - - std::size_t numElems = this->model().numGridDof(); - initialFluidStates_.resize(numElems); - for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { - auto& elemFluidState = initialFluidStates_[elemIdx]; - elemFluidState.assign(equilInitializer.initialFluidState(elemIdx)); - } - } - - void readEclRestartSolution_() - { - // Throw an exception if the grid has LGRs. Refined grid are not supported for restart. - if(this->simulator().vanguard().grid().maxLevel() > 0) { - throw std::invalid_argument("Refined grids are not yet supported for restart "); - } - - // Set the start time of the simulation - auto& simulator = this->simulator(); - const auto& schedule = simulator.vanguard().schedule(); - const auto& eclState = simulator.vanguard().eclState(); - const auto& initconfig = eclState.getInitConfig(); - const int restart_step = initconfig.getRestartStep(); - { - simulator.setTime(schedule.seconds(restart_step)); - - simulator.startNextEpisode(simulator.startTime() + simulator.time(), - schedule.stepLength(restart_step)); - simulator.setEpisodeIndex(restart_step); - } - eclWriter_->beginRestart(); - - Scalar dt = std::min(eclWriter_->restartTimeStepSize(), simulator.episodeLength()); - simulator.setTimeStepSize(dt); - - std::size_t numElems = this->model().numGridDof(); - initialFluidStates_.resize(numElems); - if constexpr (enableSolvent) { - this->solventSaturation_.resize(numElems, 0.0); - this->solventRsw_.resize(numElems, 0.0); - } - - if constexpr (enablePolymer) - this->polymer_.concentration.resize(numElems, 0.0); - - if constexpr (enablePolymerMolarWeight) { - const std::string msg {"Support of the RESTART for polymer molecular weight " - "is not implemented yet. The polymer weight value will be " - "zero when RESTART begins"}; - OpmLog::warning("NO_POLYMW_RESTART", msg); - this->polymer_.moleWeight.resize(numElems, 0.0); - } - - if constexpr (enableMICP) { - this->micp_.resize(numElems); - } - - // Initialize mixing controls before trying to set any lastRx valuesx - this->mixControls_.init(numElems, restart_step, eclState.runspec().tabdims().getNumPVTTables()); - - for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { - auto& elemFluidState = initialFluidStates_[elemIdx]; - elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx)); - eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx); - eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx); - - // Note: Function processRestartSaturations_() mutates the - // 'ssol' argument--the value from the restart file--if solvent - // is enabled. Then, store the updated solvent saturation into - // 'solventSaturation_'. Otherwise, just pass a dummy value to - // the function and discard the unchanged result. Do not index - // into 'solventSaturation_' unless solvent is enabled. - { - auto ssol = enableSolvent - ? eclWriter_->outputModule().getSolventSaturation(elemIdx) - : Scalar(0); - - processRestartSaturations_(elemFluidState, ssol); - - if constexpr (enableSolvent) { - this->solventSaturation_[elemIdx] = ssol; - this->solventRsw_[elemIdx] = eclWriter_->outputModule().getSolventRsw(elemIdx); - } - } - - // For CO2STORE and H2STORE we need to set the initial temperature for isothermal simulations - bool isThermal = eclState.getSimulationConfig().isThermal(); - bool needTemperature = (eclState.runspec().co2Storage() || eclState.runspec().h2Storage()); - if (!isThermal && needTemperature) { - const auto& fp = simulator.vanguard().eclState().fieldProps(); - elemFluidState.setTemperature(fp.get_double("TEMPI")[elemIdx]); - } - - this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv()); - - if constexpr (enablePolymer) - this->polymer_.concentration[elemIdx] = eclWriter_->outputModule().getPolymerConcentration(elemIdx); - if constexpr (enableMICP){ - this->micp_.microbialConcentration[elemIdx] = eclWriter_->outputModule().getMicrobialConcentration(elemIdx); - this->micp_.oxygenConcentration[elemIdx] = eclWriter_->outputModule().getOxygenConcentration(elemIdx); - this->micp_.ureaConcentration[elemIdx] = eclWriter_->outputModule().getUreaConcentration(elemIdx); - this->micp_.biofilmConcentration[elemIdx] = eclWriter_->outputModule().getBiofilmConcentration(elemIdx); - this->micp_.calciteConcentration[elemIdx] = eclWriter_->outputModule().getCalciteConcentration(elemIdx); - } - // if we need to restart for polymer molecular weight simulation, we need to add related here - } - - const int episodeIdx = this->episodeIndex(); - this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize()); - - // assign the restart solution to the current solution. note that we still need - // to compute real initial solution after this because the initial fluid states - // need to be correct for stuff like boundary conditions. - auto& sol = this->model().solution(/*timeIdx=*/0); - const auto& gridView = this->gridView(); - ElementContext elemCtx(simulator); - for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { - elemCtx.updatePrimaryStencil(elem); - int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0); - } - - // make sure that the ghost and overlap entities exhibit the correct - // solution. alternatively, this could be done in the loop above by also - // considering non-interior elements. Since the initial() method might not work - // 100% correctly for such elements, let's play safe and explicitly synchronize - // using message passing. - this->model().syncOverlap(); - - eclWriter_->endRestart(); - } - - void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation) - { - // each phase needs to be above certain value to be claimed to be existing - // this is used to recover some RESTART running with the defaulted single-precision format - const Scalar smallSaturationTolerance = 1.e-6; - Scalar sumSaturation = 0.0; - for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (FluidSystem::phaseIsActive(phaseIdx)) { - if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance) - elemFluidState.setSaturation(phaseIdx, 0.0); - - sumSaturation += elemFluidState.saturation(phaseIdx); - } - - } - if constexpr (enableSolvent) { - if (solventSaturation < smallSaturationTolerance) - solventSaturation = 0.0; - - sumSaturation += solventSaturation; - } - - assert(sumSaturation > 0.0); - - for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (FluidSystem::phaseIsActive(phaseIdx)) { - const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation; - elemFluidState.setSaturation(phaseIdx, saturation); - } - } - if constexpr (enableSolvent) { - solventSaturation = solventSaturation / sumSaturation; - } - } - - void readExplicitInitialCondition_() - { - const auto& simulator = this->simulator(); - const auto& vanguard = simulator.vanguard(); - const auto& eclState = vanguard.eclState(); - const auto& fp = eclState.fieldProps(); - bool has_swat = fp.has_double("SWAT"); - bool has_sgas = fp.has_double("SGAS"); - bool has_rs = fp.has_double("RS"); - bool has_rv = fp.has_double("RV"); - bool has_rvw = fp.has_double("RVW"); - bool has_pressure = fp.has_double("PRESSURE"); - bool has_salt = fp.has_double("SALT"); - bool has_saltp = fp.has_double("SALTP"); - - // make sure all required quantities are enables - if (Indices::numPhases > 1) { - if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat) - throw std::runtime_error("The ECL input file requires the presence of the SWAT keyword if " - "the water phase is active"); - if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx)) - throw std::runtime_error("The ECL input file requires the presence of the SGAS keyword if " - "the gas phase is active"); - } - if (!has_pressure) - throw std::runtime_error("The ECL input file requires the presence of the PRESSURE " - "keyword if the model is initialized explicitly"); - if (FluidSystem::enableDissolvedGas() && !has_rs) - throw std::runtime_error("The ECL input file requires the RS keyword to be present if" - " dissolved gas is enabled"); - if (FluidSystem::enableVaporizedOil() && !has_rv) - throw std::runtime_error("The ECL input file requires the RV keyword to be present if" - " vaporized oil is enabled"); - if (FluidSystem::enableVaporizedWater() && !has_rvw) - throw std::runtime_error("The ECL input file requires the RVW keyword to be present if" - " vaporized water is enabled"); - if (enableBrine && !has_salt) - throw std::runtime_error("The ECL input file requires the SALT keyword to be present if" - " brine is enabled and the model is initialized explicitly"); - if (enableSaltPrecipitation && !has_saltp) - throw std::runtime_error("The ECL input file requires the SALTP keyword to be present if" - " salt precipitation is enabled and the model is initialized explicitly"); - - std::size_t numDof = this->model().numGridDof(); - - initialFluidStates_.resize(numDof); - - std::vector waterSaturationData; - std::vector gasSaturationData; - std::vector pressureData; - std::vector rsData; - std::vector rvData; - std::vector rvwData; - std::vector tempiData; - std::vector saltData; - std::vector saltpData; - - if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1) - waterSaturationData = fp.get_double("SWAT"); - else - waterSaturationData.resize(numDof); - - if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx)) - gasSaturationData = fp.get_double("SGAS"); - else - gasSaturationData.resize(numDof); - - pressureData = fp.get_double("PRESSURE"); - if (FluidSystem::enableDissolvedGas()) - rsData = fp.get_double("RS"); - - if (FluidSystem::enableVaporizedOil()) - rvData = fp.get_double("RV"); - - if (FluidSystem::enableVaporizedWater()) - rvwData = fp.get_double("RVW"); - - // initial reservoir temperature - tempiData = fp.get_double("TEMPI"); - - // initial salt concentration data - if constexpr (enableBrine) - saltData = fp.get_double("SALT"); - - // initial precipitated salt saturation data - if constexpr (enableSaltPrecipitation) - saltpData = fp.get_double("SALTP"); - - // calculate the initial fluid states - for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { - auto& dofFluidState = initialFluidStates_[dofIdx]; - - dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx)); - - ////// - // set temperature - ////// - Scalar temperatureLoc = tempiData[dofIdx]; - if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0) - temperatureLoc = FluidSystem::surfaceTemperature; - dofFluidState.setTemperature(temperatureLoc); - - ////// - // set salt concentration - ////// - if constexpr (enableBrine) - dofFluidState.setSaltConcentration(saltData[dofIdx]); - - ////// - // set precipitated salt saturation - ////// - if constexpr (enableSaltPrecipitation) - dofFluidState.setSaltSaturation(saltpData[dofIdx]); - - ////// - // set saturations - ////// - if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) - dofFluidState.setSaturation(FluidSystem::waterPhaseIdx, - waterSaturationData[dofIdx]); - - if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){ - if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){ - dofFluidState.setSaturation(FluidSystem::gasPhaseIdx, - 1.0 - - waterSaturationData[dofIdx]); - } - else - dofFluidState.setSaturation(FluidSystem::gasPhaseIdx, - gasSaturationData[dofIdx]); - } - if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) - dofFluidState.setSaturation(FluidSystem::oilPhaseIdx, - 1.0 - - waterSaturationData[dofIdx] - - gasSaturationData[dofIdx]); - - ////// - // set phase pressures - ////// - Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase) - - // this assumes that capillary pressures only depend on the phase saturations - // and possibly on temperature. (this is always the case for ECL problems.) - std::array pc = {0}; - const auto& matParams = materialLawParams(dofIdx); - MaterialLaw::capillaryPressures(pc, matParams, dofFluidState); - Valgrind::CheckDefined(pressure); - Valgrind::CheckDefined(pc); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - if (Indices::oilEnabled) - dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx])); - else if (Indices::gasEnabled) - dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx])); - else if (Indices::waterEnabled) - //single (water) phase - dofFluidState.setPressure(phaseIdx, pressure); - } - - if (FluidSystem::enableDissolvedGas()) - dofFluidState.setRs(rsData[dofIdx]); - else if (Indices::gasEnabled && Indices::oilEnabled) - dofFluidState.setRs(0.0); - - if (FluidSystem::enableVaporizedOil()) - dofFluidState.setRv(rvData[dofIdx]); - else if (Indices::gasEnabled && Indices::oilEnabled) - dofFluidState.setRv(0.0); - - if (FluidSystem::enableVaporizedWater()) - dofFluidState.setRvw(rvwData[dofIdx]); - - ////// - // set invB_ - ////// - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx)); - dofFluidState.setInvB(phaseIdx, b); - - const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx)); - dofFluidState.setDensity(phaseIdx, rho); - - } - } - } + virtual void readEquilInitialCondition_() = 0; + virtual void readExplicitInitialCondition_() = 0; // update the hysteresis parameters of the material laws for the whole grid bool updateHysteresis_() @@ -2595,28 +1437,6 @@ protected: return true; } - void updateMaxPolymerAdsorption_() - { - // we need to update the max polymer adsoption data for all elements - this->updateProperty_("FlowProblem::updateMaxPolymerAdsorption_() failed:", - [this](unsigned compressedDofIdx, const IntensiveQuantities& iq) - { - this->updateMaxPolymerAdsorption_(compressedDofIdx,iq); - }); - } - - bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq) - { - const Scalar pa = scalarValue(iq.polymerAdsorption()); - auto& mpa = this->polymer_.maxAdsorption; - if (mpa[compressedDofIdx] < pa) { - mpa[compressedDofIdx] = pa; - return true; - } else { - return false; - } - } - Scalar getRockCompTransMultVal(std::size_t dofIdx) const { if (this->rockCompTransMultVal_.empty()) @@ -2625,7 +1445,7 @@ protected: return this->rockCompTransMultVal_[dofIdx]; } -private: +protected: struct PffDofData_ { ConditionalStorage thermalHalfTransIn; @@ -2665,6 +1485,8 @@ private: pffDofData_.update(distFn); } + virtual void updateExplicitQuantities_(int episodeIdx, int timeStepSize, bool first_step_after_restart) = 0; + void readBoundaryConditions_() { const auto& simulator = this->simulator(); @@ -2754,39 +1576,6 @@ private: return dtNext; } - void computeAndSetEqWeights_() - { - std::vector sumInvB(numPhases, 0.0); - const auto& gridView = this->gridView(); - ElementContext elemCtx(this->simulator()); - for(const auto& elem: elements(gridView, Dune::Partitions::interior)) { - elemCtx.updatePrimaryStencil(elem); - int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); - const auto& dofFluidState = initialFluidStates_[elemIdx]; - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx); - } - } - - std::size_t numDof = this->model().numGridDof(); - const auto& comm = this->simulator().vanguard().grid().comm(); - comm.sum(sumInvB.data(),sumInvB.size()); - Scalar numTotalDof = comm.sum(numDof); - - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - if (!FluidSystem::phaseIsActive(phaseIdx)) - continue; - - Scalar avgB = numTotalDof / sumInvB[phaseIdx]; - unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); - unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); - this->model().setEqWeight(activeSolventCompIdx, avgB); - } - } - int refPressurePhaseIdx_() const { if (FluidSystem::phaseIsActive(oilPhaseIdx)) { return oilPhaseIdx; @@ -2845,7 +1634,7 @@ private: // water compaction assert(!this->rockCompTransMultWc_.empty()); LhsEval SwMax = max(decay(fs.saturation(waterPhaseIdx)), this->maxWaterSaturation_[elementIdx]); - LhsEval SwDeltaMax = SwMax - initialFluidStates_[elementIdx].saturation(waterPhaseIdx); + LhsEval SwDeltaMax = SwMax - asImp_().initialFluidStates()[elementIdx].saturation(waterPhaseIdx); return this->rockCompTransMultWc_[tableIdx].eval(effectivePressure, SwDeltaMax, /*extrapolation=*/true); } @@ -2855,30 +1644,17 @@ private: std::shared_ptr materialLawManager_; std::shared_ptr thermalLawManager_; - FlowThresholdPressure thresholdPressures_; - - std::vector initialFluidStates_; - bool enableDriftCompensation_; GlobalEqVector drift_; WellModel wellModel_; AquiferModel aquiferModel_; - bool enableEclOutput_; bool enableVtkOutput_; - std::unique_ptr eclWriter_; - -#if HAVE_DAMARIS - bool enableDamarisOutput_ = false ; - std::unique_ptr damarisWriter_; -#endif + PffGridVector pffDofData_; TracerModel tracerModel_; - MixingRateControls mixControls_; - - ActionHandler actionHandler_; template struct BCData @@ -2909,12 +1685,14 @@ private: } }; + virtual void handleSolventBC(const BCProp::BCFace&, RateVector&) const = 0; + + virtual void handlePolymerBC(const BCProp::BCFace&, RateVector&) const = 0; + BCData bcindex_; bool nonTrivialBoundaryConditions_ = false; bool explicitRockCompaction_ = false; bool first_step_ = true; - ModuleParams moduleParams_; - }; diff --git a/opm/simulators/flow/FlowProblemBlackoil.hpp b/opm/simulators/flow/FlowProblemBlackoil.hpp new file mode 100644 index 000000000..862381bf7 --- /dev/null +++ b/opm/simulators/flow/FlowProblemBlackoil.hpp @@ -0,0 +1,1465 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + Copyright 2023 INRIA + Copyright 2024 SINTEF Digital + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::FlowProblemBlackoil + */ +#ifndef OPM_FLOW_PROBLEM_BLACK_HPP +#define OPM_FLOW_PROBLEM_BLACK_HPP + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include + +#include + +#include +#if HAVE_DAMARIS +#include +#endif + +#include +#include +#include +#include +#include + +namespace Opm { + +/*! + * \ingroup BlackOilSimulator + * + * \brief This problem simulates an input file given in the data format used by the + * commercial ECLiPSE simulator. + */ +template +class FlowProblemBlackoil : public FlowProblem +{ + // TODO: the naming of the Types might be able to be adjusted + using FlowProblemType = FlowProblem; + + using typename FlowProblemType::Scalar; + using typename FlowProblemType::Simulator; + using typename FlowProblemType::GridView; + using typename FlowProblemType::FluidSystem; + using typename FlowProblemType::Vanguard; + using typename FlowProblemType::GlobalEqVector; + using typename FlowProblemType::EqVector; + using FlowProblemType::dim; + using FlowProblemType::dimWorld; + using FlowProblemType::numEq; + using FlowProblemType::numPhases; + using FlowProblemType::numComponents; + + // TODO: potentially some cleaning up depending on the usage later here + using FlowProblemType::enableConvectiveMixing; + using FlowProblemType::enableBrine; + using FlowProblemType::enableDiffusion; + using FlowProblemType::enableDispersion; + using FlowProblemType::enableEnergy; + using FlowProblemType::enableExperiments; + using FlowProblemType::enableExtbo; + using FlowProblemType::enableFoam; + using FlowProblemType::enableMICP; + using FlowProblemType::enablePolymer; + using FlowProblemType::enablePolymerMolarWeight; + using FlowProblemType::enableSaltPrecipitation; + using FlowProblemType::enableSolvent; + using FlowProblemType::enableTemperature; + using FlowProblemType::enableThermalFluxBoundaries; + + using FlowProblemType::gasPhaseIdx; + using FlowProblemType::oilPhaseIdx; + using FlowProblemType::waterPhaseIdx; + + using FlowProblemType::waterCompIdx; + using FlowProblemType::oilCompIdx; + using FlowProblemType::gasCompIdx; + + using BoundaryRateVector = GetPropType; + using typename FlowProblemType::RateVector; + using typename FlowProblemType::PrimaryVariables; + using typename FlowProblemType::Indices; + using typename FlowProblemType::IntensiveQuantities; + using typename FlowProblemType::ElementContext; + + using typename FlowProblemType::MaterialLaw; + using typename FlowProblemType::DimMatrix; + + using SolventModule = BlackOilSolventModule; + using PolymerModule = BlackOilPolymerModule; + using FoamModule = BlackOilFoamModule; + using BrineModule = BlackOilBrineModule; + using ExtboModule = BlackOilExtboModule; + using MICPModule = BlackOilMICPModule; + using DispersionModule = BlackOilDispersionModule; + using DiffusionModule = BlackOilDiffusionModule; + using ConvectiveMixingModule = BlackOilConvectiveMixingModule; + using ModuleParams = typename BlackOilLocalResidualTPFA::ModuleParams; + + using InitialFluidState = typename EquilInitializer::ScalarFluidState; + using EclWriterType = EclWriter; +#if HAVE_DAMARIS + using DamarisWriterType = DamarisWriter; +#endif + + +public: + using FlowProblemType::porosity; + using FlowProblemType::pvtRegionIndex; + + /*! + * \copydoc FvBaseProblem::registerParameters + */ + static void registerParameters() + { + FlowProblemType::registerParameters(); + + EclWriterType::registerParameters(); +#if HAVE_DAMARIS + DamarisWriterType::registerParameters(); +#endif + VtkTracerModule::registerParameters(); + } + + /*! + * \copydoc Doxygen::defaultProblemConstructor + */ + explicit FlowProblemBlackoil(Simulator& simulator) + : FlowProblemType(simulator) + , thresholdPressures_(simulator) + , mixControls_(simulator.vanguard().schedule()) + , actionHandler_(simulator.vanguard().eclState(), + simulator.vanguard().schedule(), + simulator.vanguard().actionState(), + simulator.vanguard().summaryState(), + this->wellModel_, + simulator.vanguard().grid().comm()) + { + this->model().addOutputModule(new VtkTracerModule(simulator)); + + // Tell the black-oil extensions to initialize their internal data structures + const auto& vanguard = simulator.vanguard(); + + BlackOilBrineParams brineParams; + brineParams.template initFromState(vanguard.eclState()); + BrineModule::setParams(std::move(brineParams)); + + DiffusionModule::initFromState(vanguard.eclState()); + DispersionModule::initFromState(vanguard.eclState()); + + BlackOilExtboParams extboParams; + extboParams.template initFromState(vanguard.eclState()); + ExtboModule::setParams(std::move(extboParams)); + + BlackOilFoamParams foamParams; + foamParams.template initFromState(vanguard.eclState()); + FoamModule::setParams(std::move(foamParams)); + + BlackOilMICPParams micpParams; + micpParams.template initFromState(vanguard.eclState()); + MICPModule::setParams(std::move(micpParams)); + + BlackOilPolymerParams polymerParams; + polymerParams.template initFromState(vanguard.eclState()); + PolymerModule::setParams(std::move(polymerParams)); + + BlackOilSolventParams solventParams; + solventParams.template initFromState(vanguard.eclState(), vanguard.schedule()); + SolventModule::setParams(std::move(solventParams)); + + // create the ECL writer + eclWriter_ = std::make_unique(simulator); + enableEclOutput_ = Parameters::Get(); +#if HAVE_DAMARIS + // create Damaris writer + damarisWriter_ = std::make_unique(simulator); + enableDamarisOutput_ = Parameters::Get(); +#endif + } + + /*! + * \brief Called by the simulator before an episode begins. + */ + void beginEpisode() override + { + FlowProblemType::beginEpisode(); + + auto& simulator = this->simulator(); + int episodeIdx = simulator.episodeIndex(); + const auto& schedule = simulator.vanguard().schedule(); + + // Evaluate UDQ assign statements to make sure the settings are + // available as UDA controls for the current report step. + actionHandler_.evalUDQAssignments(episodeIdx, simulator.vanguard().udqState()); + + if (episodeIdx >= 0) { + const auto& oilVap = schedule[episodeIdx].oilvap(); + if (oilVap.getType() == OilVaporizationProperties::OilVaporization::VAPPARS) { + FluidSystem::setVapPars(oilVap.vap1(), oilVap.vap2()); + } else { + FluidSystem::setVapPars(0.0, 0.0); + } + } + + ConvectiveMixingModule::beginEpisode(simulator.vanguard().eclState(), simulator.vanguard().schedule(), episodeIdx, moduleParams_.convectiveMixingModuleParam); + } + + /*! + * \copydoc FvBaseProblem::finishInit + */ + void finishInit() + { + // TODO: there should be room to remove duplication for this function, + // but there is relatively complicated logic in the function calls in this function + // some refactoring is needed for this function + FlowProblemType::finishInit(); + auto& simulator = this->simulator(); + + auto finishTransmissibilities = [updated = false, this]() mutable + { + if (updated) { return; } + this->transmissibilities_.finishInit([&vg = this->simulator().vanguard()](const unsigned int it) { + return vg.gridIdxToEquilGridIdx(it); + }); + updated = true; + }; + + // calculating the TRANX, TRANY, TRANZ and NNC for output purpose + // for parallel running, it is based on global trans_ + // for serial running, it is based on the transmissibilities_ + // we try to avoid for the parallel running, has both global trans_ and transmissibilities_ allocated at the same time + if (enableEclOutput_) { + if (simulator.vanguard().grid().comm().size() > 1) { + if (simulator.vanguard().grid().comm().rank() == 0) + eclWriter_->setTransmissibilities(&simulator.vanguard().globalTransmissibility()); + } else { + finishTransmissibilities(); + eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities()); + } + + std::function equilGridToGrid = [&simulator](unsigned int i) { + return simulator.vanguard().gridEquilIdxToGridIdx(i); + }; + eclWriter_->extractOutputTransAndNNC(equilGridToGrid); + } + simulator.vanguard().releaseGlobalTransmissibilities(); + + const auto& eclState = simulator.vanguard().eclState(); + const auto& schedule = simulator.vanguard().schedule(); + + // Set the start time of the simulation + simulator.setStartTime(schedule.getStartTime()); + simulator.setEndTime(schedule.simTime(schedule.size() - 1)); + + // We want the episode index to be the same as the report step index to make + // things simpler, so we have to set the episode index to -1 because it is + // incremented by endEpisode(). The size of the initial time step and + // length of the initial episode is set to zero for the same reason. + simulator.setEpisodeIndex(-1); + simulator.setEpisodeLength(0.0); + + // the "NOGRAV" keyword from Frontsim or setting the EnableGravity to false + // disables gravity, else the standard value of the gravity constant at sea level + // on earth is used + this->gravity_ = 0.0; + if (Parameters::Get()) + this->gravity_[dim - 1] = 9.80665; + if (!eclState.getInitConfig().hasGravity()) + this->gravity_[dim - 1] = 0.0; + + if (this->enableTuning_) { + // if support for the TUNING keyword is enabled, we get the initial time + // steping parameters from it instead of from command line parameters + const auto& tuning = schedule[0].tuning(); + this->initialTimeStepSize_ = tuning.TSINIT.has_value() ? tuning.TSINIT.value() : -1.0; + this->maxTimeStepAfterWellEvent_ = tuning.TMAXWC; + } + + this->initFluidSystem_(); + + + + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && + FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) { + this->maxOilSaturation_.resize(this->model().numGridDof(), 0.0); + } + + this->readRockParameters_(simulator.vanguard().cellCenterDepths(), + [&simulator](const unsigned idx) + { + std::array coords; + simulator.vanguard().cartesianCoordinate(idx, coords); + for (auto& c : coords) { + ++c; + } + return coords; + }); + this->readMaterialParameters_(); + this->readThermalParameters_(); + + // write the static output files (EGRID, INIT) + if (enableEclOutput_) { + eclWriter_->writeInit(); + } + + finishTransmissibilities(); + + const auto& initconfig = eclState.getInitConfig(); + this->tracerModel_.init(initconfig.restartRequested()); + if (initconfig.restartRequested()) + readEclRestartSolution_(); + else + readInitialCondition_(); + + this->tracerModel_.prepareTracerBatches(); + + this->updatePffDofData_(); + + if constexpr (getPropValue()) { + const auto& vanguard = this->simulator().vanguard(); + const auto& gridView = vanguard.gridView(); + int numElements = gridView.size(/*codim=*/0); + this->polymer_.maxAdsorption.resize(numElements, 0.0); + } + + this->readBoundaryConditions_(); + + // compute and set eq weights based on initial b values + computeAndSetEqWeights_(); + + if (this->enableDriftCompensation_) { + this->drift_.resize(this->model().numGridDof()); + this->drift_ = 0.0; + } + + if (this->enableVtkOutput_ && eclState.getIOConfig().initOnly()) { + simulator.setTimeStepSize(0.0); + FlowProblemType::writeOutput(true); + } + + // after finishing the initialization and writing the initial solution, we move + // to the first "real" episode/report step + // for restart the episode index and start is already set + if (!initconfig.restartRequested()) { + simulator.startNextEpisode(schedule.seconds(1)); + simulator.setEpisodeIndex(0); + simulator.setTimeStepIndex(0); + } + + // TODO: move to the end for later refactoring of the function finishInit() + // deal with DRSDT + this->mixControls_.init(this->model().numGridDof(), + this->episodeIndex(), + eclState.runspec().tabdims().getNumPVTTables()); + } + + /*! + * \brief Called by the simulator after each time integration. + */ + void endTimeStep () override + { + FlowProblemType::endTimeStep(); + + const bool isSubStep = !this->simulator().episodeWillBeOver(); + + // For CpGrid with LGRs, ecl/vtk output is not supported yet. + const auto& grid = this->simulator().vanguard().gridView().grid(); + + using GridType = std::remove_cv_t>; + constexpr bool isCpGrid = std::is_same_v; + if (!isCpGrid || (grid.maxLevel() == 0)) { + this->eclWriter_->evalSummaryState(isSubStep); + } + + { + OPM_TIMEBLOCK(applyActions); + + const int episodeIdx = this->episodeIndex(); + auto& simulator = this->simulator(); + + // Re-ordering in case of Alugrid + this->actionHandler_ + .applyActions(episodeIdx, simulator.time() + simulator.timeStepSize(), + [this](const bool global) + { + using TransUpdateQuantities = typename Vanguard::TransmissibilityType::TransUpdateQuantities; + this->transmissibilities_ + .update(global, TransUpdateQuantities::All, [&vg = this->simulator().vanguard()] + (const unsigned int i) + { + return vg.gridIdxToEquilGridIdx(i); + }); + }); + } + + // Deal with "clogging" for the MICP model + if constexpr (enableMICP) { + auto& model = this->model(); + const auto& residual = model.linearizer().residual(); + + for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); ++globalDofIdx) { + auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx]; + MICPModule::checkCloggingMICP(model, phi, globalDofIdx); + } + } + + } + /*! + * \brief Called by the simulator after the end of an episode. + */ + void endEpisode() override + { + OPM_TIMEBLOCK(endEpisode); + const int episodeIdx = this->episodeIndex(); + // Rerun UDQ assignents following action processing on the final + // time step of this episode to make sure that any UDQ ASSIGN + // operations triggered in action blocks take effect. This is + // mainly to work around a shortcoming of the ScheduleState copy + // constructor which clears pending UDQ assignments under the + // assumption that all such assignments have been processed. If an + // action block happens to trigger on the final time step of an + // episode and that action block runs a UDQ assignment, then that + // assignment would be dropped and the rest of the simulator will + // never see its effect without this hack. + this->actionHandler_ + .evalUDQAssignments(episodeIdx, this->simulator().vanguard().udqState()); + + FlowProblemType::endEpisode(); + } + + /*! + * \brief Write the requested quantities of the current solution into the output + * files. + */ + void writeOutput(const SimulatorTimer& timer, bool verbose = true) + { + FlowProblemType::writeOutput(verbose); + + bool isSubStep = !this->simulator().episodeWillBeOver(); + + data::Solution localCellData = {}; +#if HAVE_DAMARIS + // N.B. the Damaris output has to be done before the ECL output as the ECL one + // does all kinds of std::move() relocation of data + if (enableDamarisOutput_) { + damarisWriter_->writeOutput(localCellData, isSubStep) ; + } +#endif + if (enableEclOutput_){ + eclWriter_->writeOutput(std::move(localCellData), timer, isSubStep); + } + } + + void finalizeOutput() + { + OPM_TIMEBLOCK(finalizeOutput); + // this will write all pending output to disk + // to avoid corruption of output files + eclWriter_.reset(); + } + + + /*! + * \copydoc FvBaseProblem::initialSolutionApplied() + */ + void initialSolutionApplied() override + { + FlowProblemType::initialSolutionApplied(); + + // let the object for threshold pressures initialize itself. this is done only at + // this point, because determining the threshold pressures may require to access + // the initial solution. + this->thresholdPressures_.finishInit(); + + if (this->simulator().episodeIndex() == 0) { + eclWriter_->writeInitialFIPReport(); + } + } + + void addToSourceDense(RateVector& rate, + unsigned globalDofIdx, + unsigned timeIdx) const override + { + this->aquiferModel_.addToSource(rate, globalDofIdx, timeIdx); + + // Add source term from deck + const auto& source = this->simulator().vanguard().schedule()[this->episodeIndex()].source(); + std::array ijk; + this->simulator().vanguard().cartesianCoordinate(globalDofIdx, ijk); + + if (source.hasSource(ijk)) { + const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx); + static std::array sc_map = {SourceComponent::WATER, SourceComponent::OIL, SourceComponent::GAS}; + static std::array phidx_map = {FluidSystem::waterPhaseIdx, FluidSystem::oilPhaseIdx, FluidSystem::gasPhaseIdx}; + static std::array cidx_map = {waterCompIdx, oilCompIdx, gasCompIdx}; + + for (unsigned i = 0; i < phidx_map.size(); ++i) { + const auto phaseIdx = phidx_map[i]; + const auto sourceComp = sc_map[i]; + const auto compIdx = cidx_map[i]; + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + Scalar mass_rate = source.rate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx); + if constexpr (getPropValue()) { + mass_rate /= FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx); + } + rate[Indices::canonicalToActiveComponentIndex(compIdx)] += mass_rate; + } + + if constexpr (enableSolvent) { + Scalar mass_rate = source.rate({ijk, SourceComponent::SOLVENT}) / this->model().dofTotalVolume(globalDofIdx); + if constexpr (getPropValue()) { + const auto& solventPvt = SolventModule::solventPvt(); + mass_rate /= solventPvt.referenceDensity(pvtRegionIdx); + } + rate[Indices::contiSolventEqIdx] += mass_rate; + } + if constexpr (enablePolymer) { + rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx); + } + if constexpr (enableEnergy) { + for (unsigned i = 0; i < phidx_map.size(); ++i) { + const auto phaseIdx = phidx_map[i]; + if (!FluidSystem::phaseIsActive(phaseIdx)) { + continue; + } + const auto sourceComp = sc_map[i]; + if (source.hasHrate({ijk, sourceComp})) { + rate[Indices::contiEnergyEqIdx] += source.hrate({ijk, sourceComp}) / this->model().dofTotalVolume(globalDofIdx); + } else { + const auto& intQuants = this->simulator().model().intensiveQuantities(globalDofIdx, /*timeIdx*/ 0); + auto fs = intQuants.fluidState(); + // if temperature is not set, use cell temperature as default + if (source.hasTemperature({ijk, sourceComp})) { + Scalar temperature = source.temperature({ijk, sourceComp}); + fs.setTemperature(temperature); + } + const auto& h = FluidSystem::enthalpy(fs, phaseIdx, pvtRegionIdx); + Scalar mass_rate = source.rate({ijk, sourceComp})/ this->model().dofTotalVolume(globalDofIdx); + Scalar energy_rate = getValue(h)*mass_rate; + rate[Indices::contiEnergyEqIdx] += energy_rate; + } + } + } + } + + // if requested, compensate systematic mass loss for cells which were "well + // behaved" in the last time step + if (this->enableDriftCompensation_) { + const auto& simulator = this->simulator(); + const auto& model = this->model(); + + // we use a lower tolerance for the compensation too + // assure the added drift from the last step does not + // cause convergence issues on the current step + Scalar maxCompensation = model.newtonMethod().tolerance()/10; + Scalar poro = this->porosity(globalDofIdx, timeIdx); + Scalar dt = simulator.timeStepSize(); + EqVector dofDriftRate = this->drift_[globalDofIdx]; + dofDriftRate /= dt*model.dofTotalVolume(globalDofIdx); + + // restrict drift compensation to the CNV tolerance + for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) { + Scalar cnv = std::abs(dofDriftRate[eqIdx])*dt*model.eqWeight(globalDofIdx, eqIdx)/poro; + if (cnv > maxCompensation) { + dofDriftRate[eqIdx] *= maxCompensation/cnv; + } + } + + for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) + rate[eqIdx] -= dofDriftRate[eqIdx]; + } + } + + /*! + * \brief Calculate the transmissibility multiplier due to porosity reduction. + * + * TODO: The API of this is a bit ad-hoc, it would be better to use context objects. + */ + template + LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants) const + { + OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier); + if (!enableSaltPrecipitation) + return 1.0; + + const auto& fs = intQuants.fluidState(); + unsigned tableIdx = fs.pvtRegionIndex(); + LhsEval porosityFactor = decay(1. - fs.saltSaturation()); + porosityFactor = min(porosityFactor, 1.0); + const auto& permfactTable = BrineModule::permfactTable(tableIdx); + return permfactTable.eval(porosityFactor, /*extrapolation=*/true); + } + + // temporary solution to facilitate output of initial state from flow + const InitialFluidState& initialFluidState(unsigned globalDofIdx) const + { return initialFluidStates_[globalDofIdx]; } + + std::vector& initialFluidStates() + { return initialFluidStates_; } + + const std::vector& initialFluidStates() const + { return initialFluidStates_; } + + const EclipseIO& eclIO() const + { return eclWriter_->eclIO(); } + + void setSubStepReport(const SimulatorReportSingle& report) + { return eclWriter_->setSubStepReport(report); } + + void setSimulationReport(const SimulatorReport& report) + { return eclWriter_->setSimulationReport(report); } + + InitialFluidState boundaryFluidState(unsigned globalDofIdx, const int directionId) const + { + OPM_TIMEBLOCK_LOCAL(boundaryFluidState); + const auto& bcprop = this->simulator().vanguard().schedule()[this->episodeIndex()].bcprop; + if (bcprop.size() > 0) { + FaceDir::DirEnum dir = FaceDir::FromIntersectionIndex(directionId); + + // index == 0: no boundary conditions for this + // global cell and direction + if (this->bcindex_(dir)[globalDofIdx] == 0) + return initialFluidStates_[globalDofIdx]; + + const auto& bc = bcprop[this->bcindex_(dir)[globalDofIdx]]; + if (bc.bctype == BCType::DIRICHLET ) + { + InitialFluidState fluidState; + const int pvtRegionIdx = this->pvtRegionIndex(globalDofIdx); + fluidState.setPvtRegionIndex(pvtRegionIdx); + + switch (bc.component) { + case BCComponent::OIL: + if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) + throw std::logic_error("oil is not active and you're trying to add oil BC"); + + fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0); + break; + case BCComponent::GAS: + if (!FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) + throw std::logic_error("gas is not active and you're trying to add gas BC"); + + fluidState.setSaturation(FluidSystem::gasPhaseIdx, 1.0); + break; + case BCComponent::WATER: + if (!FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) + throw std::logic_error("water is not active and you're trying to add water BC"); + + fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0); + break; + case BCComponent::SOLVENT: + case BCComponent::POLYMER: + case BCComponent::NONE: + throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC"); + } + fluidState.setTotalSaturation(1.0); + double pressure = initialFluidStates_[globalDofIdx].pressure(this->refPressurePhaseIdx_()); + const auto pressure_input = bc.pressure; + if (pressure_input) { + pressure = *pressure_input; + } + + std::array pc = {0}; + const auto& matParams = this->materialLawParams(globalDofIdx); + MaterialLaw::capillaryPressures(pc, matParams, fluidState); + Valgrind::CheckDefined(pressure); + Valgrind::CheckDefined(pc); + for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) { + const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx); + + fluidState.setPc(phaseIdx, pc[phaseIdx]); + if (Indices::oilEnabled) + fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx])); + else if (Indices::gasEnabled) + fluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx])); + else if (Indices::waterEnabled) + //single (water) phase + fluidState.setPressure(phaseIdx, pressure); + } + + double temperature = initialFluidStates_[globalDofIdx].temperature(0); // we only have one temperature + const auto temperature_input = bc.temperature; + if(temperature_input) + temperature = *temperature_input; + fluidState.setTemperature(temperature); + + if (FluidSystem::enableDissolvedGas()) { + fluidState.setRs(0.0); + fluidState.setRv(0.0); + } + if (FluidSystem::enableDissolvedGasInWater()) { + fluidState.setRsw(0.0); + } + if (FluidSystem::enableVaporizedWater()) + fluidState.setRvw(0.0); + + for (unsigned activePhaseIdx = 0; activePhaseIdx < FluidSystem::numActivePhases(); ++activePhaseIdx) { + const auto phaseIdx = FluidSystem::activeToCanonicalPhaseIdx(activePhaseIdx); + + const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState, phaseIdx, pvtRegionIdx); + fluidState.setInvB(phaseIdx, b); + + const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx); + fluidState.setDensity(phaseIdx, rho); + if (enableEnergy) { + const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx); + fluidState.setEnthalpy(phaseIdx, h); + } + } + fluidState.checkDefined(); + return fluidState; + } + } + return initialFluidStates_[globalDofIdx]; + } + + + const std::unique_ptr& eclWriter() const + { + return eclWriter_; + } + + void setConvData(const std::vector>& data) + { + eclWriter_->mutableOutputModule().setCnvData(data); + } + + /*! + * \brief Returns the maximum value of the gas dissolution factor at the current time + * for a given degree of freedom. + */ + Scalar maxGasDissolutionFactor(unsigned timeIdx, unsigned globalDofIdx) const + { + return this->mixControls_.maxGasDissolutionFactor(timeIdx, globalDofIdx, + this->episodeIndex(), + this->pvtRegionIndex(globalDofIdx)); + } + + /*! + * \brief Returns the maximum value of the oil vaporization factor at the current + * time for a given degree of freedom. + */ + Scalar maxOilVaporizationFactor(unsigned timeIdx, unsigned globalDofIdx) const + { + return this->mixControls_.maxOilVaporizationFactor(timeIdx, globalDofIdx, + this->episodeIndex(), + this->pvtRegionIndex(globalDofIdx)); + } + + /*! + * \brief Return if the storage term of the first iteration is identical to the storage + * term for the solution of the previous time step. + * + * For quite technical reasons, the storage term cannot be recycled if either DRSDT + * or DRVDT are active. Nor if the porosity is changes between timesteps + * using a pore volume multiplier (i.e., poreVolumeMultiplier() != 1.0) + */ + bool recycleFirstIterationStorage() const + { + int episodeIdx = this->episodeIndex(); + return !this->mixControls_.drsdtActive(episodeIdx) && + !this->mixControls_.drvdtActive(episodeIdx) && + this->rockCompPoroMultWc_.empty() && + this->rockCompPoroMult_.empty(); + } + + /*! + * \copydoc FvBaseProblem::initial + * + * The reservoir problem uses a constant boundary condition for + * the whole domain. + */ + template + void initial(PrimaryVariables& values, const Context& context, unsigned spaceIdx, unsigned timeIdx) const + { + unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx); + + values.setPvtRegionIndex(pvtRegionIndex(context, spaceIdx, timeIdx)); + values.assignNaive(initialFluidStates_[globalDofIdx]); + + SolventModule::assignPrimaryVars(values, + enableSolvent ? this->solventSaturation_[globalDofIdx] : 0.0, + enableSolvent ? this->solventRsw_[globalDofIdx] : 0.0); + + if constexpr (enablePolymer) + values[Indices::polymerConcentrationIdx] = this->polymer_.concentration[globalDofIdx]; + + if constexpr (enablePolymerMolarWeight) + values[Indices::polymerMoleWeightIdx]= this->polymer_.moleWeight[globalDofIdx]; + + if constexpr (enableBrine) { + if (enableSaltPrecipitation && values.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) { + values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltSaturation(); + } + else { + values[Indices::saltConcentrationIdx] = initialFluidStates_[globalDofIdx].saltConcentration(); + } + } + + if constexpr (enableMICP){ + values[Indices::microbialConcentrationIdx] = this->micp_.microbialConcentration[globalDofIdx]; + values[Indices::oxygenConcentrationIdx]= this->micp_.oxygenConcentration[globalDofIdx]; + values[Indices::ureaConcentrationIdx]= this->micp_.ureaConcentration[globalDofIdx]; + values[Indices::calciteConcentrationIdx]= this->micp_.calciteConcentration[globalDofIdx]; + values[Indices::biofilmConcentrationIdx]= this->micp_.biofilmConcentration[globalDofIdx]; + } + + values.checkDefined(); + } + + + Scalar drsdtcon(unsigned elemIdx, int episodeIdx) const + { + return this->mixControls_.drsdtcon(elemIdx, episodeIdx, + this->pvtRegionIndex(elemIdx)); + } + + /*! + * \copydoc FvBaseProblem::boundary + * + * Reservoir simulation uses no-flow conditions as default for all boundaries. + */ + template + void boundary(BoundaryRateVector& values, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) const + { + OPM_TIMEBLOCK_LOCAL(eclProblemBoundary); + if (!context.intersection(spaceIdx).boundary()) + return; + + if constexpr (!enableEnergy || !enableThermalFluxBoundaries) + values.setNoFlow(); + else { + // in the energy case we need to specify a non-trivial boundary condition + // because the geothermal gradient needs to be maintained. for this, we + // simply assume the initial temperature at the boundary and specify the + // thermal flow accordingly. in this context, "thermal flow" means energy + // flow due to a temerature gradient while assuming no-flow for mass + unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx); + unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx); + values.setThermalFlow(context, spaceIdx, timeIdx, this->initialFluidStates_[globalDofIdx] ); + } + + if (this->nonTrivialBoundaryConditions()) { + unsigned indexInInside = context.intersection(spaceIdx).indexInInside(); + unsigned interiorDofIdx = context.interiorScvIndex(spaceIdx, timeIdx); + unsigned globalDofIdx = context.globalSpaceIndex(interiorDofIdx, timeIdx); + unsigned pvtRegionIdx = pvtRegionIndex(context, spaceIdx, timeIdx); + const auto [type, massrate] = this->boundaryCondition(globalDofIdx, indexInInside); + if (type == BCType::THERMAL) + values.setThermalFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside)); + else if (type == BCType::FREE || type == BCType::DIRICHLET) + values.setFreeFlow(context, spaceIdx, timeIdx, this->boundaryFluidState(globalDofIdx, indexInInside)); + else if (type == BCType::RATE) + values.setMassRate(massrate, pvtRegionIdx); + } + } + + + /*! + * \copydoc BlackOilBaseProblem::thresholdPressure + */ + Scalar thresholdPressure(unsigned elem1Idx, unsigned elem2Idx) const + { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); } + + const FlowThresholdPressure& thresholdPressure() const + { return thresholdPressures_; } + + FlowThresholdPressure& thresholdPressure() + { return thresholdPressures_; } + + const ModuleParams& moduleParams() const + { + return moduleParams_; + } + + template + void serializeOp(Serializer& serializer) + { + serializer(static_cast(*this)); + serializer(mixControls_); + serializer(*eclWriter_); + } + +protected: + void updateExplicitQuantities_(int episodeIdx, int timeStepSize, const bool first_step_after_restart) override + { + this->updateExplicitQuantities_(first_step_after_restart); + + if constexpr (getPropValue()) + updateMaxPolymerAdsorption_(); + + mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize); + } + + void updateMaxPolymerAdsorption_() + { + // we need to update the max polymer adsoption data for all elements + this->updateProperty_("FlowProblemBlackoil::updateMaxPolymerAdsorption_() failed:", + [this](unsigned compressedDofIdx, const IntensiveQuantities& iq) + { + this->updateMaxPolymerAdsorption_(compressedDofIdx,iq); + }); + } + + bool updateMaxPolymerAdsorption_(unsigned compressedDofIdx, const IntensiveQuantities& iq) + { + const Scalar pa = scalarValue(iq.polymerAdsorption()); + auto& mpa = this->polymer_.maxAdsorption; + if (mpa[compressedDofIdx] < pa) { + mpa[compressedDofIdx] = pa; + return true; + } else { + return false; + } + } + + void computeAndSetEqWeights_() + { + std::vector sumInvB(numPhases, 0.0); + const auto& gridView = this->gridView(); + ElementContext elemCtx(this->simulator()); + for(const auto& elem: elements(gridView, Dune::Partitions::interior)) { + elemCtx.updatePrimaryStencil(elem); + int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& dofFluidState = this->initialFluidStates_[elemIdx]; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + sumInvB[phaseIdx] += dofFluidState.invB(phaseIdx); + } + } + + std::size_t numDof = this->model().numGridDof(); + const auto& comm = this->simulator().vanguard().grid().comm(); + comm.sum(sumInvB.data(),sumInvB.size()); + Scalar numTotalDof = comm.sum(numDof); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + Scalar avgB = numTotalDof / sumInvB[phaseIdx]; + unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx); + unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx); + this->model().setEqWeight(activeSolventCompIdx, avgB); + } + } + + // update the parameters needed for DRSDT and DRVDT + bool updateCompositionChangeLimits_() + { + OPM_TIMEBLOCK(updateCompositionChangeLimits); + // update the "last Rs" values for all elements, including the ones in the ghost + // and overlap regions + int episodeIdx = this->episodeIndex(); + std::array active{this->mixControls_.drsdtConvective(episodeIdx), + this->mixControls_.drsdtActive(episodeIdx), + this->mixControls_.drvdtActive(episodeIdx)}; + if (!active[0] && !active[1] && !active[2]) { + return false; + } + + this->updateProperty_("FlowProblemBlackoil::updateCompositionChangeLimits_()) failed:", + [this,episodeIdx,active](unsigned compressedDofIdx, + const IntensiveQuantities& iq) + { + const DimMatrix& perm = this->intrinsicPermeability(compressedDofIdx); + const Scalar distZ = active[0] ? this->simulator().vanguard().cellThickness(compressedDofIdx) : 0.0; + const int pvtRegionIdx = this->pvtRegionIndex(compressedDofIdx); + this->mixControls_.update(compressedDofIdx, + iq, + episodeIdx, + this->gravity_[dim - 1], + perm[dim - 1][dim - 1], + distZ, + pvtRegionIdx); + } + ); + + return true; + } + + + void readEclRestartSolution_() + { + // Throw an exception if the grid has LGRs. Refined grid are not supported for restart. + if(this->simulator().vanguard().grid().maxLevel() > 0) { + throw std::invalid_argument("Refined grids are not yet supported for restart "); + } + + // Set the start time of the simulation + auto& simulator = this->simulator(); + const auto& schedule = simulator.vanguard().schedule(); + const auto& eclState = simulator.vanguard().eclState(); + const auto& initconfig = eclState.getInitConfig(); + const int restart_step = initconfig.getRestartStep(); + { + simulator.setTime(schedule.seconds(restart_step)); + + simulator.startNextEpisode(simulator.startTime() + simulator.time(), + schedule.stepLength(restart_step)); + simulator.setEpisodeIndex(restart_step); + } + this->eclWriter_->beginRestart(); + + Scalar dt = std::min(this->eclWriter_->restartTimeStepSize(), simulator.episodeLength()); + simulator.setTimeStepSize(dt); + + std::size_t numElems = this->model().numGridDof(); + this->initialFluidStates_.resize(numElems); + if constexpr (enableSolvent) { + this->solventSaturation_.resize(numElems, 0.0); + this->solventRsw_.resize(numElems, 0.0); + } + + if constexpr (enablePolymer) + this->polymer_.concentration.resize(numElems, 0.0); + + if constexpr (enablePolymerMolarWeight) { + const std::string msg {"Support of the RESTART for polymer molecular weight " + "is not implemented yet. The polymer weight value will be " + "zero when RESTART begins"}; + OpmLog::warning("NO_POLYMW_RESTART", msg); + this->polymer_.moleWeight.resize(numElems, 0.0); + } + + if constexpr (enableMICP) { + this->micp_.resize(numElems); + } + + // Initialize mixing controls before trying to set any lastRx valuesx + this->mixControls_.init(numElems, restart_step, eclState.runspec().tabdims().getNumPVTTables()); + + for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { + auto& elemFluidState = this->initialFluidStates_[elemIdx]; + elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx)); + this->eclWriter_->outputModule().initHysteresisParams(simulator, elemIdx); + this->eclWriter_->outputModule().assignToFluidState(elemFluidState, elemIdx); + + // Note: Function processRestartSaturations_() mutates the + // 'ssol' argument--the value from the restart file--if solvent + // is enabled. Then, store the updated solvent saturation into + // 'solventSaturation_'. Otherwise, just pass a dummy value to + // the function and discard the unchanged result. Do not index + // into 'solventSaturation_' unless solvent is enabled. + { + auto ssol = enableSolvent + ? this->eclWriter_->outputModule().getSolventSaturation(elemIdx) + : Scalar(0); + + this->processRestartSaturations_(elemFluidState, ssol); + + if constexpr (enableSolvent) { + this->solventSaturation_[elemIdx] = ssol; + this->solventRsw_[elemIdx] = this->eclWriter_->outputModule().getSolventRsw(elemIdx); + } + } + + // For CO2STORE and H2STORE we need to set the initial temperature for isothermal simulations + bool isThermal = eclState.getSimulationConfig().isThermal(); + bool needTemperature = (eclState.runspec().co2Storage() || eclState.runspec().h2Storage()); + if (!isThermal && needTemperature) { + const auto& fp = simulator.vanguard().eclState().fieldProps(); + elemFluidState.setTemperature(fp.get_double("TEMPI")[elemIdx]); + } + + this->mixControls_.updateLastValues(elemIdx, elemFluidState.Rs(), elemFluidState.Rv()); + + if constexpr (enablePolymer) + this->polymer_.concentration[elemIdx] = this->eclWriter_->outputModule().getPolymerConcentration(elemIdx); + if constexpr (enableMICP){ + this->micp_.microbialConcentration[elemIdx] = this->eclWriter_->outputModule().getMicrobialConcentration(elemIdx); + this->micp_.oxygenConcentration[elemIdx] = this->eclWriter_->outputModule().getOxygenConcentration(elemIdx); + this->micp_.ureaConcentration[elemIdx] = this->eclWriter_->outputModule().getUreaConcentration(elemIdx); + this->micp_.biofilmConcentration[elemIdx] = this->eclWriter_->outputModule().getBiofilmConcentration(elemIdx); + this->micp_.calciteConcentration[elemIdx] = this->eclWriter_->outputModule().getCalciteConcentration(elemIdx); + } + // if we need to restart for polymer molecular weight simulation, we need to add related here + } + + const int episodeIdx = this->episodeIndex(); + this->mixControls_.updateMaxValues(episodeIdx, simulator.timeStepSize()); + + // assign the restart solution to the current solution. note that we still need + // to compute real initial solution after this because the initial fluid states + // need to be correct for stuff like boundary conditions. + auto& sol = this->model().solution(/*timeIdx=*/0); + const auto& gridView = this->gridView(); + ElementContext elemCtx(simulator); + for (const auto& elem : elements(gridView, Dune::Partitions::interior)) { + elemCtx.updatePrimaryStencil(elem); + int elemIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + this->initial(sol[elemIdx], elemCtx, /*spaceIdx=*/0, /*timeIdx=*/0); + } + + // make sure that the ghost and overlap entities exhibit the correct + // solution. alternatively, this could be done in the loop above by also + // considering non-interior elements. Since the initial() method might not work + // 100% correctly for such elements, let's play safe and explicitly synchronize + // using message passing. + this->model().syncOverlap(); + + this->eclWriter_->endRestart(); + } + + void readEquilInitialCondition_() + { + const auto& simulator = this->simulator(); + + // initial condition corresponds to hydrostatic conditions. + EquilInitializer equilInitializer(simulator, *(this->materialLawManager_)); + + std::size_t numElems = this->model().numGridDof(); + this->initialFluidStates_.resize(numElems); + for (std::size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) { + auto& elemFluidState = this->initialFluidStates_[elemIdx]; + elemFluidState.assign(equilInitializer.initialFluidState(elemIdx)); + } + } + + void readExplicitInitialCondition_() override + { + const auto& simulator = this->simulator(); + const auto& vanguard = simulator.vanguard(); + const auto& eclState = vanguard.eclState(); + const auto& fp = eclState.fieldProps(); + bool has_swat = fp.has_double("SWAT"); + bool has_sgas = fp.has_double("SGAS"); + bool has_rs = fp.has_double("RS"); + bool has_rv = fp.has_double("RV"); + bool has_rvw = fp.has_double("RVW"); + bool has_pressure = fp.has_double("PRESSURE"); + bool has_salt = fp.has_double("SALT"); + bool has_saltp = fp.has_double("SALTP"); + + // make sure all required quantities are enables + if (Indices::numPhases > 1) { + if (FluidSystem::phaseIsActive(waterPhaseIdx) && !has_swat) + throw std::runtime_error("The ECL input file requires the presence of the SWAT keyword if " + "the water phase is active"); + if (FluidSystem::phaseIsActive(gasPhaseIdx) && !has_sgas && FluidSystem::phaseIsActive(oilPhaseIdx)) + throw std::runtime_error("The ECL input file requires the presence of the SGAS keyword if " + "the gas phase is active"); + } + if (!has_pressure) + throw std::runtime_error("The ECL input file requires the presence of the PRESSURE " + "keyword if the model is initialized explicitly"); + if (FluidSystem::enableDissolvedGas() && !has_rs) + throw std::runtime_error("The ECL input file requires the RS keyword to be present if" + " dissolved gas is enabled"); + if (FluidSystem::enableVaporizedOil() && !has_rv) + throw std::runtime_error("The ECL input file requires the RV keyword to be present if" + " vaporized oil is enabled"); + if (FluidSystem::enableVaporizedWater() && !has_rvw) + throw std::runtime_error("The ECL input file requires the RVW keyword to be present if" + " vaporized water is enabled"); + if (enableBrine && !has_salt) + throw std::runtime_error("The ECL input file requires the SALT keyword to be present if" + " brine is enabled and the model is initialized explicitly"); + if (enableSaltPrecipitation && !has_saltp) + throw std::runtime_error("The ECL input file requires the SALTP keyword to be present if" + " salt precipitation is enabled and the model is initialized explicitly"); + + std::size_t numDof = this->model().numGridDof(); + + initialFluidStates_.resize(numDof); + + std::vector waterSaturationData; + std::vector gasSaturationData; + std::vector pressureData; + std::vector rsData; + std::vector rvData; + std::vector rvwData; + std::vector tempiData; + std::vector saltData; + std::vector saltpData; + + if (FluidSystem::phaseIsActive(waterPhaseIdx) && Indices::numPhases > 1) + waterSaturationData = fp.get_double("SWAT"); + else + waterSaturationData.resize(numDof); + + if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx)) + gasSaturationData = fp.get_double("SGAS"); + else + gasSaturationData.resize(numDof); + + pressureData = fp.get_double("PRESSURE"); + if (FluidSystem::enableDissolvedGas()) + rsData = fp.get_double("RS"); + + if (FluidSystem::enableVaporizedOil()) + rvData = fp.get_double("RV"); + + if (FluidSystem::enableVaporizedWater()) + rvwData = fp.get_double("RVW"); + + // initial reservoir temperature + tempiData = fp.get_double("TEMPI"); + + // initial salt concentration data + if constexpr (enableBrine) + saltData = fp.get_double("SALT"); + + // initial precipitated salt saturation data + if constexpr (enableSaltPrecipitation) + saltpData = fp.get_double("SALTP"); + + // calculate the initial fluid states + for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) { + auto& dofFluidState = initialFluidStates_[dofIdx]; + + dofFluidState.setPvtRegionIndex(pvtRegionIndex(dofIdx)); + + ////// + // set temperature + ////// + Scalar temperatureLoc = tempiData[dofIdx]; + if (!std::isfinite(temperatureLoc) || temperatureLoc <= 0) + temperatureLoc = FluidSystem::surfaceTemperature; + dofFluidState.setTemperature(temperatureLoc); + + ////// + // set salt concentration + ////// + if constexpr (enableBrine) + dofFluidState.setSaltConcentration(saltData[dofIdx]); + + ////// + // set precipitated salt saturation + ////// + if constexpr (enableSaltPrecipitation) + dofFluidState.setSaltSaturation(saltpData[dofIdx]); + + ////// + // set saturations + ////// + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) + dofFluidState.setSaturation(FluidSystem::waterPhaseIdx, + waterSaturationData[dofIdx]); + + if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){ + if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){ + dofFluidState.setSaturation(FluidSystem::gasPhaseIdx, + 1.0 + - waterSaturationData[dofIdx]); + } + else + dofFluidState.setSaturation(FluidSystem::gasPhaseIdx, + gasSaturationData[dofIdx]); + } + if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) + dofFluidState.setSaturation(FluidSystem::oilPhaseIdx, + 1.0 + - waterSaturationData[dofIdx] + - gasSaturationData[dofIdx]); + + ////// + // set phase pressures + ////// + Scalar pressure = pressureData[dofIdx]; // oil pressure (or gas pressure for water-gas system or water pressure for single phase) + + // this assumes that capillary pressures only depend on the phase saturations + // and possibly on temperature. (this is always the case for ECL problems.) + std::array pc = {0}; + const auto& matParams = this->materialLawParams(dofIdx); + MaterialLaw::capillaryPressures(pc, matParams, dofFluidState); + Valgrind::CheckDefined(pressure); + Valgrind::CheckDefined(pc); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + if (Indices::oilEnabled) + dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx])); + else if (Indices::gasEnabled) + dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx])); + else if (Indices::waterEnabled) + //single (water) phase + dofFluidState.setPressure(phaseIdx, pressure); + } + + if (FluidSystem::enableDissolvedGas()) + dofFluidState.setRs(rsData[dofIdx]); + else if (Indices::gasEnabled && Indices::oilEnabled) + dofFluidState.setRs(0.0); + + if (FluidSystem::enableVaporizedOil()) + dofFluidState.setRv(rvData[dofIdx]); + else if (Indices::gasEnabled && Indices::oilEnabled) + dofFluidState.setRv(0.0); + + if (FluidSystem::enableVaporizedWater()) + dofFluidState.setRvw(rvwData[dofIdx]); + + ////// + // set invB_ + ////// + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!FluidSystem::phaseIsActive(phaseIdx)) + continue; + + const auto& b = FluidSystem::inverseFormationVolumeFactor(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx)); + dofFluidState.setInvB(phaseIdx, b); + + const auto& rho = FluidSystem::density(dofFluidState, phaseIdx, pvtRegionIndex(dofIdx)); + dofFluidState.setDensity(phaseIdx, rho); + + } + } + } + + + void processRestartSaturations_(InitialFluidState& elemFluidState, Scalar& solventSaturation) + { + // each phase needs to be above certain value to be claimed to be existing + // this is used to recover some RESTART running with the defaulted single-precision format + const Scalar smallSaturationTolerance = 1.e-6; + Scalar sumSaturation = 0.0; + for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (FluidSystem::phaseIsActive(phaseIdx)) { + if (elemFluidState.saturation(phaseIdx) < smallSaturationTolerance) + elemFluidState.setSaturation(phaseIdx, 0.0); + + sumSaturation += elemFluidState.saturation(phaseIdx); + } + + } + if constexpr (enableSolvent) { + if (solventSaturation < smallSaturationTolerance) + solventSaturation = 0.0; + + sumSaturation += solventSaturation; + } + + assert(sumSaturation > 0.0); + + for (std::size_t phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (FluidSystem::phaseIsActive(phaseIdx)) { + const Scalar saturation = elemFluidState.saturation(phaseIdx) / sumSaturation; + elemFluidState.setSaturation(phaseIdx, saturation); + } + } + if constexpr (enableSolvent) { + solventSaturation = solventSaturation / sumSaturation; + } + } + + void readInitialCondition_() override + { + FlowProblemType::readInitialCondition_(); + + if constexpr (enableSolvent || enablePolymer || enablePolymerMolarWeight || enableMICP) + this->readBlackoilExtentionsInitialConditions_(this->model().numGridDof(), + enableSolvent, + enablePolymer, + enablePolymerMolarWeight, + enableMICP); + + } + + void handleSolventBC(const BCProp::BCFace& bc, RateVector& rate) const override + { + if constexpr (!enableSolvent) + throw std::logic_error("solvent is disabled and you're trying to add solvent to BC"); + + rate[Indices::solventSaturationIdx] = bc.rate; + } + + void handlePolymerBC(const BCProp::BCFace& bc, RateVector& rate) const override + { + if constexpr (!enablePolymer) + throw std::logic_error("polymer is disabled and you're trying to add polymer to BC"); + + rate[Indices::polymerConcentrationIdx] = bc.rate; + } + + void updateExplicitQuantities_(const bool first_step_after_restart) + { + OPM_TIMEBLOCK(updateExplicitQuantities); + const bool invalidateFromMaxWaterSat = this->updateMaxWaterSaturation_(); + const bool invalidateFromMinPressure = this->updateMinPressure_(); + + // update hysteresis and max oil saturation used in vappars + const bool invalidateFromHyst = this->updateHysteresis_(); + const bool invalidateFromMaxOilSat = this->updateMaxOilSaturation_(); + + // deal with DRSDT and DRVDT + const bool invalidateDRDT = !first_step_after_restart && this->updateCompositionChangeLimits_(); + + // the derivatives may have changed + const bool invalidateIntensiveQuantities + = invalidateFromMaxWaterSat || invalidateFromMinPressure || invalidateFromHyst || invalidateFromMaxOilSat || invalidateDRDT; + if (invalidateIntensiveQuantities) { + OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities); + this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } + + this->updateRockCompTransMultVal_(); + } + + FlowThresholdPressure thresholdPressures_; + + std::vector initialFluidStates_; + + bool enableEclOutput_; + std::unique_ptr eclWriter_; +#if HAVE_DAMARIS + bool enableDamarisOutput_ = false ; + std::unique_ptr damarisWriter_; +#endif + MixingRateControls mixControls_; + + ActionHandler actionHandler_; + + ModuleParams moduleParams_; +}; + +} // namespace Opm + +#endif // OPM_FLOW_PROBLEM_BLACK_HPP diff --git a/opm/simulators/flow/FlowProblemBlackoilProperties.hpp b/opm/simulators/flow/FlowProblemBlackoilProperties.hpp new file mode 100644 index 000000000..6b9e898aa --- /dev/null +++ b/opm/simulators/flow/FlowProblemBlackoilProperties.hpp @@ -0,0 +1,99 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 2 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \copydoc Opm::FlowBaseProblemBlackoil + */ +#ifndef OPM_FLOW_BASE_PROBLEM_BLACKOIL_PROPERTIES_HPP +#define OPM_FLOW_BASE_PROBLEM_BLACKOIL_PROPERTIES_HPP + + +#include + +#include + +#include +#include +#include +#include + + +#include + +namespace Opm { +template +class FlowProblemBlackoil; +} + +namespace Opm::Properties { + +namespace TTag { + +struct FlowBaseProblemBlackoil { + using InheritsFrom = std::tuple; +}; + +} + +// Set the problem property +template +struct Problem +{ using type = FlowProblemBlackoil; }; + +template +struct Model +{ using type = FIBlackOilModel; }; + +// Set the material law for fluid fluxes +template +struct MaterialLaw +{ +private: + using Scalar = GetPropType; + using FluidSystem = GetPropType; + + using Traits = ThreePhaseMaterialTraits; + +public: + using EclMaterialLawManager = ::Opm::EclMaterialLawManager; + + using type = typename EclMaterialLawManager::MaterialLaw; +}; + +// Use the "velocity module" which uses the Eclipse "NEWTRAN" transmissibilities +template +struct FluxModule +{ using type = NewTranFluxModule; }; + +// Use the dummy gradient calculator in order not to do unnecessary work. +template +struct GradientCalculator +{ using type = DummyGradientCalculator; }; + +} // namespace Opm::Properties + +#endif // OPM_FLOW_BASE_PROBLEM_BLACKOIL_PROPERTIES_HPP diff --git a/opm/simulators/flow/MixingRateControls.cpp b/opm/simulators/flow/MixingRateControls.cpp index 7ee2d8894..fd9b6b149 100644 --- a/opm/simulators/flow/MixingRateControls.cpp +++ b/opm/simulators/flow/MixingRateControls.cpp @@ -23,7 +23,7 @@ /*! * \file * - * \copydoc Opm::FlowProblem + * \copydoc Opm::FlowProblemBlackoil */ #include diff --git a/opm/simulators/flow/MixingRateControls.hpp b/opm/simulators/flow/MixingRateControls.hpp index e7ff42369..b575ea07c 100644 --- a/opm/simulators/flow/MixingRateControls.hpp +++ b/opm/simulators/flow/MixingRateControls.hpp @@ -23,7 +23,7 @@ /*! * \file * - * \copydoc Opm::FlowProblem + * \copydoc Opm::FlowProblemBlackoil */ #ifndef OPM_MIXING_RATE_CONTROLS_HPP #define OPM_MIXING_RATE_CONTROLS_HPP @@ -40,7 +40,7 @@ namespace Opm { class EclipseState; -//! \brief Class handling mixing rate controls for a FlowProblem. +//! \brief Class handling mixing rate controls for a FlowProblemBlackoil. template class MixingRateControls { diff --git a/opm/simulators/linalg/ISTLSolver.hpp b/opm/simulators/linalg/ISTLSolver.hpp index 85358d72c..b02bbd603 100644 --- a/opm/simulators/linalg/ISTLSolver.hpp +++ b/opm/simulators/linalg/ISTLSolver.hpp @@ -35,7 +35,7 @@ #include #include #include -#include +#include #include #include #include diff --git a/tests/TestTypeTag.hpp b/tests/TestTypeTag.hpp index 4fb69fc84..768a7659f 100644 --- a/tests/TestTypeTag.hpp +++ b/tests/TestTypeTag.hpp @@ -31,7 +31,7 @@ #include #include -#include +#include #include #include #include @@ -42,7 +42,7 @@ namespace TTag { struct TestTypeTag { - using InheritsFrom = std::tuple; }; @@ -51,7 +51,7 @@ struct TestTypeTag // Set the problem class template struct Problem { - using type = FlowProblem; + using type = FlowProblemBlackoil; }; // Enable experimental features for ebos: ebos is the research simulator of the OPM diff --git a/tests/test_RestartSerialization.cpp b/tests/test_RestartSerialization.cpp index c69a6b539..6989be8b0 100644 --- a/tests/test_RestartSerialization.cpp +++ b/tests/test_RestartSerialization.cpp @@ -38,7 +38,7 @@ #include #include -#include +#include #include #include #include diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 1971e5293..75c2284c9 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -37,7 +37,8 @@ #include #include -#include +#include +#include #include #include #include @@ -72,11 +73,11 @@ namespace TTag { struct TestEquilTypeTag { - using InheritsFrom = std::tuple; }; struct TestEquilVapwatTypeTag { - using InheritsFrom = std::tuple; }; } diff --git a/tests/test_glift1.cpp b/tests/test_glift1.cpp index 4c3a02dc8..62868d67b 100644 --- a/tests/test_glift1.cpp +++ b/tests/test_glift1.cpp @@ -34,7 +34,7 @@ #include #include #include -#include +#include #include #include #include @@ -123,7 +123,7 @@ BOOST_GLOBAL_FIXTURE(GliftFixture); BOOST_AUTO_TEST_CASE(G1) { - //using TypeTag = Opm::Properties::TTag::FlowProblem; + //using TypeTag = Opm::Properties::TTag::FlowProblemBlackoil; using TypeTag = Opm::Properties::TTag::TestGliftTypeTag; //using EclProblem = Opm::EclProblem; //using EclWellModel = typename EclProblem::EclWellModel; diff --git a/tests/test_wellmodel.cpp b/tests/test_wellmodel.cpp index da829a942..48c6d4697 100644 --- a/tests/test_wellmodel.cpp +++ b/tests/test_wellmodel.cpp @@ -47,7 +47,7 @@ #include #include #include -#include +#include #include