From e0fe776e2eca6c65134ea36820075728d5b2c99b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 10 Nov 2023 15:47:46 +0100 Subject: [PATCH] Move Static Dake Model Calculation to Connection Class Certain parts of the Dake model D-factor correlation relation depend only on static parameters at the cell level. Calculate the static product once, and store it as a member of Connection::CTFProperties. In particular, the new member double CTFProperties::static_dfac_corr_coeff stores the product A * (K/K0)**B * poro**C * Ke/(h*rw) with 'A', 'B', and 'C', being the Dake model parameters. We compute new values for the 'static_dfac_corr_coeff' as part of processing the WDFACCOR keyword. This static product will also be output to the restart file in follow-up commits, and will be used to simplify evaluating the full dynamic Dake model correlation. --- .../eclipse/Schedule/Well/Connection.hpp | 6 + .../eclipse/Schedule/Well/WellConnections.hpp | 5 + .../eclipse/Schedule/KeywordHandlers.cpp | 131 +++++++++++++----- .../eclipse/Schedule/Well/Connection.cpp | 9 +- .../eclipse/Schedule/Well/WellConnections.cpp | 35 +++++ tests/parser/ConnectionTests.cpp | 4 +- 6 files changed, 152 insertions(+), 38 deletions(-) diff --git a/opm/input/eclipse/Schedule/Well/Connection.hpp b/opm/input/eclipse/Schedule/Well/Connection.hpp index 4be15fbf2..293d4c11a 100644 --- a/opm/input/eclipse/Schedule/Well/Connection.hpp +++ b/opm/input/eclipse/Schedule/Well/Connection.hpp @@ -119,6 +119,10 @@ namespace Opm { /// for gas. double d_factor{}; + /// Product of certain static elements of D-factor correlation + /// law (WDFACCOR keyword). + double static_dfac_corr_coeff{}; + /// Denominator in peaceman's formula-i.e., log(r0/rw) + skin. double peaceman_denom{}; @@ -157,6 +161,7 @@ namespace Opm { serializer(this->connection_length); serializer(this->skin_factor); serializer(this->d_factor); + serializer(this->static_dfac_corr_coeff); serializer(this->peaceman_denom); } }; @@ -241,6 +246,7 @@ namespace Opm { void setKe(double Ke); void setCF(double CF); void setDefaultSatTabId(bool id); + void setStaticDFacCorrCoeff(const double c); void scaleWellPi(double wellPi); bool prepareWellPIScaling(); diff --git a/opm/input/eclipse/Schedule/Well/WellConnections.hpp b/opm/input/eclipse/Schedule/Well/WellConnections.hpp index 052a9054e..73939323b 100644 --- a/opm/input/eclipse/Schedule/Well/WellConnections.hpp +++ b/opm/input/eclipse/Schedule/Well/WellConnections.hpp @@ -39,6 +39,7 @@ namespace Opm { class FieldPropsManager; class KeywordLocation; class ScheduleGrid; + class WDFAC; } // namespace Opm namespace Opm { @@ -88,6 +89,7 @@ namespace Opm { void loadCOMPDAT(const DeckRecord& record, const ScheduleGrid& grid, const std::string& wname, + const WDFAC& wdfac, const KeywordLocation& location); void loadCOMPTRAJ(const DeckRecord& record, @@ -101,6 +103,9 @@ namespace Opm { const std::string& wname, const KeywordLocation& location); + void applyDFactorCorrelation(const ScheduleGrid& grid, + const WDFAC& wdfac); + int getHeadI() const; int getHeadJ() const; const std::vector& getMD() const; diff --git a/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp b/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp index e475ce6f6..98d83d76c 100644 --- a/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp +++ b/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp @@ -196,37 +196,51 @@ namespace { this->snapshots.back().network.update( std::move( ext_network )); } - void Schedule::handleCOMPDAT(HandlerContext& handlerContext) { + void Schedule::handleCOMPDAT(HandlerContext& handlerContext) { std::unordered_set wells; for (const auto& record : handlerContext.keyword) { - const std::string& wellNamePattern = record.getItem("WELL").getTrimmedString(0); - auto wellnames = this->wellNames(wellNamePattern, handlerContext, - isWList(handlerContext.currentStep, - wellNamePattern)); + const auto wellNamePattern = record.getItem("WELL").getTrimmedString(0); + const auto wellnames = + this->wellNames(wellNamePattern, handlerContext, + isWList(handlerContext.currentStep, wellNamePattern)); for (const auto& name : wellnames) { auto well2 = this->snapshots.back().wells.get(name); - auto connections = std::shared_ptr( new WellConnections( well2.getConnections())); - connections->loadCOMPDAT(record, handlerContext.grid, name, handlerContext.keyword.location()); - if (well2.updateConnections(connections, handlerContext.grid)) { + auto connections = std::make_shared(well2.getConnections()); + const auto origWellConnSetIsEmpty = connections->empty(); + + connections->loadCOMPDAT(record, handlerContext.grid, name, + well2.getWDFAC(), handlerContext.keyword.location()); + const auto newWellConnSetIsEmpty = connections->empty(); + + if (well2.updateConnections(std::move(connections), handlerContext.grid)) { auto wdfac = std::make_shared(well2.getWDFAC()); - wdfac->updateWDFACType(*connections); + wdfac->updateWDFACType(well2.getConnections()); + well2.updateWDFAC(std::move(wdfac)); - this->snapshots.back().wells.update( well2 ); - wells.insert( name ); + this->snapshots.back().wells.update(well2); + + wells.insert(name); } - if (connections->empty() && well2.getConnections().empty()) { + if (origWellConnSetIsEmpty && newWellConnSetIsEmpty) { const auto& location = handlerContext.keyword.location(); - auto msg = fmt::format("Problem with COMPDAT/{}\n" - "In {} line {}\n" - "Well {} is not connected to grid - will remain SHUT", name, location.filename, location.lineno, name); + + const auto msg = fmt::format(R"(Problem with COMPDAT/{} +In {} line {} +Well {} is not connected to grid - will remain SHUT)", + name, location.filename, + location.lineno, name); + OpmLog::warning(msg); } - this->snapshots.back().wellgroup_events().addEvent( name, ScheduleEvents::COMPLETION_CHANGE); + + this->snapshots.back().wellgroup_events() + .addEvent(name, ScheduleEvents::COMPLETION_CHANGE); } } + this->snapshots.back().events().addEvent(ScheduleEvents::COMPLETION_CHANGE); // In the case the wells reference depth has been defaulted in the @@ -234,9 +248,10 @@ namespace { // reference depth exactly when the COMPDAT keyword has been completely // processed. for (const auto& wname : wells) { - auto& well = this->snapshots.back().wells.get( wname ); + auto well = this->snapshots.back().wells.get(wname); well.updateRefDepth(); - this->snapshots.back().wells.update( std::move(well)); + + this->snapshots.back().wells.update(std::move(well)); } if (! wells.empty()) { @@ -268,20 +283,30 @@ namespace { this->snapshots.back().events().addEvent(ScheduleEvents::COMPLETION_CHANGE); } - void Schedule::handleCOMPTRAJ(HandlerContext& handlerContext) { + void Schedule::handleCOMPTRAJ(HandlerContext& handlerContext) + { // Keyword WELTRAJ must be read first std::unordered_set wells; external::cvf::ref cellSearchTree = nullptr; + for (const auto& record : handlerContext.keyword) { - const std::string& wellNamePattern = record.getItem("WELL").getTrimmedString(0); - auto wellnames = this->wellNames(wellNamePattern, handlerContext ); + const auto wellNamePattern = record.getItem("WELL").getTrimmedString(0); + const auto wellnames = this->wellNames(wellNamePattern, handlerContext); for (const auto& name : wellnames) { auto well2 = this->snapshots.back().wells.get(name); - auto connections = std::make_shared(WellConnections(well2.getConnections())); - // cellsearchTree is calculated only once and is used to calculated cell intersections of the perforations specified in COMPTRAJ - connections->loadCOMPTRAJ(record, handlerContext.grid, name, handlerContext.keyword.location(), cellSearchTree); - // In the case that defaults are used in WELSPECS for headI/J the headI/J are calculated based on the well trajectory data + auto connections = std::make_shared(well2.getConnections()); + + // cellsearchTree is calculated only once and is used to + // calculated cell intersections of the perforations + // specified in COMPTRAJ + connections->loadCOMPTRAJ(record, handlerContext.grid, name, + handlerContext.keyword.location(), + cellSearchTree); + + // In the case that defaults are used in WELSPECS for + // headI/J the headI/J are calculated based on the well + // trajectory data well2.updateHead(connections->getHeadI(), connections->getHeadJ()); if (well2.updateConnections(connections, handlerContext.grid)) { this->snapshots.back().wells.update( well2 ); @@ -295,19 +320,23 @@ namespace { "Well {} is not connected to grid - will remain SHUT", name, location.filename, location.lineno, name); OpmLog::warning(msg); } - this->snapshots.back().wellgroup_events().addEvent( name, ScheduleEvents::COMPLETION_CHANGE); + + this->snapshots.back().wellgroup_events() + .addEvent(name, ScheduleEvents::COMPLETION_CHANGE); } } + this->snapshots.back().events().addEvent(ScheduleEvents::COMPLETION_CHANGE); // In the case the wells reference depth has been defaulted in the // WELSPECS keyword we need to force a calculation of the wells - // reference depth exactly when the WELCOML keyword has been completely - // processed. + // reference depth exactly when the WELCOML keyword has been + // completely processed. for (const auto& wname : wells) { - auto& well = this->snapshots.back().wells.get( wname ); + auto& well = this->snapshots.back().wells.get(wname); well.updateRefDepth(); - this->snapshots.back().wells.update( std::move(well)); + + this->snapshots.back().wells.update(std::move(well)); } if (! wells.empty()) { @@ -1641,25 +1670,55 @@ File {} line {}.)", wname, location.keyword, location.filename, location.lineno) auto wdfac = std::make_shared(well.getWDFAC()); wdfac->updateWDFAC( record ); wdfac->updateTotalCF(well.getConnections()); - if (well.updateWDFAC(std::move(wdfac))) - this->snapshots.back().wells.update( std::move(well) ); + + if (well.updateWDFAC(std::move(wdfac))) { + this->snapshots.back().wells.update(std::move(well)); + + handlerContext.affected_well(well_name); + handlerContext.record_well_structure_change(); + + this->snapshots.back().events() + .addEvent(ScheduleEvents::COMPLETION_CHANGE); + + this->snapshots.back().wellgroup_events() + .addEvent(well_name, ScheduleEvents::COMPLETION_CHANGE); + } } } } void Schedule::handleWDFACCOR(HandlerContext& handlerContext) { for (const auto& record : handlerContext.keyword) { - const std::string& wellNamePattern = record.getItem("WELLNAME").getTrimmedString(0); + const std::string wellNamePattern = record.getItem("WELLNAME").getTrimmedString(0); const auto well_names = wellNames(wellNamePattern, handlerContext.currentStep); if (well_names.empty()) this->invalidNamePattern(wellNamePattern, handlerContext); for (const auto& well_name : well_names) { auto well = this->snapshots.back().wells.get(well_name); + auto conns = std::make_shared(well.getConnections()); + auto wdfac = std::make_shared(well.getWDFAC()); - wdfac->updateWDFACCOR( record ); - if (well.updateWDFAC(std::move(wdfac))) - this->snapshots.back().wells.update( std::move(well) ); + wdfac->updateWDFACCOR(record); + + conns->applyDFactorCorrelation(handlerContext.grid, *wdfac); + const auto updateConns = + well.updateConnections(std::move(conns), handlerContext.grid); + + if (well.updateWDFAC(std::move(wdfac)) || updateConns) { + this->snapshots.back().wells.update(std::move(well)); + + handlerContext.affected_well(well_name); + handlerContext.record_well_structure_change(); + + if (updateConns) { + this->snapshots.back().events() + .addEvent(ScheduleEvents::COMPLETION_CHANGE); + + this->snapshots.back().wellgroup_events() + .addEvent(well_name, ScheduleEvents::COMPLETION_CHANGE); + } + } } } } diff --git a/src/opm/input/eclipse/Schedule/Well/Connection.cpp b/src/opm/input/eclipse/Schedule/Well/Connection.cpp index 58673416b..be7712392 100644 --- a/src/opm/input/eclipse/Schedule/Well/Connection.cpp +++ b/src/opm/input/eclipse/Schedule/Well/Connection.cpp @@ -78,7 +78,8 @@ namespace Opm props.connection_length = 7.0; props.skin_factor = 8.0; props.d_factor = 9.0; - props.peaceman_denom = 10.0; + props.static_dfac_corr_coeff = 10.0; + props.peaceman_denom = 11.0; return props; } @@ -94,6 +95,7 @@ namespace Opm && (this->connection_length == that.connection_length) && (this->skin_factor == that.skin_factor) && (this->d_factor == that.d_factor) + && (this->static_dfac_corr_coeff == that.static_dfac_corr_coeff) && (this->peaceman_denom == that.peaceman_denom) ; } @@ -394,6 +396,11 @@ namespace Opm return true; } + void Connection::setStaticDFacCorrCoeff(const double c) + { + this->ctf_properties_.static_dfac_corr_coeff = c; + } + std::string Connection::str() const { std::stringstream ss; diff --git a/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp b/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp index 38ffcb58b..d6153a7e8 100644 --- a/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp +++ b/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp @@ -153,6 +153,22 @@ namespace { return peacemanDenominator(ctf_props.r0, ctf_props.rw, ctf_props.skin_factor); } + double staticForchheimerCoeffcient(const Opm::Connection::CTFProperties& ctf_props, + const double porosity, + const Opm::WDFAC& wdfac) + { + // Reference/background permeability against which to scale cell + // permeability for model evaluation. + constexpr auto Kref = 1.0*Opm::prefix::milli*Opm::unit::darcy; + + const auto& corr_coeff = wdfac.getDFactorCorrelationCoefficients(); + + return corr_coeff.coeff_a + * std::pow(ctf_props.Ke / Kref, corr_coeff.exponent_b) + * std::pow(porosity , corr_coeff.exponent_c) + * ctf_props.Ke / (ctf_props.connection_length * ctf_props.rw); + } + // Calculate permeability thickness Kh for line segment in a cell for x,y,z directions std::array permThickness(const external::cvf::Vec3d& effective_connection, @@ -339,6 +355,7 @@ namespace Opm { void WellConnections::loadCOMPDAT(const DeckRecord& record, const ScheduleGrid& grid, const std::string& wname, + const WDFAC& wdfac, const KeywordLocation& location) { const auto& itemI = record.getItem("I"); @@ -492,6 +509,9 @@ The cell ({},{},{}) in well {} is not active and the connection will be ignored) // PolymerMW module. ctf_props.re = std::sqrt(D[0] * D[1] / angle * 2); + ctf_props.static_dfac_corr_coeff = + staticForchheimerCoeffcient(ctf_props, props->poro, wdfac); + auto prev = std::find_if(this->m_connections.begin(), this->m_connections.end(), [I, J, k](const Connection& c) @@ -734,6 +754,21 @@ CF and Kh items for well {} must both be specified or both defaulted/negative)", this->md.push_back(record.getItem("MD").getSIDouble(0)); } + void WellConnections::applyDFactorCorrelation(const ScheduleGrid& grid, + const WDFAC& wdfac) + { + for (auto& conn : this->m_connections) { + const auto& complCell = grid.get_cell(conn.getI(), conn.getJ(), conn.getK()); + if (! complCell.is_active()) { + continue; + } + + conn.setStaticDFacCorrCoeff + (staticForchheimerCoeffcient(conn.ctfProperties(), + complCell.props->poro, wdfac)); + } + } + std::size_t WellConnections::size() const { return m_connections.size(); diff --git a/tests/parser/ConnectionTests.cpp b/tests/parser/ConnectionTests.cpp index 7e32848e6..ae1566b7c 100644 --- a/tests/parser/ConnectionTests.cpp +++ b/tests/parser/ConnectionTests.cpp @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -65,6 +66,7 @@ namespace { }; const auto deck = Opm::Parser{}.parseString(compdat_keyword); + const auto wdfac = Opm::WDFAC{}; const auto loc = Opm::KeywordLocation{}; const Opm::EclipseGrid grid { 10, 10, 10 }; @@ -77,7 +79,7 @@ namespace { const auto sg = Opm::ScheduleGrid { grid, field_props, completed_cells }; for (const auto& rec : deck["COMPDAT"][0]) { - connections.loadCOMPDAT(rec, sg, "WELL", loc); + connections.loadCOMPDAT(rec, sg, "WELL", wdfac, loc); } return connections;