Merge pull request #4996 from totto82/fixWDFAC

Fix wdfac and d factor in comdat

opm/simulators/wells/MultisegmentWell.hpp

@@ -188,7 +188,7 @@ namespace Opm
         template<class Value>
         void computePerfRate(const IntensiveQuantities& int_quants,
                              const std::vector<Value>& mob_perfcells,
-                             const double Tw,
+                             const std::vector<Scalar>& Tw,
                              const int seg,
                              const int perf,
                              const Value& segment_pressure,
@@ -204,7 +204,7 @@ namespace Opm
                         const Value& rv,
                         const std::vector<Value>& b_perfcells,
                         const std::vector<Value>& mob_perfcells,
-                        const double Tw,
+                        const std::vector<Scalar>& Tw,
                         const int perf,
                         const Value& segment_pressure,
                         const Value& segment_density,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -349,7 +349,7 @@ namespace Opm
                 getMobility(ebosSimulator, perf, mob, deferred_logger);
                 const double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
                 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 std::vector<Scalar> 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;
@@ -768,7 +768,7 @@ namespace Opm
                     const Value& rv,
                     const std::vector<Value>& b_perfcells,
                     const std::vector<Value>& mob_perfcells,
-                    const double Tw,
+                    const std::vector<Scalar>& Tw,
                     const int perf,
                     const Value& segment_pressure,
                     const Value& segment_density,
@@ -804,7 +804,7 @@ namespace Opm
 
             // compute component volumetric rates at standard conditions
             for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
-                const Value cq_p = - Tw * (mob_perfcells[comp_idx] * drawdown);
+                const Value cq_p = - Tw[comp_idx] * (mob_perfcells[comp_idx] * drawdown);
                 cq_s[comp_idx] = b_perfcells[comp_idx] * cq_p;
             }
 
@@ -828,9 +828,6 @@ namespace Opm
                 total_mob += mob_perfcells[comp_idx];
             }
 
-            // injection perforations total volume rates
-            const Value cqt_i = - Tw * (total_mob * drawdown);
-
             // compute volume ratio between connection and at standard conditions
             Value volume_ratio = 0.0;
             if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
@@ -871,9 +868,10 @@ namespace Opm
                 }
             }
             // injecting connections total volumerates at standard conditions
-            Value cqt_is = cqt_i / volume_ratio;
+            for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
-            for (int comp_idx = 0; comp_idx < this->numComponents(); ++comp_idx) {
+                const Value cqt_i = - Tw[componentIdx] * (total_mob * drawdown);
-                cq_s[comp_idx] = cmix_s[comp_idx] * cqt_is;
+                Value cqt_is = cqt_i / volume_ratio;
+                cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
             }
         } // end for injection perforations
 
@@ -907,7 +905,7 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     computePerfRate(const IntensiveQuantities& int_quants,
                     const std::vector<Value>& mob_perfcells,
-                    const double Tw,
+                    const std::vector<Scalar>& Tw,
                     const int seg,
                     const int perf,
                     const Value& segment_pressure,
@@ -1182,11 +1180,11 @@ namespace Opm
                 // the well index associated with the connection
                 const double trans_mult = ebos_simulator.problem().template wellTransMultiplier<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);  
+                const std::vector<Scalar> 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) {
-                    const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
+                    const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
                     ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
                     ipr_b_perf[comp_idx] += tw_mob;
                 }
@@ -1682,7 +1680,7 @@ namespace Opm
                 getMobility(ebosSimulator, perf, mob, deferred_logger);
                 const double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(int_quants, cell_idx);
                 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);
+                const std::vector<Scalar> 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;
@@ -1995,7 +1993,7 @@ namespace Opm
                 getMobility(ebosSimulator, perf, mob, deferred_logger);
                 const double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(int_quants, cell_idx);
                 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);
+                const std::vector<Scalar> 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/StandardWell.hpp

