Merge pull request #3682 from bska/refactor-flux-calc

Refactor Flux Calculator to Eliminate Warnings and Reduce Clutter

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.cpp

@@ -18,16 +18,20 @@
 */
 
 #include <opm/utility/ECLFluxCalc.hpp>
+
 #include <opm/utility/ECLPvtCommon.hpp>
-#include <opm/utility/ECLUnitHandling.hpp>
 #include <opm/utility/ECLSaturationFunc.hpp>
+#include <opm/utility/ECLUnitHandling.hpp>
+
 #include <opm/utility/imported/Units.hpp>
 
 #include <algorithm>
+#include <cstddef>
 #include <exception>
 #include <functional>
 #include <iterator>
 #include <stdexcept>
+#include <string>
 #include <utility>
 
 #include <ert/ecl/ecl_kw_magic.h>
@@ -174,7 +178,6 @@ namespace {
             };
         }
     }
-
 } // Anonymous
 
 namespace Opm
@@ -342,41 +345,18 @@ namespace Opm
         // Allocate space for storing the cell values.
         dyn_data.density.assign(this->graph_.numCells(), 0.0);
 
+        // Extract phase saturation.
+        dyn_data.saturation =
+            ::Opm::phaseSaturation(this->graph_, rstrt, phase);
+
         switch (phase) {
         case ECLPhaseIndex::Aqua:
-            dyn_data.saturation = this->graph_.rawLinearisedCellData<double>(rstrt, "SWAT");
             return this->watPVT(std::move(dyn_data));
 
         case ECLPhaseIndex::Liquid:
-            dyn_data.saturation = this->graph_.rawLinearisedCellData<double>(rstrt, "SOIL");
+            return this->oilPVT(rstrt, std::move(dyn_data));
-            if (!dyn_data.saturation.empty()) {
-                return this->oilPVT(rstrt, std::move(dyn_data));
-            } else {
-                // SOIL vector not provided. Compute from SWAT and/or SGAS.
-                // may read two times
-                auto sw = this->graph_.rawLinearisedCellData<double>(rstrt, "SWAT");
-                auto sg = this->graph_.rawLinearisedCellData<double>(rstrt, "SGAS");
-                std::vector<double>& so = dyn_data.saturation;
-                so.assign(this->graph_.numCells(), 1.0);
-                auto adjust_So_for_other_phase =
-                    [&so](const std::vector<double>& s)
-                {
-                    std::transform(std::begin(so), std::end(so),
-                                   std::begin(s) ,
-                                   std::begin(so), std::minus<double>());
-                };
-                if (sg.size() == this->graph_.numCells()) {
-                    adjust_So_for_other_phase(sg);
-                }
-
-                if (sw.size() == this->graph_.numCells()) {
-                    adjust_So_for_other_phase(sw);
-                }
-                return this->oilPVT(rstrt, std::move(dyn_data));
-            }
 
         case ECLPhaseIndex::Vapour:
-            dyn_data.saturation = this->graph_.rawLinearisedCellData<double>(rstrt, "SGAS");
             return this->gasPVT(rstrt, std::move(dyn_data));
         }
 
@@ -385,7 +365,8 @@ namespace Opm
         };
     }
 
