From dcb1e96606a6446d719641ffe13b172d1596596a Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Fri, 7 Sep 2018 10:19:00 +0200
Subject: [PATCH] Output and read PPCW when SWATINIT is enabled

---
 ebos/ecloutputblackoilmodule.hh | 24 ++++++++++++++++++++++--
 ebos/eclwriter.hh               |  4 +++-
 2 files changed, 25 insertions(+), 3 deletions(-)

diff --git a/ebos/ecloutputblackoilmodule.hh b/ebos/ecloutputblackoilmodule.hh
index 0e9a93407..44f9864da 100644
--- a/ebos/ecloutputblackoilmodule.hh
+++ b/ebos/ecloutputblackoilmodule.hh
@@ -251,6 +251,9 @@ public:
             krnSwMdcGo_.resize(bufferSize, 0.0);
         }
 
+        if (simulator_.vanguard().eclState().get3DProperties().hasDeckDoubleGridProperty("SWATINIT"))
+            ppcw_.resize(bufferSize, 0.0);
+
         if (FluidSystem::enableDissolvedGas() && rstKeywords["RSSAT"] > 0) {
             rstKeywords["RSSAT"] = 0;
             gasDissolutionFactor_.resize(bufferSize, 0.0);
@@ -495,6 +498,10 @@ public:
                 }
             }
 
+
+            if (ppcw_.size() > 0)
+                ppcw_[globalDofIdx] = matLawManager->oilWaterScaledEpsInfoDrainage(globalDofIdx).maxPcow;
+
             // hack to make the intial output of rs and rv Ecl compatible.
             // For cells with swat == 1 Ecl outputs; rs = rsSat and rv=rvSat, in all but the initial step
             // where it outputs rs and rv values calculated by the initialization. To be compatible we overwrite
@@ -713,8 +720,8 @@ public:
         if (FluidSystem::phaseIsActive(gasPhaseIdx) && saturation_[gasPhaseIdx].size() > 0) {
             sol.insert("SGAS", Opm::UnitSystem::measure::identity, std::move(saturation_[gasPhaseIdx]), Opm::data::TargetType::RESTART_SOLUTION);
         }
-        if (pcSwMdcOw_.size() > 0) {
-            sol.insert ("PPCW", Opm::UnitSystem::measure::identity, std::move(pcSwMdcOw_), Opm::data::TargetType::RESTART_SOLUTION);
+        if (ppcw_.size() > 0) {
+            sol.insert ("PPCW", Opm::UnitSystem::measure::pressure, std::move(ppcw_), Opm::data::TargetType::RESTART_SOLUTION);
         }
 
         if (gasDissolutionFactor_.size() > 0) {
@@ -771,6 +778,9 @@ public:
             sol.insert("GASKR", Opm::UnitSystem::measure::identity, std::move(relativePermeability_[gasPhaseIdx]), Opm::data::TargetType::RESTART_AUXILIARY);
         }
 
+        if (pcSwMdcOw_.size() > 0)
+            sol.insert ("PCSWM_OW", Opm::UnitSystem::measure::identity, std::move(pcSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY);
+
         if (krnSwMdcOw_.size() > 0)
             sol.insert ("KRNSW_OW", Opm::UnitSystem::measure::identity, std::move(krnSwMdcOw_), Opm::data::TargetType::RESTART_AUXILIARY);
 
@@ -945,6 +955,9 @@ public:
             pcSwMdcGo_[elemIdx] = sol.data("PCSWM_GO")[globalDofIndex];
         if (krnSwMdcGo_.size() > 0 && sol.has("KRNSW_GO"))
             krnSwMdcGo_[elemIdx] = sol.data("KRNSW_GO")[globalDofIndex];
+        if (ppcw_.size() > 0 && sol.has("PPCW"))
+            ppcw_[elemIdx] = sol.data("PPCW")[globalDofIndex];
+
     }
 
     template <class FluidState>
@@ -1004,6 +1017,12 @@ public:
             }
         }
 
+        if (simulator_.vanguard().eclState().get3DProperties().hasDeckDoubleGridProperty("SWATINIT")) {
+            auto oilWaterScaledEpsInfoDrainage = simulator.problem().materialLawManager()->oilWaterScaledEpsInfoDrainagePointerReferenceHack(elemIdx);
+            oilWaterScaledEpsInfoDrainage->maxPcow =  ppcw_[elemIdx];
+        }
+
+
     }
 
     Scalar getSolventSaturation(unsigned elemIdx) const {
@@ -1317,6 +1336,7 @@ private:
     ScalarBuffer krnSwMdcOw_;
     ScalarBuffer pcSwMdcGo_;
     ScalarBuffer krnSwMdcGo_;
+    ScalarBuffer ppcw_;
     ScalarBuffer bubblePointPressure_;
     ScalarBuffer dewPointPressure_;
     std::vector<int> failedCellsPb_;
diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh
index 512283ea3..5b04e81ba 100644
--- a/ebos/eclwriter.hh
+++ b/ebos/eclwriter.hh
@@ -319,6 +319,7 @@ public:
     void restartBegin()
     {
         bool enableHysteresis = simulator_.problem().materialLawManager()->enableHysteresis();
+        bool enableSwatinit = simulator_.vanguard().eclState().get3DProperties().hasDeckDoubleGridProperty("SWATINIT");
         std::vector<Opm::RestartKey> solutionKeys{
             {"PRESSURE" , Opm::UnitSystem::measure::pressure},
             {"SWAT" ,     Opm::UnitSystem::measure::identity, static_cast<bool>(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx))},
@@ -330,7 +331,8 @@ public:
             {"PCSWM_OW",  Opm::UnitSystem::measure::identity, enableHysteresis},
             {"KRNSW_OW",  Opm::UnitSystem::measure::identity, enableHysteresis},
             {"PCSWM_GO",  Opm::UnitSystem::measure::identity, enableHysteresis},
-            {"KRNSW_GO",  Opm::UnitSystem::measure::identity, enableHysteresis}
+            {"KRNSW_GO",  Opm::UnitSystem::measure::identity, enableHysteresis},
+            {"PPCW",      Opm::UnitSystem::measure::pressure, enableSwatinit}
         };
 
         const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();