diff --git a/opm/simulators/flow/EclWriter.hpp b/opm/simulators/flow/EclWriter.hpp index ac40c95fc..d6675754b 100644 --- a/opm/simulators/flow/EclWriter.hpp +++ b/opm/simulators/flow/EclWriter.hpp @@ -482,24 +482,30 @@ public: void beginRestart() { - bool enableHysteresis = simulator_.problem().materialLawManager()->enableHysteresis(); + bool enablePCHysteresis = simulator_.problem().materialLawManager()->enablePCHysteresis(); + bool enableNonWettingHysteresis = simulator_.problem().materialLawManager()->enableNonWettingHysteresis(); + bool enableWettingHysteresis = simulator_.problem().materialLawManager()->enableWettingHysteresis(); + bool oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + bool gasActive = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + bool waterActive = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); bool enableSwatinit = simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT"); bool opm_rst_file = Parameters::get(); bool read_temp = enableEnergy || (opm_rst_file && enableTemperature); std::vector solutionKeys{ {"PRESSURE", UnitSystem::measure::pressure}, - {"SWAT", UnitSystem::measure::identity, static_cast(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))}, - {"SGAS", UnitSystem::measure::identity, static_cast(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))}, + {"SWAT", UnitSystem::measure::identity, static_cast(waterActive)}, + {"SGAS", UnitSystem::measure::identity, static_cast(gasActive)}, {"TEMP" , UnitSystem::measure::temperature, read_temp}, {"SSOLVENT" , UnitSystem::measure::identity, enableSolvent}, {"RS", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()}, {"RV", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()}, {"RVW", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedWater()}, - {"SOMAX", UnitSystem::measure::identity, simulator_.problem().vapparsActive(simulator_.episodeIndex())}, - {"PCSWM_OW", UnitSystem::measure::identity, enableHysteresis}, - {"KRNSW_OW", UnitSystem::measure::identity, enableHysteresis}, - {"PCSWM_GO", UnitSystem::measure::identity, enableHysteresis}, - {"KRNSW_GO", UnitSystem::measure::identity, enableHysteresis}, + {"SGMAX", UnitSystem::measure::identity, (enableNonWettingHysteresis && oilActive && gasActive)}, + {"SHMAX", UnitSystem::measure::identity, (enableWettingHysteresis && oilActive && gasActive)}, + {"SOMAX", UnitSystem::measure::identity, (enableNonWettingHysteresis && oilActive && waterActive) || simulator_.problem().vapparsActive(simulator_.episodeIndex())}, + {"SOMIN", UnitSystem::measure::identity, (enablePCHysteresis && oilActive && gasActive)}, + {"SWHY1", UnitSystem::measure::identity, (enablePCHysteresis && oilActive && waterActive)}, + {"SWMAX", UnitSystem::measure::identity, (enableWettingHysteresis && oilActive && waterActive)}, {"PPCW", UnitSystem::measure::pressure, enableSwatinit} }; diff --git a/opm/simulators/flow/FlowProblem.hpp b/opm/simulators/flow/FlowProblem.hpp index 2f98a9705..e0fc892bb 100644 --- a/opm/simulators/flow/FlowProblem.hpp +++ b/opm/simulators/flow/FlowProblem.hpp @@ -1397,6 +1397,12 @@ public: if (this->simulator().episodeIndex() == 0) { eclWriter_->writeInitialFIPReport(); } + + const bool invalidateFromHyst = updateHysteresis_(); + if (invalidateFromHyst) { + OPM_TIMEBLOCK(beginTimeStepInvalidateIntensiveQuantities); + this->model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0); + } } /*! diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp index 1dcfd599f..c18022d55 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp @@ -492,8 +492,12 @@ assignToSolution(data::Solution& sol) DataEntry{"RV", UnitSystem::measure::oil_gas_ratio, rv_}, DataEntry{"RVSAT", UnitSystem::measure::oil_gas_ratio, oilVaporizationFactor_}, DataEntry{"SALT", UnitSystem::measure::salinity, cSalt_}, - DataEntry{"SOMAX", UnitSystem::measure::identity, soMax_}, + DataEntry{"SGMAX", UnitSystem::measure::identity, sgmax_}, + DataEntry{"SHMAX", UnitSystem::measure::identity, shmax_}, + DataEntry{"SOMAX", UnitSystem::measure::identity, soMax_}, + DataEntry{"SOMIN", UnitSystem::measure::identity, somin_}, DataEntry{"SSOLVENT", UnitSystem::measure::identity, sSol_}, + DataEntry{"SWHY1", UnitSystem::measure::identity, swmin_}, DataEntry{"SWMAX", UnitSystem::measure::identity, swMax_}, DataEntry{"WATKR", UnitSystem::measure::identity, relativePermeability_[waterPhaseIdx]}, DataEntry{"WAT_DEN", UnitSystem::measure::density, density_[waterPhaseIdx]}, @@ -544,13 +548,9 @@ assignToSolution(data::Solution& sol) DataEntry{"DISPY", UnitSystem::measure::length, dispY_}, DataEntry{"DISPZ", UnitSystem::measure::length, dispZ_}, DataEntry{"DRSDTCON", UnitSystem::measure::gas_oil_ratio_rate, drsdtcon_}, - DataEntry{"KRNSW_GO", UnitSystem::measure::identity, krnSwMdcGo_}, - DataEntry{"KRNSW_OW", UnitSystem::measure::identity, krnSwMdcOw_}, DataEntry{"MECHPOTF", UnitSystem::measure::pressure, mechPotentialForce_}, DataEntry{"MICROBES", UnitSystem::measure::density, cMicrobes_}, DataEntry{"OXYGEN", UnitSystem::measure::density, cOxygen_}, - DataEntry{"PCSWM_GO", UnitSystem::measure::identity, pcSwMdcGo_}, - DataEntry{"PCSWM_OW", UnitSystem::measure::identity, pcSwMdcOw_}, DataEntry{"PERMFACT", UnitSystem::measure::identity, permFact_}, DataEntry{"PORV_RC", UnitSystem::measure::identity, rockCompPorvMultiplier_}, DataEntry{"PRESPOTF", UnitSystem::measure::pressure, mechPotentialPressForce_}, @@ -786,12 +786,8 @@ setRestart(const data::Solution& sol, std::pair{"BIOFILM", &cBiofilm_}, std::pair{"CALCITE", &cCalcite_}, std::pair{"FOAM", &cFoam_}, - std::pair{"KRNSW_GO", &krnSwMdcGo_}, - std::pair{"KRNSW_OW", &krnSwMdcOw_}, std::pair{"MICROBES", &cMicrobes_}, std::pair{"OXYGEN", &cOxygen_}, - std::pair{"PCSWM_GO", &pcSwMdcGo_}, - std::pair{"PCSWM_OW", &pcSwMdcOw_}, std::pair{"PERMFACT", &permFact_}, std::pair{"POLYMER", &cPolymer_}, std::pair{"PPCW", &ppcw_}, @@ -802,7 +798,12 @@ setRestart(const data::Solution& sol, std::pair{"RVW", &rvw_}, std::pair{"SALT", &cSalt_}, std::pair{"SALTP", &pSalt_}, + std::pair{"SGMAX", &sgmax_}, + std::pair{"SHMAX", &shmax_}, std::pair{"SOMAX", &soMax_}, + std::pair{"SOMIN", &somin_}, + std::pair{"SWHY1", &swmin_}, + std::pair{"SWMAX", &swMax_}, std::pair{"TEMP", &temperature_}, std::pair{"UREA", &cUrea_}, }; @@ -855,7 +856,9 @@ doAllocBuffers(const unsigned bufferSize, const bool log, const bool isRestart, const bool vapparsActive, - const bool enableHysteresis, + const bool enablePCHysteresis, + const bool enableNonWettingHysteresis, + const bool enableWettingHysteresis, const unsigned numTracers, const std::vector& enableSolTracers, const unsigned numOutputNnc) @@ -1133,11 +1136,41 @@ doAllocBuffers(const unsigned bufferSize, soMax_.resize(bufferSize, 0.0); } - if (enableHysteresis) { - pcSwMdcOw_.resize(bufferSize, 0.0); - krnSwMdcOw_.resize(bufferSize, 0.0); - pcSwMdcGo_.resize(bufferSize, 0.0); - krnSwMdcGo_.resize(bufferSize, 0.0); + if (enableNonWettingHysteresis) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)){ + if (FluidSystem::phaseIsActive(waterPhaseIdx)){ + soMax_.resize(bufferSize, 0.0); + } + if (FluidSystem::phaseIsActive(gasPhaseIdx)){ + sgmax_.resize(bufferSize, 0.0); + } + } else { + //TODO add support for gas-water + } + } + if (enableWettingHysteresis) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)){ + if (FluidSystem::phaseIsActive(waterPhaseIdx)){ + swMax_.resize(bufferSize, 0.0); + } + if (FluidSystem::phaseIsActive(gasPhaseIdx)){ + shmax_.resize(bufferSize, 0.0); + } + } else { + //TODO add support for gas-water + } + } + if (enablePCHysteresis) { + if (FluidSystem::phaseIsActive(oilPhaseIdx)){ + if (FluidSystem::phaseIsActive(waterPhaseIdx)){ + swmin_.resize(bufferSize, 0.0); + } + if (FluidSystem::phaseIsActive(gasPhaseIdx)){ + somin_.resize(bufferSize, 0.0); + } + } else { + //TODO add support for gas-water + } } if (eclState_.fieldProps().has_double("SWATINIT")) { diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.hpp b/opm/simulators/flow/GenericOutputBlackoilModule.hpp index 914c9a4a8..36652e9ad 100644 --- a/opm/simulators/flow/GenericOutputBlackoilModule.hpp +++ b/opm/simulators/flow/GenericOutputBlackoilModule.hpp @@ -338,7 +338,9 @@ protected: const bool log, const bool isRestart, const bool vapparsActive, - const bool enableHysteresis, + const bool enablePCHysteresis, + const bool enableNonWettingHysteresis, + const bool enableWettingHysteresis, unsigned numTracers, const std::vector& enableSolTracers, unsigned numOutputNnc); @@ -449,17 +451,17 @@ protected: ScalarBuffer mFracGas_; ScalarBuffer mFracCo2_; ScalarBuffer soMax_; - ScalarBuffer pcSwMdcOw_; - ScalarBuffer krnSwMdcOw_; - ScalarBuffer pcSwMdcGo_; - ScalarBuffer krnSwMdcGo_; + ScalarBuffer swMax_; + ScalarBuffer sgmax_; + ScalarBuffer shmax_; + ScalarBuffer somin_; + ScalarBuffer swmin_; ScalarBuffer ppcw_; ScalarBuffer gasDissolutionFactor_; ScalarBuffer oilVaporizationFactor_; ScalarBuffer bubblePointPressure_; ScalarBuffer dewPointPressure_; ScalarBuffer rockCompPorvMultiplier_; - ScalarBuffer swMax_; ScalarBuffer minimumOilPressure_; ScalarBuffer saturatedOilFormationVolumeFactor_; ScalarBuffer rockCompTransMultiplier_; diff --git a/opm/simulators/flow/OutputBlackoilModule.hpp b/opm/simulators/flow/OutputBlackoilModule.hpp index 8a87873e1..4d5c9557d 100644 --- a/opm/simulators/flow/OutputBlackoilModule.hpp +++ b/opm/simulators/flow/OutputBlackoilModule.hpp @@ -238,7 +238,9 @@ public: log, isRestart, problem.vapparsActive(std::max(simulator_.episodeIndex(), 0)), - problem.materialLawManager()->enableHysteresis(), + problem.materialLawManager()->enablePCHysteresis(), + problem.materialLawManager()->enableNonWettingHysteresis(), + problem.materialLawManager()->enableWettingHysteresis(), problem.tracerModel().numTracers(), problem.tracerModel().enableSolTracers(), problem.eclWriter()->getOutputNnc().size()); @@ -558,14 +560,6 @@ public: } } - if (!this->soMax_.empty()) - this->soMax_[globalDofIdx] - = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx)); - - if (!this->swMax_.empty()) - this->swMax_[globalDofIdx] - = std::max(getValue(fs.saturation(waterPhaseIdx)), problem.maxWaterSaturation(globalDofIdx)); - if (!this->minimumOilPressure_.empty()) this->minimumOilPressure_[globalDofIdx] = std::min(getValue(fs.pressure(oilPhaseIdx)), problem.minOilPressure(globalDofIdx)); @@ -583,16 +577,67 @@ public: const auto& matLawManager = problem.materialLawManager(); if (matLawManager->enableHysteresis()) { - if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) { - matLawManager->oilWaterHysteresisParams( - this->pcSwMdcOw_[globalDofIdx], this->krnSwMdcOw_[globalDofIdx], globalDofIdx); - } - if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) { - matLawManager->gasOilHysteresisParams( - this->pcSwMdcGo_[globalDofIdx], this->krnSwMdcGo_[globalDofIdx], globalDofIdx); - } - } + if (FluidSystem::phaseIsActive(oilPhaseIdx) + && FluidSystem::phaseIsActive(waterPhaseIdx)) { + Scalar somax; + Scalar swmax; + Scalar swmin; + matLawManager->oilWaterHysteresisParams( + somax, swmax, swmin, globalDofIdx); + + if (matLawManager->enableNonWettingHysteresis()) { + if (!this->soMax_.empty()) { + this->soMax_[globalDofIdx] = somax; + } + } + if (matLawManager->enableWettingHysteresis()) { + if (!this->swMax_.empty()) { + this->swMax_[globalDofIdx] = swmax; + } + } + if (matLawManager->enablePCHysteresis()) { + if (!this->swmin_.empty()) { + this->swmin_[globalDofIdx] = swmin; + } + } + } + + if (FluidSystem::phaseIsActive(oilPhaseIdx) + && FluidSystem::phaseIsActive(gasPhaseIdx)) { + Scalar sgmax; + Scalar shmax; + Scalar somin; + matLawManager->gasOilHysteresisParams( + sgmax, shmax, somin, globalDofIdx); + + if (matLawManager->enableNonWettingHysteresis()) { + if (!this->sgmax_.empty()) { + this->sgmax_[globalDofIdx] = sgmax; + } + } + if (matLawManager->enableWettingHysteresis()) { + if (!this->shmax_.empty()) { + this->shmax_[globalDofIdx] = shmax; + } + } + if (matLawManager->enablePCHysteresis()) { + if (!this->somin_.empty()) { + this->somin_[globalDofIdx] = somin; + } + } + } + } else { + + if (!this->soMax_.empty()) + this->soMax_[globalDofIdx] + = std::max(getValue(fs.saturation(oilPhaseIdx)), problem.maxOilSaturation(globalDofIdx)); + + if (!this->swMax_.empty()) + this->swMax_[globalDofIdx] + = std::max(getValue(fs.saturation(waterPhaseIdx)), problem.maxWaterSaturation(globalDofIdx)); + + } if (!this->ppcw_.empty()) { this->ppcw_[globalDofIdx] = matLawManager->oilWaterScaledEpsInfoDrainage(globalDofIdx).maxPcow; // printf("ppcw_[%d] = %lg\n", globalDofIdx, ppcw_[globalDofIdx]); @@ -1143,14 +1188,55 @@ public: if (simulator.problem().materialLawManager()->enableHysteresis()) { auto matLawManager = simulator.problem().materialLawManager(); - if (!this->pcSwMdcOw_.empty() && !this->krnSwMdcOw_.empty()) { + if (FluidSystem::phaseIsActive(oilPhaseIdx) + && FluidSystem::phaseIsActive(waterPhaseIdx)) { + Scalar somax = 2.0; + Scalar swmax = -2.0; + Scalar swmin = 2.0; + + if (matLawManager->enableNonWettingHysteresis()) { + if (!this->soMax_.empty()) { + somax = this->soMax_[elemIdx]; + } + } + if (matLawManager->enableWettingHysteresis()) { + if (!this->swMax_.empty()) { + swmax = this->swMax_[elemIdx]; + } + } + if (matLawManager->enablePCHysteresis()) { + if (!this->swmin_.empty()) { + swmin = this->swmin_[elemIdx]; + } + } matLawManager->setOilWaterHysteresisParams( - this->pcSwMdcOw_[elemIdx], this->krnSwMdcOw_[elemIdx], elemIdx); + somax, swmax, swmin, elemIdx); } - if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) { + if (FluidSystem::phaseIsActive(oilPhaseIdx) + && FluidSystem::phaseIsActive(gasPhaseIdx)) { + Scalar sgmax = 2.0; + Scalar shmax = -2.0; + Scalar somin = 2.0; + + if (matLawManager->enableNonWettingHysteresis()) { + if (!this->sgmax_.empty()) { + sgmax = this->sgmax_[elemIdx]; + } + } + if (matLawManager->enableWettingHysteresis()) { + if (!this->shmax_.empty()) { + shmax = this->shmax_[elemIdx]; + } + } + if (matLawManager->enablePCHysteresis()) { + if (!this->somin_.empty()) { + somin = this->somin_[elemIdx]; + } + } matLawManager->setGasOilHysteresisParams( - this->pcSwMdcGo_[elemIdx], this->krnSwMdcGo_[elemIdx], elemIdx); + sgmax, shmax, somin, elemIdx); } + } if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) {