-    double ECLFluxCalc::surfaceDensity(const ECLPhaseIndex   phase) const{
+    double ECLFluxCalc::surfaceDensity(const ECLPhaseIndex phase) const
+    {
         switch (phase) {
         case ECLPhaseIndex::Aqua:
             return this->pvtWat_->surfaceMassDensity(0);
@@ -396,6 +377,11 @@ namespace Opm
         case ECLPhaseIndex::Vapour:
             return this->pvtGas_->surfaceMassDensity(0);
         }
+
+        throw std::invalid_argument {
+            "Unsupported Phase Index " +
+            std::to_string(static_cast<std::size_t>(phase))
+        };
     }
 
 

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLFluxCalc.hpp

@@ -63,13 +63,13 @@ namespace Opm
         /// graph.neighbours() \endcode.
         ///
         /// \param[in] rstrt ECL Restart data set from which to extract
-        ///            relevant data per cell.
+        ///    relevant data per cell.
         ///
         /// \param[in] phase Canonical phase for which to retrive flux.
         ///
-        /// \return Flux values corresponding to selected phase.
+        /// \return Flux values corresponding to selected phase.  Empty if
-        ///         Empty if required data is missing.
+        ///    requisite data is missing.  Numerical values in SI units
-        ///         Numerical values in SI units (rm^3/s).
+        ///    (rm^3/s).
         std::vector<double>
         flux(const ECLRestartData& rstrt,
              const ECLPhaseIndex   phase) const;
@@ -78,20 +78,19 @@ namespace Opm
         /// graph.neighbours() \endcode.
         ///
         /// \param[in] rstrt ECL Restart data set from which to extract
-        ///            relevant data per cell.
+        ///    relevant data per cell.
         ///
         /// \param[in] phase Canonical phase for which to retrive flux.
         ///
-        /// \return Mass flux values corresponding to selected phase.
+        /// \return Mass flux values corresponding to selected phase.  Empty
-        ///         Empty if required data is missing.
+        ///    if requisite data is missing.  Numerical values in SI units
-        ///         Numerical values in SI units (kg/s).
+        ///    (kg/s).
         std::vector<double>
         massflux(const ECLRestartData& rstrt,
                  const ECLPhaseIndex   phase) const;
 
-        /// Return type for the phaseProperties() method,
+        /// Return type for the phaseProperties() method, encapsulates
-        /// encapsulates dynamic properties for a single
+        /// dynamic properties for a single phase.
-        /// phase.
         struct DynamicData
         {
             std::vector<double> pressure;
@@ -103,27 +102,25 @@ namespace Opm
         /// Retrive dynamical properties of a single phase on all cells.
         ///
         /// \param[in] rstrt ECL Restart data set from which to extract
-        ///            relevant data per cell.
+        ///    relevant data per cell.
         ///
         /// \param[in] phase Canonical phase for which to retrive properties.
         ///
-        /// \return DynamicData struct containing cell-values for phase properties.
+        /// \return DynamicData struct containing cell-values for phase
-        ///         Numerical values in SI units (kg/s).
+        ///    properties.  Numerical values in SI units (kg/s).
         DynamicData phaseProperties(const ECLRestartData& rstrt,
                                     const ECLPhaseIndex   phase) const;
 
         /// Retrive the constant surface density of a phase.
         ///
-        /// \param[in] phase Canonical phase for which to retrive the surface density.
+        /// \param[in] phase Canonical phase for which to retrive the
+        ///    surface density.
         ///
         /// \return Density of given phase at surface conditions.
-        ///         Numerical value in SI units (kg/m^3).
+        ///    Numerical value in SI units (kg/m^3).
         double surfaceDensity(const ECLPhaseIndex phase) const;
 
-
     private:
-
-
         double singleFlux(const int connection,
                           const DynamicData& dyn_data) const;
 

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp

@@ -121,6 +121,15 @@ namespace {
         return std::make_pair(ToSI::disGas(*u),
                               deadOilUnitConverter(*u));
     }
+
+    std::vector<bool>::size_type const_compr_index()
+    {
+#if defined(LOGIHEAD_CONSTANT_OILCOMPR_INDEX)
+        return LOGIHEAD_CONSTANT_OILCOMPR_INDEX;
+#else
+        return (39 - 1);        // Reverse engineering...
+#endif  // LOGIHEAD_CONSTANT_OILCOMPR_INDEX
+    }
 }
 
 // ---------------------------------------------------------------------
@@ -202,12 +211,12 @@ public:
     }
 
 private:
-    double po_ref_       { 1.0 };
+    double po_ref_ { 1.0 };
-    double fvf_     { 1.0 }; // B
+    double fvf_    { 1.0 }; // B
-    double visc_ { 1.0 }; // mu
+    double visc_   { 1.0 }; // mu
-    double Co_           { 1.0 };
+    double Co_     { 1.0 };
     double cv_     { 0.0 }; // Cv
