Merge pull request #335 from totto82/fix_thpress_init

WIP Fix thpress restart

ebos/eclproblem.hh

@@ -446,9 +446,6 @@ public:
             // create the ECL writer
             eclWriter_.reset(new EclWriterType(simulator));
 
-        // Loading the solution from a restart file is done recursively, we need this
-        // bool variable to signal the stop condition.
-        restartApplied = false;
     }
 
     /*!
@@ -519,19 +516,31 @@ public:
         readMaterialParameters_();
         readThermalParameters_();
         transmissibilities_.finishInit();
-        readInitialCondition_();
 
-        // Set the start time of the simulation
+        const auto& eclState = simulator.vanguard().eclState();
+        const auto& initconfig = eclState.getInitConfig();
         const auto& timeMap = simulator.vanguard().schedule().getTimeMap();
-        simulator.setStartTime( timeMap.getStartTime(/*timeStepIdx=*/0) );
+        if(initconfig.restartRequested()) {
+            // Set the start time of the simulation
+            simulator.setStartTime( timeMap.getStartTime(/*timeStepIdx=*/initconfig.getRestartStep()) );
+            simulator.setEpisodeIndex(initconfig.getRestartStep());
+            simulator.setEpisodeLength(0.0);
+            simulator.setTimeStepSize(0.0);
 
-        // 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 inside beginEpisode(). 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);
-        simulator.setTimeStepSize(0.0);
+            readEclRestartSolution_();
+        } else {
+            readInitialCondition_();
+            // Set the start time of the simulation
+            simulator.setStartTime( timeMap.getStartTime(/*timeStepIdx=*/0) );
+
+            // 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 inside beginEpisode(). 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);
+            simulator.setTimeStepSize(0.0);
+        }
 
         updatePffDofData_();
 
@@ -1228,23 +1237,11 @@ public:
             wellManager_.init(this->simulator().vanguard().eclState(), this->simulator().vanguard().schedule());
         }
 
-        // The initialSolutionApplied is called recursively by readEclRestartSolution_().
-        if (restartApplied)
-            return;
-
         // 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();
 
-        const auto& eclState = this->simulator().vanguard().eclState();
-        const auto& initconfig = eclState.getInitConfig();
-        if(initconfig.restartRequested()) {
-            restartApplied = true;
-            this->simulator().setEpisodeIndex(initconfig.getRestartStep());
-            readEclRestartSolution_();
-        }
-
         // release the memory of the EQUIL grid since it's no longer needed after this point
         this->simulator().vanguard().releaseEquilGrid();
 
@@ -1707,6 +1704,7 @@ private:
 
         for (size_t elemIdx = 0; elemIdx < numElems; ++elemIdx) {
             auto& elemFluidState = initialFluidStates_[elemIdx];
+            elemFluidState.setPvtRegionIndex(pvtRegionIndex(elemIdx));
             eclWriter_->eclOutputModule().initHysteresisParams(this->simulator(), elemIdx);
             eclWriter_->eclOutputModule().assignToFluidState(elemFluidState, elemIdx);
             if (enableSolvent)
@@ -1714,7 +1712,6 @@ private:
             if (enablePolymer)
                  polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);
         }
-        this->model().applyInitialSolution();
     }
 
     void readExplicitInitialCondition_()
@@ -2072,8 +2069,6 @@ private:
 
     PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
 
-    bool restartApplied;
-
 };
 } // namespace Ewoms
 

ebos/eclthresholdpressure.hh

@@ -116,15 +116,6 @@ public:
         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
@@ -132,10 +123,6 @@ public:
             throw std::runtime_error("The maximum number of supported equilibration regions is 255!");
         }
 
-        // allocate the array which specifies the threshold pressures
-        thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
-        thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
-
         // internalize the data specified using the EQLNUM keyword
         const std::vector<int>& equilRegionData =
             eclState.get3DProperties().getIntGridProperty("EQLNUM").getData();
@@ -147,6 +134,18 @@ public:
             elemEquilRegion_[elemIdx] = equilRegionData[cartElemIdx] - 1;
         }
 
+        /*
+          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;
+
+        // allocate the array which specifies the threshold pressures
+        thpres_.resize(numEquilRegions_*numEquilRegions_, 0.0);
+        thpresDefault_.resize(numEquilRegions_*numEquilRegions_, 0.0);
+
         computeDefaultThresholdPressures_();
         applyExplicitThresholdPressures_();
     }

ebos/eclwriter.hh

@@ -253,7 +253,9 @@ public:
             {"KRNSW_GO",  Opm::UnitSystem::measure::identity, enableHysteresis}
         };
 
-        std::vector<Opm::RestartKey> extraKeys = {{"OPMEXTRA", Opm::UnitSystem::measure::identity, false}};
+        const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
+        std::vector<Opm::RestartKey> extraKeys = {{"OPMEXTRA", Opm::UnitSystem::measure::identity, false},
+                                                  {"THPRES", Opm::UnitSystem::measure::pressure, inputThpres.active()}};
 
         unsigned episodeIdx = simulator_.episodeIndex();
         const auto& gridView = simulator_.vanguard().gridView();
@@ -265,7 +267,6 @@ public:
             unsigned globalIdx = collectToIORank_.localIdxToGlobalIdx(elemIdx);
             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();