move relaxationFactorFractionsForProducer to StandardWellPrimaryVariables

opm/simulators/wells/StandardWell.hpp

@@ -122,7 +122,6 @@ namespace Opm
         using Eval = typename StdWellEval::Eval;
         using EvalWell = typename StdWellEval::EvalWell;
         using BVectorWell = typename StdWellEval::BVectorWell;
-        using DiagMatrixBlockWellType = typename StdWellEval::DiagMatrixBlockWellType;
 
         StandardWell(const Well& well,
                      const ParallelWellInfo& pw_info,

opm/simulators/wells/StandardWellEval.cpp

@@ -69,53 +69,6 @@ extendEval(const Eval& in) const
     return out;
 }
 
-template<class FluidSystem, class Indices, class Scalar>
-double
-StandardWellEval<FluidSystem,Indices,Scalar>::
-relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
-                                  const BVectorWell& dwells)
-{
-    // TODO: not considering solvent yet
-    // 0.95 is a experimental value, which remains to be optimized
-    double relaxation_factor = 1.0;
-
-    if (FluidSystem::numActivePhases() > 1) {
-        if constexpr (has_wfrac_variable) {
-            const double relaxation_factor_w = StandardWellGeneric<Scalar>::
-                                               relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
-            relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
-        }
-
-        if constexpr (has_gfrac_variable) {
-            const double relaxation_factor_g = StandardWellGeneric<Scalar>::
-                                               relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
-            relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
-        }
-
-
-        if constexpr (has_wfrac_variable && has_gfrac_variable) {
-            // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will
-            // not be negative oil fraction later
-            const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
-            const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
-            const double possible_updated_sum = original_sum - relaxed_update;
-            // We only relax if fraction is above 1.
-            // The newton solver should handle the rest
-            const double epsilon = 0.001;
-            if (possible_updated_sum > 1.0 + epsilon) {
-                // since the orignal sum <= 1.0 the epsilon asserts that
-                // the relaxed_update is non trivial.
-                assert(relaxed_update != 0.);
-
-                const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
-                relaxation_factor *= further_relaxation_factor;
-            }
-        }
-        assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
-    }
-    return relaxation_factor;
-}
-
 template<class FluidSystem, class Indices, class Scalar>
 void
 StandardWellEval<FluidSystem,Indices,Scalar>::
@@ -166,7 +119,7 @@ updatePrimaryVariablesNewton(const BVectorWell& dwells,
     // 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()) ? relaxationFactorFractionsProducer(old_primary_variables, dwells)
+        (baseif_.isProducer()) ? this->primary_variables_.relaxationFactorFractionsProducer(old_primary_variables, dwells)
                                : 1.0;
 
     // update the second and third well variable (The flux fractions)

opm/simulators/wells/StandardWellEval.hpp

@@ -62,7 +62,7 @@ protected:
 public:
     using EvalWell = typename PrimaryVariables::EvalWell;
     using Eval = DenseAd::Evaluation<Scalar, Indices::numEq>;
-    using BVectorWell = typename StandardWellGeneric<Scalar>::BVectorWell;
+    using BVectorWell = typename StandardWellEquations<Scalar,Indices::numEq>::BVectorWell;
 
     //! \brief Returns a const reference to equation system.
     const StandardWellEquations<Scalar,Indices::numEq>& linSys() const
@@ -85,11 +85,6 @@ protected:
 
     EvalWell extendEval(const Eval& in) const;
 
-    // calculate a relaxation factor to avoid overshoot of the fractions for producers
-    // which might result in negative rates
-    static double relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
-                                                    const BVectorWell& dwells);
-
     // computing the accumulation term for later use in well mass equations
     void computeAccumWell();
 

opm/simulators/wells/StandardWellGeneric.hpp