-    double rhoS_         { 0.0 };
+    double rhoS_   { 0.0 };
 
     double recipFvf(const double po) const
     {
@@ -245,20 +254,20 @@ private:
 };
 
 DeadOilConstCompr::DeadOilConstCompr(ElemIt               xBegin,
-                     ElemIt               xEnd,
+                                     ElemIt               xEnd,
-                     const ConvertUnits&  convert,
+                                     const ConvertUnits&  convert,
-                     std::vector<ElemIt>& colIt)
+                                     std::vector<ElemIt>& colIt)
 {
     // Recall: Table is
     //
     //    [ Po, Bo, Co, mu_o, Cv ]
     //
-    // xBegin is Pw, colIt is remaining four columns.
+    // xBegin is Po, colIt is remaining four columns.
 
-    this->fvf_     = convert.column[0](*colIt[0]); // Bo
+    this->fvf_  = convert.column[0](*colIt[0]); // Bo
-    this->Co_           = convert.column[1](*colIt[1]); // Co
+    this->Co_   = convert.column[1](*colIt[1]); // Co
     this->visc_ = convert.column[2](*colIt[2]); // mu_o
-    this->cv_     = convert.column[3](*colIt[3]); // Cw - Cv
+    this->cv_   = convert.column[3](*colIt[3]); // Cv
 
     // Honour requirement that constructor advances column iterators.
     const auto N = std::distance(xBegin, xEnd);
@@ -280,7 +289,6 @@ DeadOilConstCompr::DeadOilConstCompr(ElemIt               xBegin,
     }
 }
 
-
 // =====================================================================
 
 class DeadOil : public PVxOBase
