Implement Forchheimer term in wellIndex

Add output of CDFAC

Add effect of compaction on CTFAC

ebos/eclgenericvanguard.cc

@@ -53,6 +53,7 @@
 #include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
 #include <opm/input/eclipse/Schedule/UDQ/UDQState.hpp>
 #include <opm/input/eclipse/Schedule/Well/NameOrder.hpp>
+#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
 #include <opm/input/eclipse/Schedule/Well/Well.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellBrineProperties.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>

opm/simulators/utils/ParallelSerialization.cpp

@@ -46,6 +46,7 @@
 #include <opm/input/eclipse/Schedule/UDQ/UDQASTNode.hpp>
 #include <opm/input/eclipse/Schedule/UDQ/UDQConfig.hpp>
 #include <opm/input/eclipse/Schedule/Well/NameOrder.hpp>
+#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellConnections.hpp>
 #include <opm/input/eclipse/Schedule/Well/Well.hpp>
 #include <opm/input/eclipse/Schedule/Well/WellBrineProperties.hpp>
@@ -61,6 +62,7 @@
 #include <opm/input/eclipse/Schedule/Well/WVFPDP.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPEXP.hpp>
 
+
 #include <ebos/eclmpiserializer.hh>
 
 #include <dune/common/parallel/mpihelper.hh>

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -685,8 +685,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"WCUTBACK", {true, std::nullopt}},
         {"WCUTBACT", {true, std::nullopt}},
         {"WCYCLE", {true, std::nullopt}},
-        {"WDFACCOR", {true, std::nullopt}},
-        {"WDFAC", {true, std::nullopt}},
         {"WDRILTIM", {true, std::nullopt}},
         {"WDRILPRI", {true, std::nullopt}},
         {"WDRILRES", {true, std::nullopt}},

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -853,6 +853,7 @@ assignShutConnections(data::Wells& wsrpt,
 
             xc.effective_Kh = conn.Kh();
             xc.trans_factor = conn.CF();
+            xc.d_factor = conn.dFactor();
         }
 
         ++wellID;

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -133,6 +133,14 @@ public:
         return this->active_wgstate_.well_state;
     }
 
+    /*
+      Will return the currently active nupcolWellState; must initialize
+      the internal nupcol wellstate with initNupcolWellState() first.
+    */
+    const WellState& nupcolWellState() const
+    {
+        return this->nupcol_wgstate_.well_state;
+    }
     GroupState& groupState() { return this->active_wgstate_.group_state; }
 
     WellTestState& wellTestState() { return this->active_wgstate_.well_test_state; }
@@ -276,18 +284,13 @@ protected:
         return this->last_valid_wgstate_.well_state;
     }
 
+
     const WGState& prevWGState() const
     {
         return this->last_valid_wgstate_;
     }
