Well primaryvariable scaling

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -156,6 +156,7 @@ template<class TypeTag, class MyTypeTag>
 struct MaxInnerIterWells {
     using type = UndefinedProperty;
 };
+
 template<class TypeTag, class MyTypeTag>
 struct AlternativeWellRateInit {
     using type = UndefinedProperty;
@@ -165,6 +166,24 @@ struct MaximumNumberOfWellSwitches {
     using type = UndefinedProperty;
 };
 
+
+template<class TypeTag, class MyTypeTag>
+struct WellBhpScaling {
+  using type = double;
+  static constexpr type value = 1.0;
+};
+template<class TypeTag, class MyTypeTag>
+struct WellRateScaling {
+  using type = double;
+  static constexpr type value = 1.0;
+};
+
+template<class TypeTag, class MyTypeTag>
+struct BhpControlScaling {
+  using type = double;
+  static constexpr type value = 1.0;
+};    
+
 template<class TypeTag>
 struct DbhpMaxRel<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
@@ -495,6 +514,9 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheck, "Enable the well operability checking");
             EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations");
             EWOMS_REGISTER_PARAM(TypeTag, int, MaximumNumberOfWellSwitches, "Maximum number of times a well can switch to the same control");
+	    EWOMS_REGISTER_PARAM(TypeTag, double, WellBhpScaling, "Scaling of the Bhp variable");
+	    EWOMS_REGISTER_PARAM(TypeTag, double, WellRateScaling, "Scaling of the Rate variable");
+            EWOMS_REGISTER_PARAM(TypeTag, double, BhpControlScaling, "Scaling of bhp control equation");
         }
     };
 } // namespace Opm

opm/simulators/wells/BlackoilWellModel.hpp

@@ -350,7 +350,9 @@ namespace Opm {
             double gravity_{};
             std::vector<double> depth_{};
             bool alternative_well_rate_init_{};
-
+	    double bhp_scaling_;
+	    double rate_scaling_;
+            double bhp_control_scaling_;
             std::unique_ptr<RateConverterType> rateConverter_{};
             std::unique_ptr<AverageRegionalPressureType> regionalAveragePressureCalculator_{};
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -84,7 +84,11 @@ namespace Opm {
         extractLegacyDepth_();
 
         gravity_ = ebosSimulator_.problem().gravity()[2];
-
+	// this parameters is registered in BlackoilModelParamtersEbos.hpp
+	// but no need of putting them in to the parameters.
+	bhp_scaling_ = EWOMS_GET_PARAM(TypeTag, double, WellBhpScaling);
+	rate_scaling_ = EWOMS_GET_PARAM(TypeTag, double, WellRateScaling);
+        bhp_control_scaling_ = EWOMS_GET_PARAM(TypeTag, double, BhpControlScaling);  
         initial_step_ = true;
 
         // add the eWoms auxiliary module for the wells to the list
@@ -106,6 +110,7 @@ namespace Opm {
             const bool well_opened_this_step = report_step_starts_ && events.hasEvent(wellPtr->name(), effective_events_mask);
             wellPtr->init(&this->phase_usage_, this->depth_, this->gravity_,
                           this->local_num_cells_, this->B_avg_, well_opened_this_step);
+	    wellPtr->setScalings(bhp_scaling_,rate_scaling_,bhp_control_scaling_);
         }
     }
 

opm/simulators/wells/MultisegmentWell.hpp

@@ -166,7 +166,9 @@ namespace Opm
             const SummaryState& summary_state,
             DeferredLogger& deferred_logger,
             double alq_value) const override;
-
+
+        void setScalings(double bhp_scaling, double rate_scaling, double bhp_control_scaling);
+        double bhpControlScaling() const{return this->bhp_control_scaling_;}
     protected:
         int number_segments_;
 

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -214,7 +214,7 @@ checkConvergenceControlEq(const WellState& well_state,
         }
     }
 
-    const double well_control_residual = std::abs(resWell_[0][SPres]);
+    const double well_control_residual = std::abs(resWell_[0][SPres])/this->bhp_control_scaling_;
     const int dummy_component = -1;
     if (std::isnan(well_control_residual)) {
         report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()});
