thresholdPressures(): move the code which calculates the maximum gravity corrected pressure difference between EQLNUM regions to its own function

requested by [at]atgeirr

opm/core/utility/thresholdPressures.hpp

@@ -32,27 +32,19 @@
 
 namespace Opm
 {
-
+/// \brief Compute the maximum gravity corrected pressure difference of all
-    /// \brief Get a vector of pressure thresholds from EclipseState.
+///        equilibration regions given a reservoir state.
-    /// This function looks at EQLOPTS, THPRES and EQLNUM to determine
-    /// pressure thresholds.  It does not consider the case where the
-    /// threshold values are defaulted, in which case they should be
-    /// determined from the initial, equilibrated simulation state.
 /// \tparam    Grid           Type of grid object (UnstructuredGrid or CpGrid).
+/// \param[out] maxDp         The resulting pressure difference between equilibration regions
 /// \param[in] deck           Input deck, EQLOPTS and THPRES are accessed from it.
 /// \param[in] eclipseState   Processed eclipse state, EQLNUM is accessed from it.
 /// \param[in] grid           The grid to which the thresholds apply.
-    /// \return                   A vector of pressure thresholds, one
+/// \param[in] initialState   The state of the reservoir
-    ///                           for each face in the grid. A value
+/// \param[in] props          The object which calculates fluid properties
-    ///                           of zero means no threshold for that
+/// \param[in] gravity        The gravity constant
-    ///                           particular face. An empty vector is
-    ///                           returned if there is no THPRES
-    ///                           feature used in the deck.
-
-
-
 template <class Grid>
-    std::vector<double> thresholdPressures(const DeckConstPtr& deck,
+void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
+                  const DeckConstPtr& deck,
                   EclipseStateConstPtr eclipseState,
                   const Grid& grid,
                   const BlackoilState& initialState,
@@ -61,10 +53,7 @@ namespace Opm
 {
 
     const PhaseUsage& pu = props.phaseUsage();
-        SimulationConfigConstPtr simulationConfig = eclipseState->getSimulationConfig();
+
-        std::vector<double> thpres_vals;
-        if (simulationConfig->hasThresholdPressure()) {
-            std::shared_ptr<const ThresholdPressure> thresholdPressure = simulationConfig->getThresholdPressure();
     std::shared_ptr<GridProperty<int>> eqlnum = eclipseState->getIntGridProperty("EQLNUM");
     const auto& eqlnumData = eqlnum->getData();
 
@@ -255,9 +244,7 @@ namespace Opm
 
     // Calculate the maximum pressure potential difference between all PVT region
     // transitions of the initial solution.
-            std::map<std::pair<int, int>, double> maxDp;
     const int num_faces = UgGridHelpers::numFaces(grid);
-            thpres_vals.resize(num_faces, 0.0);
     const auto& fc = UgGridHelpers::faceCells(grid);
     for (int face = 0; face < num_faces; ++face) {
         const int c1 = fc(face, 0);
@@ -304,8 +291,45 @@ namespace Opm
                 maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));
         }
     }
+}
+
+    /// \brief Get a vector of pressure thresholds from EclipseState.
+    /// This function looks at EQLOPTS, THPRES and EQLNUM to determine
+    /// pressure thresholds.  It does not consider the case where the
+    /// threshold values are defaulted, in which case they should be
+    /// determined from the initial, equilibrated simulation state.
+    /// \tparam    Grid           Type of grid object (UnstructuredGrid or CpGrid).
+    /// \param[in] deck           Input deck, EQLOPTS and THPRES are accessed from it.
+    /// \param[in] eclipseState   Processed eclipse state, EQLNUM is accessed from it.
+    /// \param[in] maxDp          The maximum gravity corrected pressure differences between
+    ///                           the equilibration regions.
+    /// \param[in] grid           The grid to which the thresholds apply.
+    /// \return                   A vector of pressure thresholds, one
+    ///                           for each face in the grid. A value
+    ///                           of zero means no threshold for that
+    ///                           particular face. An empty vector is
+    ///                           returned if there is no THPRES
+    ///                           feature used in the deck.
+
+
+
+    template <class Grid>
+    std::vector<double> thresholdPressures(const DeckConstPtr& deck,
+                                           EclipseStateConstPtr eclipseState,
+                                           const Grid& grid,
+                                           const std::map<std::pair<int, int>, double>& maxDp)
+    {
+        SimulationConfigConstPtr simulationConfig = eclipseState->getSimulationConfig();
+        std::vector<double> thpres_vals;
+        if (simulationConfig->hasThresholdPressure()) {
+            std::shared_ptr<const ThresholdPressure> thresholdPressure = simulationConfig->getThresholdPressure();
+            std::shared_ptr<GridProperty<int>> eqlnum = eclipseState->getIntGridProperty("EQLNUM");
+            const auto& eqlnumData = eqlnum->getData();
 
             // Set threshold pressure values for each cell face.
+            const int num_faces = UgGridHelpers::numFaces(grid);
+            const auto& fc = UgGridHelpers::faceCells(grid);
+            const int* gc = UgGridHelpers::globalCell(grid);
             thpres_vals.resize(num_faces, 0.0);
             for (int face = 0; face < num_faces; ++face) {
                 const int c1 = fc(face, 0);
@@ -328,7 +352,7 @@ namespace Opm
                         // has been defaulted to the maximum pressure potential difference between
                         // these regions
                         const auto barrierId = std::make_pair(eq1, eq2);
-                        thpres_vals[face] = maxDp[barrierId];
+                        thpres_vals[face] = maxDp.at(barrierId);
                     }
                 }
             }