-    /*
+
-      Will return the currently active nupcolWellState; must initialize
+
-      the internal nupcol wellstate with initNupcolWellState() first.
-    */
-    const WellState& nupcolWellState() const
-    {
-        return this->nupcol_wgstate_.well_state;
-    }
 
     /*
       Will store a copy of the input argument well_state in the

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -597,6 +597,11 @@ namespace Opm {
                 well->updateWaterThroughput(dt, this->wellState());
             }
         }
+        // update connection transmissibility factor and d factor (if applicable) in the wellstate
+        for (const auto& well : well_container_) {
+            well->updateConnectionTransmissibilityFactor(ebosSimulator_, this->wellState().well(well->indexOfWell()));
+            well->updateConnectionDFactor(ebosSimulator_, this->wellState().well(well->indexOfWell()));
+        }
 
         if (Indices::waterEnabled) {
             this->updateFiltrationParticleVolume(dt, FluidSystem::waterPhaseIdx);

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -347,9 +347,9 @@ namespace Opm
                 // flux for each perforation
                 std::vector<Scalar> mob(this->num_components_, 0.);
                 getMobility(ebosSimulator, perf, mob, deferred_logger);
-                double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
+                const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
-                const double Tw = this->well_index_[perf] * trans_mult;
+                const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
-
+                const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
                 const Scalar seg_pressure = segment_pressure[seg];
                 std::vector<Scalar> cq_s(this->num_components_, 0.);
                 Scalar perf_press = 0.0;
@@ -1180,8 +1180,9 @@ namespace Opm
                 }
 
                 // the well index associated with the connection
-                const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
+                const double trans_mult = ebos_simulator.problem().template rockCompTransMultiplier<double>(int_quantities, cell_idx);
-
+                const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+                const double tw_perf = this->wellIndex(perf, int_quantities, trans_mult, wellstate_nupcol);  
                 std::vector<double> ipr_a_perf(this->ipr_a_.size());
                 std::vector<double> ipr_b_perf(this->ipr_b_.size());
                 for (int comp_idx = 0; comp_idx < this->num_components_; ++comp_idx) {
@@ -1680,7 +1681,8 @@ namespace Opm
                 std::vector<EvalWell> mob(this->num_components_, 0.0);
                 getMobility(ebosSimulator, perf, mob, deferred_logger);
                 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
-                const double Tw = this->well_index_[perf] * trans_mult;
+                const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+                const double Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
                 std::vector<EvalWell> cq_s(this->num_components_, 0.0);
                 EvalWell perf_press;
                 PerforationRates perfRates;
@@ -1992,7 +1994,8 @@ namespace Opm
                 std::vector<Scalar> mob(this->num_components_, 0.0);
                 getMobility(ebosSimulator, perf, mob, deferred_logger);
                 const double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(int_quants, cell_idx);
-                const double Tw = this->well_index_[perf] * trans_mult;
+                const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+                const double Tw = this->wellIndex(perf, int_quants, trans_mult, wellstate_nupcol);
                 std::vector<Scalar> cq_s(this->num_components_, 0.0);
                 Scalar perf_press = 0.0;
                 PerforationRates perf_rates;

opm/simulators/wells/PerfData.cpp

@@ -40,6 +40,7 @@ PerfData::PerfData(std::size_t num_perf, double pressure_first_connection_, bool
     , micp_rates(num_perf)
     , cell_index(num_perf)
     , connection_transmissibility_factor(num_perf)
+    , connection_d_factor(num_perf)
     , satnum_id(num_perf)
     , ecl_index(num_perf)
 {
@@ -65,6 +66,7 @@ PerfData PerfData::serializationTestObject()
     result.micp_rates = {16.0};
     result.cell_index = {17, 18, 19, 20};
     result.connection_transmissibility_factor = {21.0};
+    result.connection_d_factor = {21.5};
     result.satnum_id = {22, 23};
     result.ecl_index = {24};
     result.water_throughput = {25.0, 26.0};
@@ -119,6 +121,7 @@ bool PerfData::operator==(const PerfData& rhs) const
            this->micp_rates == rhs.micp_rates &&
            this->cell_index == rhs.cell_index &&
            this->connection_transmissibility_factor == rhs.connection_transmissibility_factor &&
+           this->connection_d_factor == rhs.connection_d_factor &&
            this->satnum_id == rhs.satnum_id &&
            this->ecl_index == rhs.ecl_index &&
            this->water_throughput == rhs.water_throughput &&

opm/simulators/wells/PerfData.hpp

@@ -57,6 +57,7 @@ public:
         serializer(micp_rates);
         serializer(cell_index);
         serializer(connection_transmissibility_factor);
+        serializer(connection_d_factor);
         serializer(satnum_id);
         serializer(ecl_index);
         serializer(water_throughput);
@@ -78,6 +79,7 @@ public:
     std::vector<double> micp_rates;
     std::vector<std::size_t> cell_index;
     std::vector<double> connection_transmissibility_factor;
+    std::vector<double> connection_d_factor;
     std::vector<int> satnum_id;
     std::vector<std::size_t> ecl_index;
 

opm/simulators/wells/PerforationData.hpp

@@ -30,6 +30,7 @@ struct PerforationData
 {
     int cell_index;
     double connection_transmissibility_factor;
+    double connection_d_factor;
     int satnum_id;
     /// \brief The original index of the perforation in ECL Schedule
     std::size_t ecl_index;

opm/simulators/wells/SingleWellState.cpp

@@ -51,6 +51,7 @@ SingleWellState::SingleWellState(const std::string& name_,
     for (std::size_t perf = 0; perf < perf_input.size(); perf++) {
         this->perf_data.cell_index[perf] = perf_input[perf].cell_index;
         this->perf_data.connection_transmissibility_factor[perf] = perf_input[perf].connection_transmissibility_factor;
+        this->perf_data.connection_d_factor[perf] = perf_input[perf].connection_d_factor;
         this->perf_data.satnum_id[perf] = perf_input[perf].satnum_id;
         this->perf_data.ecl_index[perf] = perf_input[perf].ecl_index;
     }

opm/simulators/wells/StandardWell_impl.hpp

@@ -369,9 +369,11 @@ namespace Opm
         auto& ws = well_state.well(this->index_of_well_);
         ws.phase_mixing_rates.fill(0.0);
 
+
         const int np = this->number_of_phases_;
 
         std::vector<RateVector> connectionRates = this->connectionRates_; // Copy to get right size.
+
         auto& perf_data = ws.perf_data;
         auto& perf_rates = perf_data.phase_rates;
         for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
@@ -493,7 +495,8 @@ namespace Opm
 
         PerforationRates perf_rates;
         double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants,  cell_idx);
-        const double Tw = this->well_index_[perf] * trans_mult;
+        const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+        const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
         computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
                         cq_s, perf_rates, deferred_logger);
 
@@ -1362,7 +1365,8 @@ namespace Opm
             std::vector<Scalar> mob(this->num_components_, 0.);
             getMobility(ebosSimulator, perf, mob, deferred_logger);
             double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
-            const double Tw = this->well_index_[perf] * trans_mult;
+            const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+            const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
 
             std::vector<Scalar> cq_s(this->num_components_, 0.);
             PerforationRates perf_rates;
@@ -2269,7 +2273,8 @@ namespace Opm
             getMobility(ebosSimulator, perf, mob, deferred_logger);
             std::vector<Scalar> cq_s(this->num_components_, 0.);
             double trans_mult = ebosSimulator.problem().template rockCompTransMultiplier<double>(intQuants,  cell_idx);
-            const double Tw = this->well_index_[perf] * trans_mult;
+            const auto& wellstate_nupcol = ebosSimulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+            const double Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
             PerforationRates perf_rates;
             computePerfRate(intQuants, mob, bhp.value(), Tw, perf, allow_cf,
                             cq_s, perf_rates, deferred_logger);

opm/simulators/wells/WellInterface.hpp

@@ -349,6 +349,13 @@ public:
         return 0;
     }
 
+    double wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const;
+
+    void updateConnectionDFactor(const Simulator& simulator, SingleWellState& ws) const;
+
+    void updateConnectionTransmissibilityFactor(const Simulator& simulator, SingleWellState& ws) const;
+
+
 protected:
     // simulation parameters
     const ModelParameters& param_;
@@ -439,6 +446,8 @@ protected:
                                 const std::vector<Scalar>& mobility,
                                 double* connII,
                                 DeferredLogger& deferred_logger) const;
+
+
 };
 
 }

opm/simulators/wells/WellInterface_impl.hpp

@@ -22,6 +22,7 @@
 #include <opm/common/Exceptions.hpp>
 
 #include <opm/input/eclipse/Schedule/ScheduleTypes.hpp>
+#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/simulators/wells/GroupState.hpp>
 #include <opm/simulators/wells/TargetCalculator.hpp>
@@ -1317,6 +1318,88 @@ namespace Opm
         }
     }
 
+    template <typename TypeTag>
+    double
+    WellInterface<TypeTag>::
+    wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const {
+
+        const auto& wdfac = this->well_ecl_.getWDFAC();
+        if (!wdfac.useDFactor()) {
+            return this->well_index_[perf] * trans_mult;
+        }
+
+        if constexpr (! Indices::gasEnabled) {
+            return this->well_index_[perf] * trans_mult;
+        }
+
+        // closed connection are still closed
+        if (this->well_index_[perf] == 0)
+            return 0.0;
+
+        // for gas wells we may want to add a Forchheimer term if the WDFAC or WDFACCOR keyword is used
+        const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
+        // viscosity is evaluated at connection pressure
+        const double connection_pressure = ws.perf_data.pressure[perf];
+        const double mu = FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, getValue(intQuants.fluidState().Rv()), getValue(intQuants.fluidState().Rvw()));
+        const double phi = getValue(intQuants.porosity());
+        //double k = connection.Kh()/h * trans_mult;
+        double Kh = connection.Kh()* trans_mult;
+        double Ke = connection.Ke()* trans_mult;
+        double h = Kh / Ke;
+        double rw = connection.rw();
+        double rho = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, this->pvtRegionIdx());
+        double scaling = 3.141592653589 * Kh;
+        double d = wdfac.useConnectionDFactor()? connection.dFactor() : wdfac.getDFactor(rho, mu, Ke, phi, rw, h);
+        const PhaseUsage& pu = this->phaseUsage();
+        double Q = std::abs(ws.perf_data.phase_rates[perf*pu.num_phases + pu.phase_pos[Gas]]);
+        return 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (Q/2 * d / scaling));
+    }
+
+    template <typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    updateConnectionDFactor(const Simulator& simulator, SingleWellState& ws) const {
+
+        const auto& wdfac = this->well_ecl_.getWDFAC();
+        if (!wdfac.useDFactor()) {
+            return;
+        }
+        double rho = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, this->pvtRegionIdx());
+        auto& perf_data = ws.perf_data;
+
+        for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
+            const int cell_idx = this->well_cells_[perf];
+            const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
+            const double trans_mult = simulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
+            // viscosity is evaluated at connection pressure
+            const double connection_pressure = ws.perf_data.pressure[perf];
+            const double mu = FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, getValue(intQuants.fluidState().Rv()), getValue(intQuants.fluidState().Rvw()));
+            const double phi = getValue(intQuants.porosity());
+            const auto& connection = this->well_ecl_.getConnections()[perf_data.ecl_index[perf]];
+            double Kh = connection.Kh()* trans_mult;
+            double Ke = connection.Ke()* trans_mult;
+            double h = Kh / Ke;
+            double rw = connection.rw();
+            double d = wdfac.useConnectionDFactor()? connection.dFactor() : wdfac.getDFactor(rho, mu, Ke, phi, rw, h);
+            perf_data.connection_d_factor[perf] = d;
+        }
+    }
+
+    template <typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    updateConnectionTransmissibilityFactor(const Simulator& simulator, SingleWellState& ws) const {
+        auto& perf_data = ws.perf_data;
+        for (int perf = 0; perf < this->number_of_perforations_; ++perf) {
+            const int cell_idx = this->well_cells_[perf];
+            const auto& intQuants = simulator.model().intensiveQuantities(cell_idx, /*timeIdx=*/ 0);
+            const double trans_mult = simulator.problem().template rockCompTransMultiplier<double>(intQuants, cell_idx);
+            const auto& connection = this->well_ecl_.getConnections()[perf_data.ecl_index[perf]];
+            perf_data.connection_transmissibility_factor[perf] = connection.CF() * trans_mult;
+        }
+    }
+
+
     template<typename TypeTag>
     typename WellInterface<TypeTag>::Eval
     WellInterface<TypeTag>::getPerfCellPressure(const typename WellInterface<TypeTag>::FluidState& fs) const