@@ -48,25 +48,7 @@ class WellState;
 template<class Scalar>
 class StandardWellGeneric
 {
-protected:
+public:
-    // sparsity pattern for the matrices
-    //[A C^T    [x       =  [ res
-    // B  D ]   x_well]      res_well]
-
-    // the vector type for the res_well and x_well
-    using VectorBlockWellType = Dune::DynamicVector<Scalar>;
-    using BVectorWell = Dune::BlockVector<VectorBlockWellType>;
-
-    // the matrix type for the diagonal matrix D
-    using DiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
-    using DiagMatWell = Dune::BCRSMatrix<DiagMatrixBlockWellType>;
-
-    // the matrix type for the non-diagonal matrix B and C^T
-    using OffDiagMatrixBlockWellType = Dune::DynamicMatrix<Scalar>;
-    using OffDiagMatWell = Dune::BCRSMatrix<OffDiagMatrixBlockWellType>;
-
-    StandardWellGeneric(const WellInterfaceGeneric& baseif);
-
     // calculate a relaxation factor to avoid overshoot of total rates
     static double relaxationFactorRate(const std::vector<double>& primary_variables,
                                        const double newton_update);
@@ -75,6 +57,9 @@ protected:
     static double relaxationFactorFraction(const double old_value,
                                            const double dx);
 
+protected:
+    StandardWellGeneric(const WellInterfaceGeneric& baseif);
+
     void computeConnectionPressureDelta();
 
     // Base interface reference

opm/simulators/wells/StandardWellPrimaryVariables.cpp

@@ -34,9 +34,12 @@
 
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 
+#include <opm/simulators/wells/StandardWellGeneric.hpp>
 #include <opm/simulators/wells/WellInterfaceIndices.hpp>
 #include <opm/simulators/wells/WellState.hpp>
 
+#include <algorithm>
+
 namespace Opm {
 
 template<class FluidSystem, class Indices, class Scalar>
@@ -524,6 +527,52 @@ processFractions()
     }
 }
 
+template<class FluidSystem, class Indices, class Scalar>
+double StandardWellPrimaryVariables<FluidSystem,Indices,Scalar>::
+relaxationFactorFractionsProducer(const std::vector<double>& primary_variables,
+                                  const BVectorWell& dwells) const
+{
+    // TODO: not considering solvent yet
+    // 0.95 is a experimental value, which remains to be optimized
+    double relaxation_factor = 1.0;
+
+    if (FluidSystem::numActivePhases() > 1) {
+        if constexpr (has_wfrac_variable) {
+            const double relaxation_factor_w = StandardWellGeneric<Scalar>::
+                                               relaxationFactorFraction(primary_variables[WFrac], dwells[0][WFrac]);
+            relaxation_factor = std::min(relaxation_factor, relaxation_factor_w);
+        }
+
+        if constexpr (has_gfrac_variable) {
+            const double relaxation_factor_g = StandardWellGeneric<Scalar>::
+                                               relaxationFactorFraction(primary_variables[GFrac], dwells[0][GFrac]);
+            relaxation_factor = std::min(relaxation_factor, relaxation_factor_g);
+        }
+
+
+        if constexpr (has_wfrac_variable && has_gfrac_variable) {
+            // We need to make sure the even with the relaxation_factor, the sum of F_w and F_g is below one, so there will
+            // not be negative oil fraction later
+            const double original_sum = primary_variables[WFrac] + primary_variables[GFrac];
+            const double relaxed_update = (dwells[0][WFrac] + dwells[0][GFrac]) * relaxation_factor;
+            const double possible_updated_sum = original_sum - relaxed_update;
+            // We only relax if fraction is above 1.
+            // The newton solver should handle the rest
+            const double epsilon = 0.001;
+            if (possible_updated_sum > 1.0 + epsilon) {
+                // since the orignal sum <= 1.0 the epsilon asserts that
+                // the relaxed_update is non trivial.
+                assert(relaxed_update != 0.);
+
+                const double further_relaxation_factor = std::abs((1. - original_sum) / relaxed_update) * 0.95;
+                relaxation_factor *= further_relaxation_factor;
+            }
+        }
+        assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0);
+    }
+    return relaxation_factor;
+}
+
 #define INSTANCE(...) \
 template class StandardWellPrimaryVariables<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>;
 

opm/simulators/wells/StandardWellPrimaryVariables.hpp

@@ -132,6 +132,11 @@ public:
     //! \brief Handle non-reasonable fractions due to numerical overshoot.
     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;