@@ -407,7 +415,7 @@ private:
 namespace {
     std::vector<std::unique_ptr<PVxOBase>>
     createDeadOil(const ::Opm::ECLPropTableRawData& raw,
-                  const bool   const_compr,
+                  const bool                        const_compr,
                   const int                         usys)
     {
         using PVTInterp = std::unique_ptr<PVxOBase>;
@@ -504,7 +512,7 @@ namespace {
 
     std::vector<std::unique_ptr<PVxOBase>>
     createPVTFunction(const ::Opm::ECLPropTableRawData& raw,
-                      const bool   const_compr,
+                      const bool                        const_compr,
                       const int                         usys)
     {
         if (raw.numPrimary == 0) {
@@ -613,7 +621,7 @@ private:
 Opm::ECLPVT::Oil::Impl::
 Impl(const ECLPropTableRawData& raw,
      const int                  usys,
-     const bool     const_compr,
+     const bool                 const_compr,
      std::vector<double>        rhoS)
     : eval_(createPVTFunction(raw, const_compr, usys))
     , rhoS_(std::move(rhoS))
@@ -778,9 +786,8 @@ fromECLOutput(const ECLInitFileData& init)
         return OPtr{};
     }
 
-    const auto& lh = init.keywordData<bool>(LOGIHEAD_KW);
+    const auto& lh             = init.keywordData<bool>(LOGIHEAD_KW);
-    const int LOGIHEAD_CONST_COMPR_INDEX = 38;
+    const auto  is_const_compr = static_cast<bool>(lh[const_compr_index()]);
-    const bool is_const_compr = lh[LOGIHEAD_CONST_COMPR_INDEX];
 
     auto raw = ::Opm::ECLPropTableRawData{};
 

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp

@@ -48,6 +48,13 @@
 #include <ert/ecl/ecl_kw_magic.h>
 
 namespace {
+    std::vector<double>
+    gas_saturation(const ::Opm::ECLGraph&       G,
+                   const ::Opm::ECLRestartData& rstrt)
+    {
+        return G.rawLinearisedCellData<double>(rstrt, "SGAS");
+    }
+
     std::vector<double>
     oil_saturation(const std::vector<double>&   sg,
                    const std::vector<double>&   sw,
@@ -84,6 +91,13 @@ namespace {
         return so;
     }
 
+    std::vector<double>
+    water_saturation(const ::Opm::ECLGraph&       G,
+                     const ::Opm::ECLRestartData& rstrt)
+    {
+        return G.rawLinearisedCellData<double>(rstrt, "SWAT");
+    }
+
     std::vector<int>
     satnumVector(const ::Opm::ECLGraph&        G,
                  const ::Opm::ECLInitFileData& init)
@@ -1881,8 +1895,8 @@ kro(const ECLGraph&       G,
         return kr;
     }
 
-    const auto& sg = G.rawLinearisedCellData<double>(rstrt, "SGAS");
+    const auto& sg = gas_saturation(G, rstrt);
-    const auto& sw = G.rawLinearisedCellData<double>(rstrt, "SWAT");
+    const auto& sw = water_saturation(G, rstrt);
 
     auto so_g = oil_saturation(sg, sw, G, rstrt);
     auto so_w = so_g;
@@ -2010,7 +2024,7 @@ krg(const ECLGraph&       G,
         return kr;
     }
 
-    auto sg = G.rawLinearisedCellData<double>(rstrt, "SGAS");
+    auto sg = gas_saturation(G, rstrt);
 
     if (enableHorizontalEPS(scaling) && this->eps_) {
         this->eps_->scaleKrGas(this->rmap_, sg);
@@ -2162,7 +2176,7 @@ krw(const ECLGraph&       G,
         return kr;
     }
 
-    auto sw = G.rawLinearisedCellData<double>(rstrt, "SWAT");
+    auto sw = water_saturation(G, rstrt);
 
     if (enableHorizontalEPS(scaling) && this->eps_) {
         this->eps_->scaleKrWat(this->rmap_, sw);
@@ -2564,3 +2578,31 @@ getSatFuncCurve(const std::vector<RawCurve>& func,
 {
     return this->pImpl_->getSatFuncCurve(func, activeCell, scaling);
 }
+
+// =====================================================================
+
+std::vector<double>
+Opm::phaseSaturation(const ECLGraph&       G,
+                     const ECLRestartData& rstrt,
+                     const ECLPhaseIndex   phase)
+{
+    switch (phase) {
+    case ECLPhaseIndex::Aqua:
+        return water_saturation(G, rstrt);
+
+    case ECLPhaseIndex::Liquid: {
+        const auto sg = gas_saturation(G, rstrt);
+        const auto sw = water_saturation(G, rstrt);
+
+        return oil_saturation(sg, sw, G, rstrt);
+    }
+
+    case ECLPhaseIndex::Vapour:
+        return gas_saturation(G, rstrt);
+    }
+
+    throw std::invalid_argument {
+        "Unsupported Phase Index " +
+        std::to_string(static_cast<std::size_t>(phase))
+    };
+}

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp

@@ -41,6 +41,34 @@ namespace Opm {
     class ECLRestartData;
     class ECLInitFileData;
 
+    /// Extract phase saturation of single phase for all active cells in all
+    /// grids.
+    ///
+    /// Handles the case of oil saturation being explicitly stored in a
+    /// result set or implicitly defined from the gas and/or water
+    /// saturations.
+    ///
+    /// \param[in] G Connected topology of current model's active cells.
+    ///    Needed to linearise phase saturations (e.g., SOIL) that are
+    ///    distributed on local grids to all of the model's active cells
+    ///    (\code member function G.rawLinearisedCellData() \endcode).
+    ///
+    /// \param[in] rstrt ECLIPSE restart vectors.  Result set view
+    ///    assumed to be positioned at a particular report step of
+    ///    interest.
+    ///
+    /// \param[in] phase Phase for which to extract the phase saturation
+    ///    values.
+    ///
+    /// \return Phase saturation values of active phase \p phase for all
+    ///    active cells in model \p G.  Empty if phase \p phase is not
+    ///    actually active in the current result set or if the saturation
+    ///    values are not stored on the current report/restart step.
+    std::vector<double>
+    phaseSaturation(const ECLGraph&       G,
+                    const ECLRestartData& rstrt,
+                    const ECLPhaseIndex   phase);
+
     /// Gateway to engine for computing relative permeability values based
     /// on tabulated saturation functions in ECL output.
     class ECLSaturationFunc