refactoring wdfac getDfactor

This commit is contained in:
Tor Harald Sandve 2023-11-22 09:17:59 +01:00
parent 98220372eb
commit e362e1018d
3 changed files with 57 additions and 39 deletions

View File

@ -43,18 +43,14 @@ namespace Opm {
public:
static WDFAC serializationTestObject();
double getDFactor(double rho, double mu, double k, double phi, double rw, double h) const;
double getDFactor(const Connection& connection, double mu, double rho, double phi, double trans_mult) 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);
void updateWDFACType(const WellConnections& connections);
bool useDFactor() const;
bool useConnectionDFactor() const;
bool useWellDFactor() const;
WDFACTYPE getType() const;
bool operator==(const WDFAC& other) const;
bool operator!=(const WDFAC& other) const;
@ -66,6 +62,7 @@ namespace Opm {
serializer(m_b);
serializer(m_c);
serializer(m_d);
serializer(m_total_cf);
serializer(m_type);
}
@ -74,6 +71,7 @@ namespace Opm {
double m_b{0.0};
double m_c{0.0};
double m_d{0.0};
double m_total_cf{-1.0};
WDFACTYPE m_type = WDFACTYPE::NONE;
};

View File

@ -34,6 +34,7 @@
#include <cmath>
#include <iostream>
#include <cassert>
#include <numeric>
namespace Opm {
@ -44,6 +45,7 @@ namespace Opm {
result.m_b = 0.456;
result.m_c = 0.457;
result.m_d = 0.458;
result.m_total_cf = 1.0;
result.m_type = WDFACTYPE::NONE;
return result;
@ -54,6 +56,7 @@ namespace Opm {
&& (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);
}
@ -79,26 +82,43 @@ namespace Opm {
if (non_trivial_dfactor) {
m_type = WDFACTYPE::CON_DFACTOR;
}
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, double trans_mult) const {
double WDFAC::getDFactor(double rho, double mu, double k, double phi, double rw, double h) const {
if (m_total_cf < 0.0) {
throw std::invalid_argument { "Total connection factor is not set" };
}
switch (m_type)
{
case WDFACTYPE::NONE:
return 0.0;
case WDFACTYPE::DFACTOR:
return m_d;
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
return d * m_total_cf / connection.CF();
}
case WDFACTYPE::DAKEMODEL:
{
const auto k_md = unit::convert::to(k, prefix::milli*unit::darcy);
double Kh = connection.Kh() * trans_mult;
double Ke = connection.Ke() * trans_mult;
double h = Kh / Ke;
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 specific_gravity = rho / 1.225; // divide by density of air at standard conditions.
return beta * specific_gravity * k / (h * mu * rw );
return beta * specific_gravity * Ke / (h * mu * rw );
}
case WDFACTYPE::CON_DFACTOR:
throw std::invalid_argument { "getDFactor should not be called if D factor is given by COMPDAT item 12." };
default:
break;
}
@ -109,9 +129,6 @@ namespace Opm {
return m_type != WDFACTYPE::NONE;
}
WDFACTYPE WDFAC::getType() const {
return m_type;
}
bool WDFAC::operator!=(const WDFAC& other) const {
return !(*this == other);
}

View File

@ -5711,13 +5711,13 @@ DEPTHZ
121*2000.0 /
PORO
1000*0.1 /
1000*0.3 /
PERMX
1000*1 /
1000*10 /
PERMY
1000*0.1 /
1000*10 /
PERMZ
1000*0.01 /
1000*10 /
SCHEDULE
@ -5730,9 +5730,10 @@ WELSPECS
/
COMPDAT
'W1' 3 3 1 1 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 /
'W1' 3 3 2 2 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 /
'W1' 3 3 3 3 'OPEN' 1* 46.825 0.311 4332.346 1* 1* 'X' 22.123 /
'W1' 3 3 1 1 'OPEN' 1* 1 0.216 200 1* 1* 'X' /
'W1' 3 3 2 2 'OPEN' 1* 2 0.216 200 1* 1* 'X' /
'W1' 3 3 3 3 'OPEN' 1* 3 0.216 200 1* 1* 'X' /
'W2' 3 3 3 3 'OPEN' 1* 1 0.216 200 1* 11 'X' /
/
WDFAC
@ -5744,6 +5745,13 @@ DATES -- 2
10 NOV 2008 /
/
COMPDAT
'W1' 3 3 1 1 'OPEN' 1* 1* 0.216 200 1* 1* 'X' /
'W1' 3 3 2 2 'OPEN' 1* 1* 0.216 200 1* 1* 'X' /
'W1' 3 3 3 3 'OPEN' 1* 1* 0.216 200 1* 1* 'X' /
'W2' 3 3 3 3 'OPEN' 1* 1 0.216 200 1* 11 'X' /
/
WDFACCOR
-- 'W1' 8.957e10 1.1045 0.0 /
'W1' 1.984e-7 -1.1045 0.0 /
@ -5754,9 +5762,10 @@ DATES -- 3
/
COMPDAT
'W1' 3 3 1 1 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 /
'W1' 3 3 2 2 'OPEN' 1* 32.948 0.311 3047.839 1* 0 'X' 22.100 /
'W1' 3 3 3 3 'OPEN' 1* 46.825 0.311 4332.346 1* 11 'X' 22.123 /
'W1' 3 3 1 1 'OPEN' 1* 1 0.216 200 1* 1* 'X' /
'W1' 3 3 2 2 'OPEN' 1* 2 0.216 200 1* 0 'X' /
'W1' 3 3 3 3 'OPEN' 1* 3 0.216 200 1* 11 'X' /
'W2' 3 3 3 3 'OPEN' 1* 1 0.216 200 1* 11 'X' /
/
@ -5774,43 +5783,37 @@ END
double rho = 1.0;
double mu = 0.01*Opm::prefix::centi * Opm::unit::Poise;
double k = 10.0 * Opm::prefix::milli * Opm::unit::darcy;
double h = 20;
double rw = 0.108;
double phi = 0.3;
double trans_mult = 1.0;
// WDFAC overwrites D factor in COMDAT
BOOST_CHECK(wdfac11.getType() != WDFACTYPE::CON_DFACTOR);
BOOST_CHECK(wdfac11.useDFactor());
BOOST_CHECK(wdfac21.getType() != WDFACTYPE::CON_DFACTOR);
BOOST_CHECK_CLOSE(wdfac11.getDFactor(rho, mu, k, phi, rw, h), 1*Opm::unit::day, 1e-12);
BOOST_CHECK_CLOSE(wdfac21.getDFactor(rho, mu, k, phi, rw, h), 2*Opm::unit::day, 1e-12);
// well d factor scaled by connection CF.
BOOST_CHECK_CLOSE(wdfac11.getDFactor(well11.getConnections()[0], mu, rho, phi, trans_mult), 6*1*Opm::unit::day, 1e-12);
BOOST_CHECK_CLOSE(wdfac21.getDFactor(well21.getConnections()[0], mu, rho, phi, trans_mult), 2*Opm::unit::day, 1e-12);
const auto& well12 = sched.getWell("W1", 2);
const auto& well22 = sched.getWell("W2", 2);
const auto& wdfac12 = well12.getWDFAC();
const auto& wdfac22 = well22.getWDFAC();
BOOST_CHECK_CLOSE(wdfac12.getDFactor(rho, mu, k, phi, rw, h), 5.19e-1, 3);
BOOST_CHECK_CLOSE(wdfac22.getDFactor(rho, mu, k, phi, rw, h), 2*Opm::unit::day, 1e-12);
BOOST_CHECK_CLOSE(wdfac12.getDFactor(well12.getConnections()[0], mu, rho, phi, trans_mult), 5.19e-1, 3);
BOOST_CHECK_CLOSE(wdfac22.getDFactor(well22.getConnections()[0], mu, rho, phi, trans_mult), 2*Opm::unit::day, 1e-12);
const auto& well13 = sched.getWell("W1", 3);
const auto& well23 = sched.getWell("W2", 3);
const auto& wdfac13 = well13.getWDFAC();
const auto& wdfac23 = well23.getWDFAC();
BOOST_CHECK(wdfac13.getType() == WDFACTYPE::CON_DFACTOR);
BOOST_CHECK(wdfac13.useDFactor());
BOOST_CHECK(wdfac23.getType() != WDFACTYPE::CON_DFACTOR);
BOOST_CHECK_CLOSE(well13.getConnections()[0].dFactor(), 0*Opm::unit::day, 1e-12);
BOOST_CHECK_CLOSE(well13.getConnections()[1].dFactor(), 0*Opm::unit::day, 1e-12);
BOOST_CHECK_CLOSE(well13.getConnections()[2].dFactor(), 11*Opm::unit::day, 1e-12);
BOOST_CHECK_CLOSE(wdfac23.getDFactor(rho, mu, k, phi, rw, h), 2*Opm::unit::day, 1e-12);
BOOST_CHECK_THROW(wdfac13.getDFactor(rho, mu, k, phi, rw, h ), std::exception);
BOOST_CHECK_CLOSE(wdfac13.getDFactor(well13.getConnections()[2], mu, rho, phi, trans_mult), 6/3*11*Opm::unit::day, 1e-12);
BOOST_CHECK_CLOSE(wdfac23.getDFactor(well23.getConnections()[0], mu, rho, phi, trans_mult), 2*Opm::unit::day, 1e-12);
}
BOOST_AUTO_TEST_CASE(createDeckWithBC) {