Pull Dake Model Correlation Parameters Out to Helper Structure

Some of these parameters must be output to the restart file so make
it easier to refer to them as a group.

opm/input/eclipse/Schedule/Well/WDFAC.hpp

@@ -21,6 +21,7 @@
 #define WDFAC_HPP_HEADER_INCLUDED
 
 namespace Opm {
+    class Connection;
     class DeckRecord;
     class WellConnections;
 } // namespace Opm
@@ -31,49 +32,177 @@ namespace Opm { namespace RestartIO {
 
 namespace Opm {
 
-    enum class WDFACTYPE {
-        NONE = 1,
-        DFACTOR = 2,
-        DAKEMODEL = 3,
-        CON_DFACTOR = 4
-    };
-
     class WDFAC
     {
     public:
+        /// Parameters for Dake's D-factor correlation model.
+        ///
+        /// In particular, holds the coefficient 'A' and the exponents 'B'
+        /// and 'C' of the correlation relation
+        ///
+        ///   D = A * (Ke/K0)**B * porosity**C * Ke / (h * rw) * (sg_g/mu_g)
+        ///
+        /// in which
+        ///
+        ///   * Ke is the connection's effective permeability (sqrt(Kx*Ky)
+        ///     in the case of a vertical connection)
+        ///
+        ///   * K0 is a reference/background permeability scale (1mD)
+        ///
+        ///   * h is the effective length of the connection's perforation
+        ///     interval (dz*ntg in the case of a vertical connection)
+        ///
+        ///   * rw is the connection's wellbore radius
+        ///
+        ///   * sg_g is the specific gravity of surface condition gas
+        ///     relative to surface condition air
+        ///
+        ///   * mu_g is the reservoir condition viscosity of the free gas phase.
+        struct Correlation
+        {
+            /// Multiplicative coefficient 'A'.
+            double coeff_a{0.0};
+
+            /// Power coefficient 'B' for the effective permeability.
+            double exponent_b{0.0};
+
+            /// Power coefficient 'C' for the porosity term.
+            double exponent_c{0.0};
+
+            /// Serialisation test object.
+            static Correlation serializationTestObject();
+
+            /// Equality operator
+            ///
+            /// \param[in] other Object to which \c *this will be compared.
+            bool operator==(const Correlation& other) const;
+
+            /// Inequality operator
+            ///
+            /// \param[in] other Object to which \c *this will be compared.
+            bool operator!=(const Correlation& other) const
+            {
+                return ! (*this == other);
+            }
+
+            /// Serialisation operator
+            ///
+            /// \tparam Serializer Protocol for serialising and
+            ///   deserialising objects between memory and character
+            ///   buffers.
+            ///
+            /// \param[in,out] serializer Serialisation object.
+            template <class Serializer>
+            void serializeOp(Serializer& serializer)
+            {
+                serializer(this->coeff_a);
+                serializer(this->exponent_b);
+                serializer(this->exponent_c);
+            }
+        };
+
+        /// Serialisation test object
         static WDFAC serializationTestObject();
 
         double getDFactor(const Connection& connection, double mu, double rho, double phi) const;
-        void updateWDFAC(const DeckRecord& record);
-        //void updateWDFAC(const RestartIO::RstWell& rst_well);
-        void updateWDFACCOR(const DeckRecord& record);
-        //void updateWDFACOR(const RestartIO::RstWell& rst_well);
 
+        /// Configure D-factor calculation from well-level D-factor
+        /// description (keyword 'WDFAC')
+        ///
+        /// \param[in] record Well-level D-factor description.  Single
+        ///   record from WDFAC keyword.
+        void updateWDFAC(const DeckRecord& record);
+
+        /// Configure D-factor calculation from Dake correlation model
+        /// (keyword WDFACCOR).
+        ///
+        /// \param[in] record Dake correlation model description.  Single
+        ///   record from WDFACCOR keyword.
+        void updateWDFACCOR(const DeckRecord& record);
+
+        /// Check if any input-level connctions have a non-trivial D-factor
+        /// and update this well's D-factor category accordingly.
+        ///
+        /// \param[in] connections Connection set as defined by keyword
+        ///   COMPDAT.  This function will detect if any of the connections
+        ///   created from COMPDAT define a non-trivial D-factor at the
+        ///   connection level (item 12 of COMPDAT) and update the D-factor
+        ///   category if so.
         void updateWDFACType(const WellConnections& connections);
+
+        /// Capture sum of all CTFs for the purpose of translating
+        /// well-level D-factors to connection-level D-factors.
+        ///
+        /// \param[in] connections Connection set as defined by keyword
+        ///   COMPDAT.
         void updateTotalCF(const WellConnections& connections);
+
+        /// Retrieve current D-factor correlation model coefficients.
+        const Correlation& getDFactorCorrelationCoefficients() const
+        {
+            return this->m_corr;
+        }
+
+        /// Whether or not a flow-dependent skin factor ('D') has been
+        /// configured for the current well.
         bool useDFactor() const;
 
+        /// Equality operator
+        ///
+        /// \param[in] other Object to which \c *this will be compared.
         bool operator==(const WDFAC& other) const;
-        bool operator!=(const WDFAC& other) const;
 
-        template<class Serializer>
+        /// Inequality operator
+        ///
+        /// \param[in] other Object to which \c *this will be compared.
+        bool operator!=(const WDFAC& other) const
+        {
+            return ! (*this == other);
+        }
+
+        /// Serialisation operator
+        ///
+        /// \tparam Serializer Protocol for serialising and deserialising
+        ///   objects between memory and character buffers.
+        ///
+        /// \param[in,out] serializer Serialisation object.
+        template <class Serializer>
         void serializeOp(Serializer& serializer)
         {
-            serializer(m_a);
+            serializer(this->m_type);
-            serializer(m_b);
+            serializer(this->m_d);
-            serializer(m_c);
+            serializer(this->m_total_cf);
-            serializer(m_d);
+            serializer(this->m_corr);
-            serializer(m_total_cf);
-            serializer(m_type);
         }
 
     private:
-        double m_a{0.0};
+        /// D-factor categories.
-        double m_b{0.0};
+        enum class WDFacType
-        double m_c{0.0};
+        {
+            /// No flow-dependent skin factor is configured for this well.
+            NONE = 1,
+
+            /// Well-level D-factor.
+            DFACTOR = 2,
+
+            /// Use Dake's D-factor correlation model.
+            DAKEMODEL = 3,
+
+            /// Connection-level D-factor.
+            CON_DFACTOR = 4,
+        };
+
+        /// D-factor category for this well.
+        WDFacType m_type { WDFacType::NONE };
+
+        /// Well-level D-factor for this well.
         double m_d{0.0};
+
+	/// Total CTF sum for this well.
         double m_total_cf{-1.0};
-        WDFACTYPE m_type = WDFACTYPE::NONE;
+
+        /// Coefficients for Dake's correlation model.
+        Correlation m_corr{};
     };
 
 } // namespace Opm

src/opm/input/eclipse/Schedule/Well/WDFAC.cpp

@@ -17,126 +17,159 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
 
 #include <opm/io/eclipse/rst/well.hpp>
 #include <opm/output/eclipse/VectorItems/well.hpp>
-#include <opm/input/eclipse/Parser/ParserKeywords/W.hpp>
+
-#include <opm/input/eclipse/Deck/DeckRecord.hpp>
+#include <opm/input/eclipse/Schedule/Well/Connection.hpp>
+#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
+
 #include <opm/input/eclipse/Units/Dimension.hpp>
 #include <opm/input/eclipse/Units/UnitSystem.hpp>
 #include <opm/input/eclipse/Units/Units.hpp>
-#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
-#include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
 
+#include <opm/input/eclipse/Deck/DeckRecord.hpp>
 
-#include <string>
+#include <opm/input/eclipse/Parser/ParserKeywords/W.hpp>
-#include <vector>
+
-#include <cmath>
+#include <algorithm>
-#include <iostream>
 #include <cassert>
+#include <cmath>
 #include <numeric>
+#include <string>
 
 namespace Opm {
 
+    WDFAC::Correlation WDFAC::Correlation::serializationTestObject()
+    {
+        return { 1.23, 0.456, 0.457 };
+    }
+
+    bool WDFAC::Correlation::operator==(const Correlation& other) const
+    {
+        return (this->coeff_a == other.coeff_a)
+            && (this->exponent_b == other.exponent_b)
+            && (this->exponent_c == other.exponent_c)
+            ;
+    }
+
+    // -------------------------------------------------------------------------
+
     WDFAC WDFAC::serializationTestObject()
     {
         WDFAC result;
-        result.m_a = 1.23;
+
-        result.m_b = 0.456;
+        result.m_type = WDFacType::DAKEMODEL;
-        result.m_c = 0.457;
         result.m_d = 0.458;
         result.m_total_cf = 1.0;
-        result.m_type = WDFACTYPE::NONE;
+        result.m_corr = Correlation::serializationTestObject();
 
         return result;
     }
 
-    bool WDFAC::operator==(const WDFAC& other) const {
+    void WDFAC::updateWDFAC(const DeckRecord& record)
-        return (m_a == other.m_a)
+    {
-            && (m_b == other.m_b)
-            && (m_c == other.m_c)
-            && (m_d == other.m_d)
-            && (m_total_cf == other.m_total_cf)
-            && (m_type == other.m_type);
-    }
-
-    void WDFAC::updateWDFAC(const DeckRecord& record) {
         m_d = record.getItem<ParserKeywords::WDFAC::DFACTOR>().getSIDouble(0);
-        m_type = WDFACTYPE::DFACTOR;
+
+        m_type = WDFacType::DFACTOR;
     }
 
-    void WDFAC::updateWDFACCOR(const DeckRecord& record) {
+    void WDFAC::updateWDFACCOR(const DeckRecord& record)
-        m_a = record.getItem<ParserKeywords::WDFACCOR::A>().getSIDouble(0);
+    {
-        m_b = record.getItem<ParserKeywords::WDFACCOR::B>().getSIDouble(0);
+        this->m_corr.coeff_a = record.getItem<ParserKeywords::WDFACCOR::A>().getSIDouble(0);
-        m_c = record.getItem<ParserKeywords::WDFACCOR::C>().getSIDouble(0);
+        this->m_corr.exponent_b = record.getItem<ParserKeywords::WDFACCOR::B>().getSIDouble(0);
-        m_type = WDFACTYPE::DAKEMODEL;
+        this->m_corr.exponent_c = record.getItem<ParserKeywords::WDFACCOR::C>().getSIDouble(0);
+
+        m_type = WDFacType::DAKEMODEL;
     }
 
-    void WDFAC::updateWDFACType(const WellConnections& connections) {
+    void WDFAC::updateWDFACType(const WellConnections& connections)
-
+    {
         const auto non_trivial_dfactor =
             std::any_of(connections.begin(), connections.end(),
                 [](const auto& conn) { return conn.dFactor() != 0.0; });
 
-        // non-trivial dfactors detected use connection D factors
         if (non_trivial_dfactor) {
-            m_type = WDFACTYPE::CON_DFACTOR;
+            // Non-trivial D-factors detected.  Use connection D-factors.
-            updateTotalCF(connections);
+            m_type = WDFacType::CON_DFACTOR;
+            this->updateTotalCF(connections);
         }
     }
 
-    void WDFAC::updateTotalCF(const WellConnections& connections) {
+    void WDFAC::updateTotalCF(const WellConnections& connections)
-        m_total_cf = std::accumulate(connections.begin(), connections.end(), 0.0,
+    {
-                [](const double tot_cf, const auto& conn) { return tot_cf + conn.CF(); });
+        this->m_total_cf = std::accumulate(connections.begin(), connections.end(), 0.0,
+            [](const double tot_cf, const auto& conn) { return tot_cf + conn.CF(); });
     }
 
-    double WDFAC::getDFactor(const Connection& connection, double mu, double rho, double phi) const {
+    double WDFAC::getDFactor(const Connection& connection, double mu, double rho, double phi) const
-
+    {
-        switch (m_type)
+        switch (this->m_type) {
-        {
+        case WDFacType::NONE:
-        case WDFACTYPE::NONE:
             return 0.0;
-        case WDFACTYPE::DFACTOR: {
+
-            if (m_total_cf < 0.0) {
+        case WDFacType::DFACTOR:
-                throw std::invalid_argument { "Total connection factor is not set" };
+        {
-            }
-            return  m_d * m_total_cf / connection.CF();
-        }
-        case WDFACTYPE::CON_DFACTOR: {
-            double d =  connection.dFactor();
-            // If a negative d factor is set in COMPDAT individual connection d factors should be used directly.
-            if (d < 0)
-                return -d;
-            // If a positive d factor is set in COMPDAT the connection d factors is treated like a well d factor.
-            // and thus scaled with the connection index
             if (m_total_cf < 0.0) {
                 throw std::invalid_argument { "Total connection factor is not set" };
             }
 
-            return d * m_total_cf / connection.CF();
+            return this->m_d * this->m_total_cf / connection.CF();
         }
-        case WDFACTYPE::DAKEMODEL:
+
+        case WDFacType::DAKEMODEL:
         {
-            double Kh = connection.Kh();
+            const double Kh = connection.Kh();
-            double Ke = connection.Ke();
+            const double Ke = connection.Ke();
-            double h = Kh / Ke;
+            const double h = Kh / Ke;
-            double rw = connection.rw();
+            const double rw = connection.rw();
+
             const auto k_md = unit::convert::to(Ke, prefix::milli*unit::darcy);
-            double beta = m_a * (std::pow(k_md, m_b) * std::pow(phi, m_c));
+            double beta = m_corr.coeff_a * (std::pow(k_md, m_corr.exponent_b) * std::pow(phi, m_corr.exponent_c));
             double specific_gravity = rho / 1.225; // divide by density of air at standard conditions. 
             return beta * specific_gravity * Ke / (h * mu * rw );
         }
+        break;
+
+        case WDFacType::CON_DFACTOR:
+        {
+            const double d = connection.dFactor();
+
+            // If a negative D-factor is set in COMPDAT, then the individual
+            // connection D-factor should be used directly.
+            if (d < 0.0) {
+                return -d;
+            }
+
+            if (this->m_total_cf < 0.0) {
+                throw std::invalid_argument { "Total connection factor is not set" };
+            }
+
+            // If a positive D-factor is set in COMPDAT, then the connection
+            // D-factor is treated as a well-level D-factor and thus scaled
+            // with the connection transmissibility factor.
+            return d * m_total_cf / connection.CF();
+        }
+
         default:
             break;
         }
+
         return 0.0;
     }
 
-    bool WDFAC::useDFactor() const {
+    bool WDFAC::useDFactor() const
-        return m_type != WDFACTYPE::NONE;
+    {
+        return m_type != WDFacType::NONE;
     }
 
-    bool WDFAC::operator!=(const WDFAC& other) const {
+    bool WDFAC::operator==(const WDFAC& other) const
-        return !(*this == other);
+    {
+        return (this->m_type == other.m_type)
+            && (this->m_d == other.m_d)
+            && (this->m_total_cf == other.m_total_cf)
+            && (this->m_corr == other.m_corr)
+            ;
     }
 }

tests/test_Serialization.cpp

@@ -316,6 +316,8 @@ TEST_FOR_TYPE(VFPInjTable)
 TEST_FOR_TYPE(VFPProdTable)
 TEST_FOR_TYPE(ViscrefTable)
 TEST_FOR_TYPE(WatdentTable)
+TEST_FOR_TYPE_NAMED(WDFAC::Correlation, Correlation)
+TEST_FOR_TYPE(WDFAC)
 TEST_FOR_TYPE(Well)
 TEST_FOR_TYPE(Welldims)
 TEST_FOR_TYPE(WellBrineProperties)
@@ -336,11 +338,14 @@ TEST_FOR_TYPE(WListManager)
 TEST_FOR_TYPE(WriteRestartFileEvents)
 
 
+namespace {
+
 bool init_unit_test_func()
 {
     return true;
 }
 
+} // Anonymous namespace
 
 int main(int argc, char** argv)
 {