Use Static D-Factor Correlation Coefficient

This commit removes the porosity parameter from the getDFactor()
member function.  We then reimplement the Dake model correlation in
terms of the Connection's 'static_dfac_corr_coeff' since this
bypasses the exponentiation and, in the future, will simplify
restarting a simulation run using the WDFACCOR keyword.

While here, also split out the implementation of the connection
level D-factor to a helper function and make another helper function
for scaling a well-level D-factor to a connection-level value.

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

@@ -104,8 +104,6 @@ namespace Opm {
         /// Serialisation test object
         static WDFAC serializationTestObject();
 
-        double getDFactor(const Connection& connection, double mu, double rho, double phi) const;
-
         /// Configure D-factor calculation from well-level D-factor
         /// description (keyword 'WDFAC')
         ///
@@ -137,6 +135,20 @@ namespace Opm {
         ///   COMPDAT.
         void updateTotalCF(const WellConnections& connections);
 
+        /// Retrieve currently configured D-factor for single connection
+        ///
+        /// \param[in] rhoGS Surface condition mass density of gas
+        ///
+        /// \param[in] gas_visc Reservoir condition gas viscosity
+        ///
+        /// \param[in] conn Reservoir connection for which to retrieve the
+        ///   D-factor.
+        ///
+        /// \return D-factor for connection \p conn.
+        double getDFactor(const double      rhoGS,
+                          const double      gas_visc,
+                          const Connection& conn) const;
+
         /// Retrieve current D-factor correlation model coefficients.
         const Correlation& getDFactorCorrelationCoefficients() const
         {
@@ -203,6 +215,27 @@ namespace Opm {
 
         /// Coefficients for Dake's correlation model.
         Correlation m_corr{};
+
+        /// Retrieve connection-level D-Factor from COMPDAT entries
+        ///
+        /// Possibly translated from well-level values.
+        ///
+        /// \param[in] conn Reservoir connection for which to retrieve the
+        ///   D-factor.
+        ///
+        /// \return Connection-level D-factor.
+        double connectionLevelDFactor(const Connection& conn) const;
+
+        /// Translate well-level D-factor to connection level D-factor
+        ///
+        /// \param[in] dfac Well-level D-factor.
+        ///
+        /// \param[in] conn Reservoir connection for which to retrieve the
+        ///   D-factor.
+        ///
+        /// \return Connection-level D-factor, translated from well level.
+        double scaledWellLevelDFactor(const double      dfac,
+                                      const Connection& conn) const;
     };
 
 } // namespace Opm

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

@@ -39,6 +39,22 @@
 #include <numeric>
 #include <string>
 
+namespace {
+
+    double dakeModelDFactor(const double                          rhoGS,
+                            const double                          gas_visc,
+                            const Opm::Connection::CTFProperties& ctf_props)
+    {
+        using namespace Opm::unit;
+
+        // Specific gravity of gas relative to air at standard conditions.
+        constexpr auto rho_air = 1.22*kilogram / cubic(meter);
+        const auto specific_gravity = rhoGS / rho_air;
+
+        return ctf_props.static_dfac_corr_coeff * specific_gravity / gas_visc;
+    }
+}
+
 namespace Opm {
 
     WDFAC::Correlation WDFAC::Correlation::serializationTestObject()
@@ -103,57 +119,22 @@ namespace Opm {
             [](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 double      rhoGS,
+                             const double      gas_visc,
+                             const Connection& conn) const
     {
         switch (this->m_type) {
         case WDFacType::NONE:
             return 0.0;
 
         case WDFacType::DFACTOR:
-        {
+            return this->scaledWellLevelDFactor(this->m_d, conn);
-            if (m_total_cf < 0.0) {
-                throw std::invalid_argument { "Total connection factor is not set" };
-            }
-
-            return this->m_d * this->m_total_cf / connection.CF();
-        }
 
         case WDFacType::DAKEMODEL:
-        {
+            return dakeModelDFactor(rhoGS, gas_visc, conn.ctfProperties());
-            const double Kh = connection.Kh();
-            const double Ke = connection.Ke();
-            const double h = Kh / Ke;
-            const double rw = connection.rw();
-
-            const auto k_md = unit::convert::to(Ke, prefix::milli*unit::darcy);
-            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:
-        {
+            return this->connectionLevelDFactor(conn);
-            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;
@@ -172,4 +153,32 @@ namespace Opm {
             && (this->m_corr == other.m_corr)
             ;
     }
+
+    double WDFAC::connectionLevelDFactor(const Connection& conn) const
+    {
+        const double d = conn.dFactor();
+
+        // Negative D-factor values in COMPDAT should be used directly as
+        // connection-level D-factors.
+        if (d < 0.0) {
+            return -d;
+        }
+
+        // Positive D-factor values in COMPDAT are treated as well-level
+        // values and scaled with the CTF for translation to connection
+        // level.
+        return this->scaledWellLevelDFactor(d, conn);
+    }
+
+    double WDFAC::scaledWellLevelDFactor(const double      dfac,
+                                         const Connection& conn) const
+    {
+        if (this->m_total_cf < 0.0) {
+            throw std::invalid_argument {
+                "Total well-level connection factor is not set"
+            };
+        }
+
+        return dfac * this->m_total_cf / conn.CF();
+    }
 }

tests/parser/ScheduleTests.cpp

@@ -5866,7 +5866,6 @@ END
 
     const double rho = 1.0;
     const double mu = 0.01*prefix::centi*unit::Poise;
-    const double phi = 0.3;
 
     {
         const auto& well11 = sched.getWell("W1", 1);
@@ -5879,8 +5878,8 @@ END
                             R"(Well "W1" must use D-Factors at step 1)");
 
         // Well-level D-factor scaled by connection transmissibility factor.
-        BOOST_CHECK_CLOSE(wdfac11.getDFactor(well11.getConnections()[0], mu, rho, phi), 6*1.0*dFacUnit, 1e-12);
+        BOOST_CHECK_CLOSE(wdfac11.getDFactor(rho, mu, well11.getConnections()[0]), 6*1.0*dFacUnit, 1e-12);
-        BOOST_CHECK_CLOSE(wdfac21.getDFactor(well21.getConnections()[0], mu, rho, phi),   2.0*dFacUnit, 1e-12);
+        BOOST_CHECK_CLOSE(wdfac21.getDFactor(rho, mu, well21.getConnections()[0]),   2.0*dFacUnit, 1e-12);
     }
 
     {
@@ -5889,8 +5888,29 @@ END
         const auto& wdfac12 = well12.getWDFAC();
         const auto& wdfac22 = well22.getWDFAC();
 
-        BOOST_CHECK_CLOSE(wdfac12.getDFactor(well12.getConnections()[0], mu, rho, phi), 5.19e-1, 3);
+        // Intentional copy.
-        BOOST_CHECK_CLOSE(wdfac22.getDFactor(well22.getConnections()[0], mu, rho, phi), 2.0*dFacUnit, 1e-12);
+        auto conn0 = well12.getConnections()[0];
+        {
+            const auto& corr = wdfac12.getDFactorCorrelationCoefficients();
+
+            const auto k = 10.0*prefix::milli*unit::darcy;
+            const auto k0 = 1.0*prefix::milli*unit::darcy;
+            const auto h = 20.0*unit::meter;
+            const auto rw = 0.108*unit::meter;
+            const auto poro = 0.3;
+
+            const auto static_dfac_corr_coeff = corr.coeff_a
+                * std::pow(k / k0, corr.exponent_b)
+                * std::pow(poro  , corr.exponent_c)
+                * k / (h * rw);
+
+            BOOST_CHECK_CLOSE(static_dfac_corr_coeff, 6.238808556951547e-06, 1.0e-8);
+
+            conn0.setStaticDFacCorrCoeff(static_dfac_corr_coeff);
+        }
+
+        BOOST_CHECK_CLOSE(wdfac12.getDFactor(rho, mu, conn0), 5.19e-1, 3);
+        BOOST_CHECK_CLOSE(wdfac22.getDFactor(rho, mu, conn0), 2.0*dFacUnit, 1e-12);
     }
 
     {
@@ -5906,8 +5926,8 @@ END
         BOOST_CHECK_CLOSE(well13.getConnections()[1].dFactor(),  0.0*dFacUnit, 1e-12);
         BOOST_CHECK_CLOSE(well13.getConnections()[2].dFactor(), 11.0*dFacUnit, 1e-12);
 
-        BOOST_CHECK_CLOSE(wdfac13.getDFactor(well13.getConnections()[2], mu, rho, phi), 6.0/3.0*11.0*dFacUnit, 1e-12);
+        BOOST_CHECK_CLOSE(wdfac13.getDFactor(rho, mu, well13.getConnections()[2]), 6.0/3.0*11.0*dFacUnit, 1e-12);
-        BOOST_CHECK_CLOSE(wdfac23.getDFactor(well23.getConnections()[0], mu, rho, phi),          2.0*dFacUnit, 1e-12);
+        BOOST_CHECK_CLOSE(wdfac23.getDFactor(rho, mu, well23.getConnections()[0]),          2.0*dFacUnit, 1e-12);
     }
 }