using effective polymer viscosity for polymer mobility calculation.

opm/polymer/PolymerProperties.cpp

@@ -243,6 +243,52 @@ namespace Opm
         }
     }
 
+    void PolymerProperties::effectiveInvPolyVisc(const double c,
+                                                 const double* visc,
+                                                 double& inv_mu_p_eff) const
+    {
+        double dummy;
+        effectiveInvPolyViscBoth(c, visc, inv_mu_p_eff, dummy, false);
+
+    }
+
+    void PolymerProperties::effectiveInvPolyViscWithDer(const double c,
+                                                        const double* visc,
+                                                        double& inv_mu_p_eff,
+                                                        double& d_inv_mu_p_eff_dc) const
+    {
+        effectiveInvPolyViscBoth(c, visc, inv_mu_p_eff, d_inv_mu_p_eff_dc, true);
+
+    }
+
+    void PolymerProperties::effectiveInvPolyViscBoth(const double c,
+                                                     const double* visc,
+                                                     double& inv_mu_p_eff,
+                                                     double& dinv_mu_p_eff_dc,
+                                                     const bool if_with_der) const
+    {
+        const double omega = mix_param_;
+        const double mu_w = visc[0];
+
+        double mu_m = 0.0;
+        double dmu_m_dc = 0.0;
+
+        if (if_with_der) {
+            mu_m = viscMultWithDer(c, &dmu_m_dc)*mu_w;
+            dmu_m_dc *= mu_w;
+        } else {
+            mu_m = viscMult(c)*mu_w;
+        }
+
+        const double inv_mu_m_omega = std::pow(mu_m, -omega);
+        const double mu_p = viscMult(c_max_) * mu_w;
+        inv_mu_p_eff = inv_mu_m_omega * std::pow(mu_p, omega - 1.);
+
+        if (if_with_der) {
+            dinv_mu_p_eff_dc = -omega * dmu_m_dc * std::pow(mu_m, -omega - 1) * std::pow(mu_p, omega - 1);
+        }
+    }
+
     void PolymerProperties::effectiveRelperm(const double c,
                                              const double cmax,
                                              const double* relperm,

opm/polymer/PolymerProperties.hpp

@@ -292,6 +292,15 @@ namespace Opm
                                                         const double* visc,
                                                         double& inv_mu_w_eff,
                                                         double& dinv_mu_w_eff_dc) const;
+        void effectiveInvPolyVisc(const double c,
+                                  const double* visc,
+                                  double& inv_mu_p_eff) const;
+
+        void effectiveInvPolyViscWithDer(const double c,
+                                         const double* visc,
+                                         double& inv_mu_p_eff,
+                                         double& d_inv_mu_p_eff_dc) const;
+
         void effectiveRelperm(const double c,
                               const double cmax,
                               const double* relperm,
@@ -401,6 +410,13 @@ namespace Opm
         void effectiveInvViscBoth(const double c, const double* visc,
                                   double& inv_mu_w_eff,
                                   double& dinv_mu_w_eff_dc, bool if_with_der) const;
+
+        void effectiveInvPolyViscBoth(const double c,
+                                      const double* visc,
+                                      double& inv_mu_p_eff,
+                                      double& dinv_mu_p_eff_dc,
+                                      const bool if_with_der) const;
+
         void effectiveRelpermBoth(const double c,
                                   const double cmax,
                                   const double* relperm,

opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp

@@ -464,7 +464,8 @@ namespace Opm {
                 // Reduce mobility of water phase by relperm reduction and effective viscosity increase.
                 rq_[actph].mob = tr_mult * krw_eff * inv_wat_eff_visc;
                 // Compute polymer mobility.
-                rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_wat_eff_visc;
+                const ADB inv_poly_eff_visc = polymer_props_ad_.effectiveInvPolymerVisc(state.concentration, mu.value().data());
+                rq_[poly_pos_].mob = tr_mult * mc * krw_eff * inv_poly_eff_visc;
                 rq_[poly_pos_].b = rq_[actph].b;
                 rq_[poly_pos_].dh = rq_[actph].dh;
                 UpwindSelector<double> upwind(grid_, ops_, rq_[poly_pos_].dh.value());

opm/polymer/fullyimplicit/PolymerPropsAd.cpp

@@ -192,6 +192,34 @@ namespace Opm {
 
 
 
+    ADB PolymerPropsAd::effectiveInvPolymerVisc(const ADB& c, const double* visc) const
+    {
+        const int nc = c.size();
+        V inv_mu_p_eff(nc);
+        V dinv_mu_p_eff(nc);
+        for (int i = 0; i < nc; ++i) {
+            double im = 0;
+            double dim = 0;
+            // TODO: the usage of visc can be likely wrong, while more investigation will be requried.
+            polymer_props_.effectiveInvPolyViscWithDer(c.value()(i), visc, im, dim);
+            inv_mu_p_eff(i) = im;
+            dinv_mu_p_eff(i) = dim;
+        }
+
+        ADB::M dim_diag(dinv_mu_p_eff.matrix().asDiagonal());
+        const int num_blocks = c.numBlocks();
+        std::vector<ADB::M> jacs(num_blocks);
+
+        for (int block = 0; block < num_blocks; ++block) {
+            jacs[block] = dim_diag * c.derivative()[block];
+        }
+        return ADB::function(std::move(inv_mu_p_eff), std::move(jacs));
+    }
+
+
+
+
+
     V PolymerPropsAd::polymerWaterVelocityRatio(const V& c) const
     {
         const int nc = c.size();

opm/polymer/fullyimplicit/PolymerPropsAd.hpp

@@ -92,6 +92,12 @@ namespace Opm {
         ADB 
         effectiveInvWaterVisc(const ADB& c,const double* visc) const;
 
+		/// \param[in] c		ADB of polymer concentraion values.
+		/// \param[in] visc		Array of water viscosity values
+		/// \return 			ADB of inverse effective polymer viscosity.
+        ADB
+        effectiveInvPolymerVisc(const ADB& c, const double* visc) const;
+		
 		/// \param[in] c		Array of n polymer concentraion values.
 		/// \return 			Array of n mc values, here mc means m(c) * c.
         V