@@ -411,8 +411,8 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
         // update the segment pressure
         {
             const int sign = dwells[seg][SPres] > 0.? 1 : -1;
-            const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change);
-            primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5);
+            const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change/this->bhp_scaling_);
+            primary_variables_[seg][SPres] = std::max( old_primary_variables[seg][SPres] - dx_limited, 1e5/this->bhp_scaling_);
         }
 
         // update the total rate // TODO: should we have a limitation of the total rate change?
@@ -456,7 +456,7 @@ updatePrimaryVariables(const WellState& well_state) const
         // calculate the total rate for each segment
         double total_seg_rate = 0.0;
         // the segment pressure
-        primary_variables_[seg][SPres] = segment_pressure[seg];
+        primary_variables_[seg][SPres] = segment_pressure[seg]/this->bhp_scaling_;
         // TODO: under what kind of circustances, the following will be wrong?
         // the definition of g makes the gas phase is always the last phase
         for (int p = 0; p < baseif_.numPhases(); p++) {
@@ -470,7 +470,7 @@ updatePrimaryVariables(const WellState& well_state) const
                 total_seg_rate = std::min(total_seg_rate, 0.);
             }
         }
-        primary_variables_[seg][WQTotal] = total_seg_rate;
+        primary_variables_[seg][WQTotal] = total_seg_rate/this->rate_scaling_;
         if (std::abs(total_seg_rate) > 0.) {
             if (has_wfrac_variable) {
                 const int water_pos = pu.phase_pos[Water];
@@ -607,18 +607,18 @@ getSegmentRateUpwinding(const int seg,
         if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)
                 && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == comp_idx
                 && phase == InjectorType::WATER)
-            return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
+            return primary_variables_evaluation_[seg][WQTotal]*this->rate_scaling_ / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
 
 
         if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
                 && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == comp_idx
                 && phase == InjectorType::OIL)
-            return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
+            return primary_variables_evaluation_[seg][WQTotal]*this->rate_scaling_ / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
 
         if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)
                 && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == comp_idx
                 && phase == InjectorType::GAS)
-            return primary_variables_evaluation_[seg][WQTotal] / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
+            return primary_variables_evaluation_[seg][WQTotal]*this->rate_scaling_ / baseif_.scalingFactor(baseif_.ebosCompIdxToFlowCompIdx(comp_idx));
 
         return 0.0;
     }
@@ -627,7 +627,7 @@ getSegmentRateUpwinding(const int seg,
 
     assert(segment_rate.derivative(SPres + Indices::numEq) == 0.);
 
-    return segment_rate;
+    return segment_rate*this->rate_scaling_;
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -821,7 +821,7 @@ typename MultisegmentWellEval<FluidSystem,Indices,Scalar>::EvalWell
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 getSegmentPressure(const int seg) const
 {
-    return primary_variables_evaluation_[seg][SPres];
+    return primary_variables_evaluation_[seg][SPres]*this->bhp_scaling_;
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -838,7 +838,7 @@ MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 getSegmentRate(const int seg,
                const int comp_idx) const
 {
-    return primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx);
+    return primary_variables_evaluation_[seg][WQTotal] * volumeFractionScaled(seg, comp_idx)*this->rate_scaling_;
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
@@ -1543,9 +1543,9 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
         }
 
         // update the segment pressure
-        segment_pressure[seg] = primary_variables_[seg][SPres];
+        segment_pressure[seg] = primary_variables_[seg][SPres]*this->bhp_scaling_;
         if (seg == 0) { // top segment
-            ws.bhp = segment_pressure[seg];
+            ws.bhp = segment_pressure[seg]*this->bhp_scaling_;
         }
     }
     updateThp(well_state, rho, deferred_logger);

opm/simulators/wells/MultisegmentWellGeneric.cpp

@@ -47,6 +47,9 @@ MultisegmentWellGeneric(WellInterfaceGeneric& baseif)
     , segment_inlets_(numberOfSegments())
     , segment_depth_diffs_(numberOfSegments(), 0.0)
     , perforation_segment_depth_diffs_(baseif_.numPerfs(), 0.0)
