move updatePrimaryVariablesNewton to StandardWellPrimaryVariables

opm/simulators/wells/StandardWellEval.cpp

@@ -107,57 +107,6 @@ computeAccumWell()
     }
 }
 
-template<class FluidSystem, class Indices, class Scalar>
-void
-StandardWellEval<FluidSystem,Indices,Scalar>::
-updatePrimaryVariablesNewton(const BVectorWell& dwells,
-                             [[maybe_unused]] const double dFLimit,
-                             const double dBHPLimit) const
-{
-    const std::vector<double> old_primary_variables = primary_variables_.value_;
-
-    // for injectors, very typical one of the fractions will be one, and it is easy to get zero value
-    // fractions. not sure what is the best way to handle it yet, so we just use 1.0 here
-    [[maybe_unused]] const double relaxation_factor_fractions =
-        (baseif_.isProducer()) ? this->primary_variables_.relaxationFactorFractionsProducer(old_primary_variables, dwells)
-                               : 1.0;
-
-    // update the second and third well variable (The flux fractions)
-
-    if constexpr (has_wfrac_variable) {
-        const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
-        const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
-        // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
-        primary_variables_.value_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
-    }
-
-    if constexpr (has_gfrac_variable) {
-        const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
-        const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
-        primary_variables_.value_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
-    }
-
-    if constexpr (Indices::enableSolvent) {
-        const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
-        const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
-        primary_variables_.value_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
-    }
-
-    this->primary_variables_.processFractions();
-
-    // updating the total rates Q_t
-    const double relaxation_factor_rate = this->relaxationFactorRate(old_primary_variables, dwells[0][WQTotal]);
-    primary_variables_.value_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
-
-    // 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);
-        // 1e5 to make sure bhp will not be below 1bar
-        primary_variables_.value_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
-    }
-}
-
 template<class FluidSystem, class Indices, class Scalar>
 ConvergenceReport
 StandardWellEval<FluidSystem,Indices,Scalar>::

opm/simulators/wells/StandardWellEval.hpp

@@ -118,10 +118,6 @@ protected:
     void updateWellStateFromPrimaryVariables(WellState& well_state,
                                              DeferredLogger& deferred_logger) const;
 
-    void updatePrimaryVariablesNewton(const BVectorWell& dwells,
-                                      const double dFLimit,
-                                      const double dBHPLimit) const;
-
     mutable PrimaryVariables primary_variables_; //!< Primary variables for well
 
     // the saturations in the well bore under surface conditions at the beginning of the time step

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -215,6 +215,55 @@ updatePolyMW(const BVectorWell& dwells)
     }
 }
 
+template<class FluidSystem, class Indices, class Scalar>
+void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
+updateNewton(const BVectorWell& dwells,
+             [[maybe_unused]] const double dFLimit,
+             const double dBHPLimit)
+{
+    const std::vector<double> old_primary_variables = value_;
+
+    // for injectors, very typical one of the fractions will be one, and it is easy to get zero value
+    // fractions. not sure what is the best way to handle it yet, so we just use 1.0 here
+    [[maybe_unused]] const double relaxation_factor_fractions =
+        well_.isProducer() ? relaxationFactorFractionsProducer(old_primary_variables, dwells)
+                           : 1.0;
+
+    // update the second and third well variable (The flux fractions)
+
+    if constexpr (has_wfrac_variable) {
+        const int sign2 = dwells[0][WFrac] > 0 ? 1: -1;
+        const double dx2_limited = sign2 * std::min(std::abs(dwells[0][WFrac] * relaxation_factor_fractions), dFLimit);
+        // primary_variables_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
+        value_[WFrac] = old_primary_variables[WFrac] - dx2_limited;
+    }
+
+    if constexpr (has_gfrac_variable) {
+        const int sign3 = dwells[0][GFrac] > 0 ? 1: -1;
+        const double dx3_limited = sign3 * std::min(std::abs(dwells[0][GFrac] * relaxation_factor_fractions), dFLimit);
+        value_[GFrac] = old_primary_variables[GFrac] - dx3_limited;
+    }
+
+    if constexpr (Indices::enableSolvent) {
+        const int sign4 = dwells[0][SFrac] > 0 ? 1: -1;
+        const double dx4_limited = sign4 * std::min(std::abs(dwells[0][SFrac]) * relaxation_factor_fractions, dFLimit);
+        value_[SFrac] = old_primary_variables[SFrac] - dx4_limited;
+    }
+
+    this->processFractions();
+
+    // updating the total rates Q_t
+    const double relaxation_factor_rate = StandardWellGeneric<Scalar>::relaxationFactorRate(old_primary_variables, dwells[0][WQTotal]);
+    value_[WQTotal] = old_primary_variables[WQTotal] - dwells[0][WQTotal] * relaxation_factor_rate;
+
+    // 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);
+    // 1e5 to make sure bhp will not be below 1bar
+    value_[Bhp] = std::max(old_primary_variables[Bhp] - dx1_limited, 1e5);
+}
+
 template<class FluidSystem, class Indices, class Scalar>
 void StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
 copyToWellState(WellState& well_state,

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -111,7 +111,12 @@ public:
     //! \brief Copy values from well state.
     void update(const WellState& well_state, DeferredLogger& deferred_logger);
 
-    //! \brief Update polymer molecular weight values from solution vector.
+    //! \brief Update values from newton update vector.
+    void updateNewton(const BVectorWell& dwells,
+                      const double dFLimit,
+                      const double dBHPLimit);
+
+    //! \brief Update polymer molecular weight values from newton update vector.
     void updatePolyMW(const BVectorWell& dwells);
 
     //! \brief Copy values to well state.
@@ -129,18 +134,18 @@ public:
     //! \brief Returns scaled rate for a component.
     EvalWell getQs(const int compIdx) const;
 
-    //! \brief Handle non-reasonable fractions due to numerical overshoot.
+private:
-    void processFractions();
-
     //! \brief Calculate a relaxation factor for producers.
     //! \details To avoid overshoot of the fractions which might result in negative rates.
     double relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
                                              const BVectorWell& dwells) const;
 
-private:
     //! \brief Returns volume fraction for a component.
     EvalWell volumeFraction(const unsigned compIdx) const;
 
+    //! \brief Handle non-reasonable fractions due to numerical overshoot.
+    void processFractions();
+
     const WellInterfaceIndices<FluidSystem,Indices,Scalar>& well_; //!< Reference to well interface
 
     //! \brief Total number of the well equations and primary variables.

opm/simulators/wells/StandardWell_impl.hpp

@@ -957,7 +957,7 @@ namespace Opm
     {
         const double dFLimit = this->param_.dwell_fraction_max_;
         const double dBHPLimit = this->param_.dbhp_max_rel_;
-        this->StdWellEval::updatePrimaryVariablesNewton(dwells, dFLimit, dBHPLimit);
+        this->primary_variables_.updateNewton(dwells, dFLimit, dBHPLimit);
 
         updateExtraPrimaryVariables(dwells);