@@ -286,7 +286,7 @@ namespace Opm
         void computePerfRate(const IntensiveQuantities& intQuants,
                              const std::vector<Value>& mob,
                              const Value& bhp,
-                             const double Tw,
+                             const std::vector<Scalar>& Tw,
                              const int perf,
                              const bool allow_cf,
                              std::vector<Value>& cq_s,
@@ -302,7 +302,7 @@ namespace Opm
                              const Value& rvw,
                              const Value& rsw,
                              std::vector<Value>& b_perfcells_dense,
-                             const double Tw,
+                             const std::vector<Scalar>& Tw,
                              const int perf,
                              const bool allow_cf,
                              const Value& skin_pressure,

opm/simulators/wells/StandardWell_impl.hpp

@@ -119,7 +119,7 @@ namespace Opm
     computePerfRate(const IntensiveQuantities& intQuants,
                     const std::vector<Value>& mob,
                     const Value& bhp,
-                    const double Tw,
+                    const std::vector<Scalar>& Tw,
                     const int perf,
                     const bool allow_cf,
                     std::vector<Value>& cq_s,
@@ -224,7 +224,7 @@ namespace Opm
                     const Value& rvw,
                     const Value& rsw,
                     std::vector<Value>& b_perfcells_dense,
-                    const double Tw,
+                    const std::vector<Scalar>& Tw,
                     const int perf,
                     const bool allow_cf,
                     const Value& skin_pressure,
@@ -249,7 +249,7 @@ namespace Opm
 
             // compute component volumetric rates at standard conditions
             for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
-                const Value cq_p = - Tw * (mob[componentIdx] * drawdown);
+                const Value cq_p = - Tw[componentIdx] * (mob[componentIdx] * drawdown);
                 cq_s[componentIdx] = b_perfcells_dense[componentIdx] * cq_p;
             }
 
@@ -270,9 +270,6 @@ namespace Opm
                 total_mob_dense += mob[componentIdx];
             }
 
-            // injection perforations total volume rates
-            const Value cqt_i = - Tw * (total_mob_dense * drawdown);
-
             // compute volume ratio between connection at standard conditions
             Value volumeRatio = bhp * 0.0; // initialize it with the correct type
 ;
@@ -306,8 +303,9 @@ namespace Opm
             }
 
             // injecting connections total volumerates at standard conditions
-            Value cqt_is = cqt_i / volumeRatio;
             for (int componentIdx = 0; componentIdx < this->numComponents(); ++componentIdx) {
+                const Value cqt_i = - Tw[componentIdx] * (total_mob_dense * drawdown);
+                Value cqt_is = cqt_i / volumeRatio;
                 cq_s[componentIdx] = cmix_s[componentIdx] * cqt_is;
             }
 
@@ -496,7 +494,7 @@ namespace Opm
         PerforationRates perf_rates;
         double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants,  cell_idx);
         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 std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
         computePerfRate(intQuants, mob, bhp, Tw, perf, allow_cf,
                         cq_s, perf_rates, deferred_logger);
 
@@ -790,12 +788,13 @@ namespace Opm
             }
 
             // the well index associated with the connection
-            const double tw_perf = this->well_index_[perf]*ebos_simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
+            double trans_mult = ebos_simulator.problem().template wellTransMultiplier<double>(int_quantities, cell_idx);
-
+            const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+            const std::vector<Scalar> 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) {
-                const double tw_mob = tw_perf * mob[comp_idx] * b_perf[comp_idx];
+                const double tw_mob = tw_perf[comp_idx] * mob[comp_idx] * b_perf[comp_idx];
                 ipr_a_perf[comp_idx] += tw_mob * pressure_diff;
                 ipr_b_perf[comp_idx] += tw_mob;
             }
@@ -1366,7 +1365,7 @@ namespace Opm
             getMobility(ebosSimulator, perf, mob, deferred_logger);
             double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants, cell_idx);
             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 std::vector<Scalar> Tw = this->wellIndex(perf, intQuants, trans_mult, wellstate_nupcol);
 
             std::vector<Scalar> cq_s(this->num_components_, 0.);
             PerforationRates perf_rates;
@@ -1669,7 +1668,8 @@ namespace Opm
             std::vector<EvalWell> cq_s(this->num_components_, {this->primary_variables_.numWellEq() + Indices::numEq, 0.});
             PerforationRates perf_rates;
             double trans_mult = ebos_simulator.problem().template wellTransMultiplier<double>(int_quant, cell_idx);
-            const double Tw = this->well_index_[perf] * trans_mult;
+            const auto& wellstate_nupcol = ebos_simulator.problem().wellModel().nupcolWellState().well(this->index_of_well_);
+            const std::vector<Scalar> Tw = this->wellIndex(perf, int_quant, trans_mult, wellstate_nupcol);  
             computePerfRate(int_quant, mob, bhp, Tw, perf, allow_cf, cq_s,
                             perf_rates, deferred_logger);
             // TODO: make area a member