+    , bhp_scaling_(1.0)
+    , rate_scaling_(1.0)
+    , bhp_control_scaling_(1.0)  
 {
     // since we decide to use the WellSegments from the well parser. we can reuse a lot from it.
     // for other facilities needed but not available from parser, we need to process them here

opm/simulators/wells/MultisegmentWellGeneric.hpp

@@ -126,6 +126,9 @@ class MultisegmentWellGeneric
     // or in another way, the depth difference between the perforation and
     // the segment the perforation belongs to
     std::vector<double> perforation_segment_depth_diffs_;
+    double bhp_scaling_;
+    double rate_scaling_;
+    double bhp_control_scaling_;
 };
 
 }

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -2028,6 +2028,12 @@ namespace Opm
     }
 
 
+    template<typename TypeTag>
+    void MultisegmentWell<TypeTag>::setScalings(double bhp_scaling, double rate_scaling, double bhp_control_scaling){
+      this->bhp_scaling_ = bhp_scaling;
+      this->rate_scaling_ = rate_scaling;
+      this->bhp_control_scaling_ = bhp_control_scaling;
+    }
 
 
 

opm/simulators/wells/StandardWell.hpp

@@ -249,7 +249,8 @@ namespace Opm
                                     double* connII,
                                     DeferredLogger& deferred_logger) const;
 
-
+        void setScalings(double bhp_scaling, double rate_scaling, double bhp_control_scaling);
+        double bhpControlScaling() const {return this->bhp_control_scaling_;}
     protected:
         bool regularize_;
 

opm/simulators/wells/StandardWellEval.cpp

@@ -216,9 +216,9 @@ getQs(const int comp_idx) const
             //                         "Multi phase injectors are not supported, requested for well " + name());
             break;
         }
-        return inj_frac * primary_variables_evaluation_[WQTotal];
+        return inj_frac * this->getWQTotal();//primary_variables_evaluation_[WQTotal]*this->rate_scaling_;
     } else { // producers
-        return primary_variables_evaluation_[WQTotal] * wellVolumeFractionScaled(comp_idx);
+        return this->getWQTotal() * wellVolumeFractionScaled(comp_idx);
     }
 }
 
@@ -277,13 +277,13 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
     if (baseif_.isInjector()) {
         switch (baseif_.wellEcl().injectorType()) {
         case InjectorType::WATER:
-            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]];
+            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]]/this->rate_scaling_;
             break;
         case InjectorType::GAS:
-            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]];
+            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]]/this->rate_scaling_;
             break;
         case InjectorType::OIL:
-            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]];
+            primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]]/this->rate_scaling_;
             break;
         case InjectorType::MULTI:
             // Not supported.
@@ -292,7 +292,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
             break;
         }
     } else {
-            primary_variables_[WQTotal] = total_well_rate;
+            primary_variables_[WQTotal] = total_well_rate/this->rate_scaling_;
     }
 
     if (std::abs(total_well_rate) > 0.) {
@@ -353,7 +353,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
 
 
     // BHP
-    primary_variables_[Bhp] = ws.bhp;
+    primary_variables_[Bhp] = ws.bhp/this->bhp_scaling_;
 }
 
 template<class FluidSystem, class Indices, class Scalar>
@@ -665,12 +665,12 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
     }
 
     auto& ws = well_state.well(baseif_.indexOfWell());
-    ws.bhp = primary_variables_[Bhp];
+    ws.bhp = primary_variables_[Bhp]*this->bhp_scaling_;
 
     // calculate the phase rates based on the primary variables
     // for producers, this is not a problem, while not sure for injectors here
     if (baseif_.isProducer()) {
-        const double g_total = primary_variables_[WQTotal];
+        const double g_total = primary_variables_[WQTotal]*this->rate_scaling_;
         for (int p = 0; p < baseif_.numPhases(); ++p) {
             ws.surface_rates[p] = g_total * F[p];
         }
@@ -680,13 +680,13 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
         }
         switch (baseif_.wellEcl().injectorType()) {
         case InjectorType::WATER:
-            ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal];
+            ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal]*this->rate_scaling_;
             break;
         case InjectorType::GAS:
