Merge pull request #5402 from totto82/output_wetting_hyst

Output maximum/minimum saturations directly for restart hysteresis

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<TypeTag, Properties::EnableOpmRstFile>();
         bool read_temp = enableEnergy || (opm_rst_file && enableTemperature);
         std::vector<RestartKey> solutionKeys{
             {"PRESSURE", UnitSystem::measure::pressure},
-            {"SWAT", UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))},
+            {"SWAT", UnitSystem::measure::identity, static_cast<bool>(waterActive)},
-            {"SGAS", UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx))},
+            {"SGAS", UnitSystem::measure::identity, static_cast<bool>(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())},
+            {"SGMAX", UnitSystem::measure::identity, (enableNonWettingHysteresis && oilActive && gasActive)},
-            {"PCSWM_OW", UnitSystem::measure::identity, enableHysteresis},
+            {"SHMAX", UnitSystem::measure::identity, (enableWettingHysteresis && oilActive && gasActive)},
-            {"KRNSW_OW", UnitSystem::measure::identity, enableHysteresis},
+            {"SOMAX", UnitSystem::measure::identity, (enableNonWettingHysteresis && oilActive && waterActive) || simulator_.problem().vapparsActive(simulator_.episodeIndex())},
-            {"PCSWM_GO", UnitSystem::measure::identity, enableHysteresis},
+            {"SOMIN", UnitSystem::measure::identity, (enablePCHysteresis && oilActive && gasActive)},
-            {"KRNSW_GO", UnitSystem::measure::identity, enableHysteresis},
+            {"SWHY1", UnitSystem::measure::identity, (enablePCHysteresis && oilActive && waterActive)},
+            {"SWMAX", UnitSystem::measure::identity, (enableWettingHysteresis && oilActive && waterActive)},
             {"PPCW", UnitSystem::measure::pressure, enableSwatinit}
         };
 

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);
+        }
     }
 
     /*!

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{"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<bool>& enableSolTracers,
                const unsigned numOutputNnc)
@@ -1133,11 +1136,41 @@ doAllocBuffers(const unsigned bufferSize,
         soMax_.resize(bufferSize, 0.0);
     }
 
-    if (enableHysteresis) {
+    if (enableNonWettingHysteresis) {
-        pcSwMdcOw_.resize(bufferSize, 0.0);
+        if (FluidSystem::phaseIsActive(oilPhaseIdx)){
-        krnSwMdcOw_.resize(bufferSize, 0.0);
+            if (FluidSystem::phaseIsActive(waterPhaseIdx)){
-        pcSwMdcGo_.resize(bufferSize, 0.0);
+                soMax_.resize(bufferSize, 0.0);
-        krnSwMdcGo_.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")) {

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<bool>& enableSolTracers,
                         unsigned numOutputNnc);
@@ -449,17 +451,17 @@ protected:
     ScalarBuffer mFracGas_;
     ScalarBuffer mFracCo2_;
     ScalarBuffer soMax_;
-    ScalarBuffer pcSwMdcOw_;
+    ScalarBuffer swMax_;
-    ScalarBuffer krnSwMdcOw_;
+    ScalarBuffer sgmax_;
-    ScalarBuffer pcSwMdcGo_;
+    ScalarBuffer shmax_;
-    ScalarBuffer krnSwMdcGo_;
+    ScalarBuffer somin_;
+    ScalarBuffer swmin_;
     ScalarBuffer ppcw_;
     ScalarBuffer gasDissolutionFactor_;
     ScalarBuffer oilVaporizationFactor_;
     ScalarBuffer bubblePointPressure_;
     ScalarBuffer dewPointPressure_;
     ScalarBuffer rockCompPorvMultiplier_;
-    ScalarBuffer swMax_;
     ScalarBuffer minimumOilPressure_;
     ScalarBuffer saturatedOilFormationVolumeFactor_;
     ScalarBuffer rockCompTransMultiplier_;

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()) {
+                if (FluidSystem::phaseIsActive(oilPhaseIdx) 
+                    && FluidSystem::phaseIsActive(waterPhaseIdx)) {
+                        Scalar somax;
+                        Scalar swmax;
+                        Scalar swmin;
+
                         matLawManager->oilWaterHysteresisParams(
-                        this->pcSwMdcOw_[globalDofIdx], this->krnSwMdcOw_[globalDofIdx], globalDofIdx);
+                            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 (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) {
-                    matLawManager->gasOilHysteresisParams(
-                        this->pcSwMdcGo_[globalDofIdx], this->krnSwMdcGo_[globalDofIdx], globalDofIdx);
                     }
                 }
 
+                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 (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];
+                    }
                 }
-            if (!this->pcSwMdcGo_.empty() && !this->krnSwMdcGo_.empty()) {
                 matLawManager->setGasOilHysteresisParams(
-                    this->pcSwMdcGo_[elemIdx], this->krnSwMdcGo_[elemIdx], elemIdx);
+                        sgmax, shmax, somin, elemIdx);
             }
+
         }
 
         if (simulator_.vanguard().eclState().fieldProps().has_double("SWATINIT")) {