opm/simulators/wells/WellState.cpp

@@ -598,6 +598,7 @@ void WellState::reportConnections(std::vector<data::Connection>& connections,
         connection.pressure = perf_pressure[i];
         connection.reservoir_rate = perf_rates[i];
         connection.trans_factor = perf_data.connection_transmissibility_factor[i];
+        connection.d_factor = perf_data.connection_d_factor[i];
         if (!ws.producer) {
             const auto& filtrate_data = perf_data.filtrate_data;
             auto& filtrate = connection.filtrate;

tests/test_ParallelSerialization.cpp

@@ -111,6 +111,7 @@
 #include <opm/input/eclipse/Schedule/Well/WListManager.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPDP.hpp>
 #include <opm/input/eclipse/Schedule/Well/WVFPEXP.hpp>
+#include <opm/input/eclipse/Schedule/Well/WDFAC.hpp>
 #include <opm/input/eclipse/Schedule/WriteRestartFileEvents.hpp>
 #include <opm/input/eclipse/EclipseState/SimulationConfig/BCConfig.hpp>
 #include <opm/input/eclipse/EclipseState/SimulationConfig/RockConfig.hpp>

tests/test_parallel_wbp_calculation.cpp

@@ -420,7 +420,7 @@ namespace {
                                // 0.03, 0.0, 0.01, 0.02, 0.03, ...
                                ((k + 3 - topConn) % 4) / 100.0,
 
-                               1.0, 1.0, 0.5, 0.5, 1.0, 0.0, 0,
+                               1.0, 1.0, 0.5, 0.5, 1.0, 0.0, 0, 0.0, 0.0,
                                Opm::Connection::Direction::Z,
                                Opm::Connection::CTFKind::DeckValue, k - topConn, false);
         }

tests/test_wellstate.cpp

@@ -109,6 +109,7 @@ struct Setup
                         Opm::PerforationData pd;
                         pd.cell_index = active_index;
                         pd.connection_transmissibility_factor = completion.CF();
+                        pd.connection_d_factor = completion.dFactor();
                         pd.satnum_id = completion.satTableId();
                         well_perf_data[well_index].push_back(pd);
                     }
@@ -580,7 +581,7 @@ BOOST_AUTO_TEST_CASE(TESTPerfData) {
 
 BOOST_AUTO_TEST_CASE(TestSingleWellState) {
     Opm::ParallelWellInfo pinfo;
-    std::vector<Opm::PerforationData> connections = {{0,1,1,0},{1,1,1,1},{2,1,1,2}};
+    std::vector<Opm::PerforationData> connections = {{0,1,1,0,0},{1,1,1,0,1},{2,1,1,0,2}};
     Opm::PhaseUsage pu;
 
     // This is totally bonkers, but the pu needs a complete deck to initialize properly