-            ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal];
+            ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal]*this->rate_scaling_;
             break;
         case InjectorType::OIL:
-            ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal];
+            ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal]*this->rate_scaling_;
             break;
         case InjectorType::MULTI:
             // Not supported.
@@ -771,9 +771,9 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
     // updating the bottom hole pressure
     {
         const int sign1 = dwells[0][Bhp] > 0 ? 1: -1;
-        const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit);
+        const double dx1_limited = sign1 * std::min(std::abs(dwells[0][Bhp]), std::abs(old_primary_variables[Bhp]) * dBHPLimit);//NB dBHPLims ser param_.dbhp_max_rel_
         // 1e5 to make sure bhp will not be below 1bar
-        primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
+        primary_variables_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5/this->bhp_scaling_);
     }
 }
 

opm/simulators/wells/StandardWellEval.hpp

@@ -108,14 +108,14 @@ class StandardWellEval : public StandardWellGeneric<Scalar>
 
     void initPrimaryVariablesEvaluation() const;
 
-    const EvalWell& getBhp() const
+    const EvalWell getBhp() const
     {
-        return primary_variables_evaluation_[Bhp];
+        return primary_variables_evaluation_[Bhp]*this->bhp_scaling_;
     }
 
-    const EvalWell& getWQTotal() const
+    const EvalWell getWQTotal() const
     {
-        return primary_variables_evaluation_[WQTotal];
+        return primary_variables_evaluation_[WQTotal]*this->rate_scaling_;
     }
 
     EvalWell extendEval(const Eval& in) const;

opm/simulators/wells/StandardWellGeneric.cpp

@@ -54,6 +54,10 @@ StandardWellGeneric(int Bhp,
     , perf_pressure_diffs_(baseif_.numPerfs())
     , parallelB_(duneB_, baseif_.parallelWellInfo())
     , Bhp_(Bhp)
+    , bhp_scaling_(1.0)
+    , rate_scaling_(1.0)
+    , bhp_control_scaling_(1.0)  
+
 {
     duneB_.setBuildMode(OffDiagMatWell::row_wise);
     duneC_.setBuildMode(OffDiagMatWell::row_wise);
@@ -68,12 +72,13 @@ relaxationFactorRate(const std::vector<double>& primary_variables,
                      const BVectorWell& dwells)
 {
     double relaxation_factor = 1.0;
-    static constexpr int WQTotal = 0;
+    static constexpr int WQTotal = 0;//NB this may be inconsiten with other definitions
 
     // For injector, we only check the total rates to avoid sign change of rates
-    const double original_total_rate = primary_variables[WQTotal];
+    // we only need scaled rates
+    const double original_total_rate = primary_variables[WQTotal];//*this->rate_scaling_;
     const double newton_update = dwells[0][WQTotal];
-    const double possible_update_total_rate = primary_variables[WQTotal] - newton_update;
+    const double possible_update_total_rate = (primary_variables[WQTotal] - newton_update);;//*this->rate_scaling_;
 
     // 0.8 here is a experimental value, which remains to be optimized
     // if the original rate is zero or possible_update_total_rate is zero, relaxation_factor will
@@ -466,7 +471,7 @@ checkConvergenceControlEq(const WellState& well_state,
         }
     }
 
-    const double well_control_residual = std::abs(this->resWell_[0][Bhp_]);
+    const double well_control_residual = std::abs(this->resWell_[0][Bhp_])/this->bhp_control_scaling_;
     const int dummy_component = -1;
     if (std::isnan(well_control_residual)) {
         report.setWellFailed({ctrltype, CR::Severity::NotANumber, dummy_component, baseif_.name()});

opm/simulators/wells/StandardWellGeneric.hpp

@@ -32,7 +32,6 @@
 
 #include <optional>
 #include <vector>
-
 namespace Opm
 {
 
@@ -137,6 +136,11 @@ class StandardWellGeneric
 
 private:
     int Bhp_; // index of Bhp
+protected:
+     double bhp_scaling_;
+     double rate_scaling_;
+     double bhp_control_scaling_;
+
 };
 
 }