From 3b5e746702518e4def2ef66544d3448ab60c4458 Mon Sep 17 00:00:00 2001
From: Joakim Hove <joakim.hove@gmail.com>
Date: Thu, 31 May 2018 11:58:27 +0200
Subject: [PATCH] Load the THPRES values from the restart file

---
 ebos/eclproblem.hh           | 14 ++++++++------
 ebos/eclthresholdpressure.hh | 15 ++++++++++++++-
 ebos/eclwriter.hh            | 13 ++++++++++---
 3 files changed, 32 insertions(+), 10 deletions(-)

diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index ff9dcd46f..105f3f178 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -399,7 +399,8 @@ public:
             // create the ECL writer
             eclWriter_.reset(new EclWriterType(simulator));
 
-        // Hack to compute the initial thpressure values for restarts
+        // Loading the solution from a restart file is done recursively, we need this
+        // bool variable to signal the stop condition.
         restartApplied = false;
     }
 
@@ -847,9 +848,11 @@ public:
     { return thresholdPressures_.thresholdPressure(elem1Idx, elem2Idx); }
 
 
-    const EclThresholdPressure<TypeTag>& thresholdPressure() const {
-        return thresholdPressures_;
-    }
+    const EclThresholdPressure<TypeTag>& thresholdPressure() const
+    { return thresholdPressures_; }
+
+    EclThresholdPressure<TypeTag>& thresholdPressure()
+    { return thresholdPressures_; }
 
     /*!
      * \copydoc FvBaseMultiPhaseProblem::porosity
@@ -1178,8 +1181,7 @@ public:
             wellManager_.init(this->simulator().vanguard().eclState(), this->simulator().vanguard().schedule());
         }
 
-        // the initialSolutionApplied is called recursively by readEclRestartSolution_()
-        // in order to setup the inital threshold pressures correctly
+        // The initialSolutionApplied is called recursively by readEclRestartSolution_().
         if (restartApplied)
             return;
 
diff --git a/ebos/eclthresholdpressure.hh b/ebos/eclthresholdpressure.hh
index 83fcf1026..44b422526 100644
--- a/ebos/eclthresholdpressure.hh
+++ b/ebos/eclthresholdpressure.hh
@@ -86,6 +86,10 @@ public:
         enableThresholdPressure_ = false;
     }
 
+
+    void setFromRestart(const std::vector<Scalar>& values)
+    { thpres_ = values; }
+
     /*!
      * \brief Actually compute the threshold pressures over a face as a pre-compute step.
      */
@@ -105,10 +109,19 @@ public:
         const auto& eclState = vanguard.eclState();
         const auto& simConfig = eclState.getSimulationConfig();
 
-        enableThresholdPressure_ = simConfig.hasThresholdPressure();
+        enableThresholdPressure_ = simConfig.useThresholdPressure();
         if (!enableThresholdPressure_)
             return;
 
+        /*
+          If this is a restart run the ThresholdPressure object will be active,
+          but it will *not* be properly initialized with numerical values. The
+          values must instead come from the THPRES vector in the restart file.
+        */
+        if (simConfig.getThresholdPressure().restart())
+            return;
+
+
         numEquilRegions_ = eclState.getTableManager().getEqldims().getNumEquilRegions();
         if (numEquilRegions_ > 0xff) {
             // make sure that the index of an equilibration region can be stored in a
diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh
index 8c8e43031..5433d2d0a 100644
--- a/ebos/eclwriter.hh
+++ b/ebos/eclwriter.hh
@@ -210,7 +210,7 @@ public:
             if (!isSubStep)
                 restartValue.addExtra("OPMEXTRA", std::vector<double>(1, nextStepSize));
 
-            if (simConfig.hasThresholdPressure())
+            if (simConfig.useThresholdPressure())
                 restartValue.addExtra("THPRES", Opm::UnitSystem::measure::pressure, simulator_.problem().thresholdPressure().data());
 
             // first, create a tasklet to write the data for the current time step to disk
@@ -255,10 +255,17 @@ public:
         unsigned numElements = gridView.size(/*codim=*/0);
         eclOutputModule_.allocBuffers(numElements, episodeIdx, /*isSubStep=*/false, /*log=*/false);
 
-        auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
+        auto restartValues = eclIO_->loadRestart(solution_keys, extra_keys);
         for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
             unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
-            eclOutputModule_.setRestart(restart_values.solution, elemIdx, globalIdx);
+            eclOutputModule_.setRestart(restartValues.solution, elemIdx, globalIdx);
+        }
+        const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
+        if (inputThpres.active()) {
+            Simulator& mutableSimulator = const_cast<Simulator&>(simulator_);
+            auto& thpres = mutableSimulator.problem().thresholdPressure();
+            const auto& thpresValues = restartValues.getExtra("THPRES");
+            thpres.setFromRestart(thpresValues);
         }
     }