@@ -2274,7 +2274,7 @@ namespace Opm
             std::vector<Scalar> cq_s(this->num_components_, 0.);
             double trans_mult = ebosSimulator.problem().template wellTransMultiplier<double>(intQuants,  cell_idx);
             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 std::vector<Scalar> 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,7 +349,7 @@ public:
         return 0;
     }
 
-    double wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const;
+    std::vector<double> wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const;
 
     void updateConnectionDFactor(const Simulator& simulator, SingleWellState& ws) const;
 
@@ -447,6 +447,8 @@ protected:
                                 double* connII,
                                 DeferredLogger& deferred_logger) const;
 
+    double computeConnectionDFactor(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const double total_tw, const SingleWellState& ws) const;
+
 
 };
 

opm/simulators/wells/WellInterface_impl.hpp

@@ -1324,40 +1324,37 @@ namespace Opm
     }
 
     template <typename TypeTag>
-    double
+    std::vector<double>
     WellInterface<TypeTag>::
     wellIndex(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const SingleWellState& ws) const {
 
+        std::vector<Scalar> wi(this->num_components_, this->well_index_[perf] * trans_mult);
         const auto& wdfac = this->well_ecl_.getWDFAC();
         if (!wdfac.useDFactor()) {
-            return this->well_index_[perf] * trans_mult;
+            return wi;
         }
-
+        // for gas wells we may want to add a Forchheimer term if the WDFAC or WDFACCOR keyword is used        
         if constexpr (! Indices::gasEnabled) {
-            return this->well_index_[perf] * trans_mult;
+            return wi;
         }
-
         // closed connection are still closed
         if (this->well_index_[perf] == 0)
-            return 0.0;
+            return std::vector<Scalar>(this->num_components_, 0.0);
 
-        // for gas wells we may want to add a Forchheimer term if the WDFAC or WDFACCOR keyword is used
+        double tot_tw = 0.0;
-        const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
+        for (const auto& c : this->well_ecl_.getConnections()) {
-        // viscosity is evaluated at connection pressure
+            tot_tw += c.CF();
-        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 d = computeConnectionDFactor(perf, intQuants, trans_mult, tot_tw, ws);
-        //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));
+        const auto& connection = this->well_ecl_.getConnections()[ws.perf_data.ecl_index[perf]];
+        double Kh = connection.Kh()* trans_mult;
+        double scaling = 3.141592653589 * Kh;
+        const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+        wi[gas_comp_idx]  = 1.0/(1.0/(trans_mult * this->well_index_[perf]) + (Q/2 * d / scaling));
+        return wi;
     }
 
     template <typename TypeTag>
@@ -1369,27 +1366,54 @@ namespace Opm
         if (!wdfac.useDFactor()) {
             return;
         }
-        double rho = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, this->pvtRegionIdx());
         auto& perf_data = ws.perf_data;
-
+        double tot_tw = 0.0;
+        for (const auto& c : this->well_ecl_.getConnections()) {
+            tot_tw += c.CF();
+        }
         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
+            perf_data.connection_d_factor[perf] = computeConnectionDFactor(perf, intQuants, trans_mult, tot_tw, ws);
-            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>
+    double 
+    WellInterface<TypeTag>::
+    computeConnectionDFactor(const int perf, const IntensiveQuantities& intQuants, const double trans_mult, const double total_tw, const SingleWellState& ws) const {
+        const double connection_pressure = ws.perf_data.pressure[perf];
+        // viscosity is evaluated at connection pressure
+        const auto& rv = getValue(intQuants.fluidState().Rv());
+        const double psat = FluidSystem::gasPvt().saturationPressure(this->pvtRegionIdx(), ws.temperature, rv);
+        const double mu = connection_pressure < psat ?
+                                FluidSystem::gasPvt().saturatedViscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure) :
+                                FluidSystem::gasPvt().viscosity(this->pvtRegionIdx(), ws.temperature, connection_pressure, rv, getValue(intQuants.fluidState().Rvw()));
+        double rho = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, this->pvtRegionIdx());
+        const double phi = getValue(intQuants.porosity());
+        const auto& connection = this->well_ecl_.getConnections()[ws.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();
+        const auto& wdfac = this->well_ecl_.getWDFAC();
+        if (wdfac.useConnectionDFactor()) {
+            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 well index
+            return d * total_tw / connection.CF();
+        } else {
+
+            double d = wdfac.getDFactor(rho, mu, Ke, phi, rw, h);
+            return d * total_tw / connection.CF();
+        }
+    }
+
+
     template <typename TypeTag>
     void
     WellInterface<TypeTag>::