From d1727d018346a934ad4e5d8c167acf16777295ac Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Mon, 31 Jul 2017 15:04:25 +0200
Subject: [PATCH] correcting the injectivity with polymer injection.

---
 opm/autodiff/StandardWell.hpp      |  1 +
 opm/autodiff/StandardWell_impl.hpp | 44 ++++++++++++++++++++++++++++++
 opm/autodiff/WellInterface.hpp     |  2 ++
 3 files changed, 47 insertions(+)

diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp
index a8797df79..ca7868ea0 100644
--- a/opm/autodiff/StandardWell.hpp
+++ b/opm/autodiff/StandardWell.hpp
@@ -62,6 +62,7 @@ namespace Opm
         using typename WellInterface<TypeTag>::Mat;
         using typename WellInterface<TypeTag>::BVector;
         using typename WellInterface<TypeTag>::Eval;
+        using typename WellInterface<TypeTag>::PolymerModule;
 
         using WellInterface<TypeTag>::numEq;
         static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? numEq-1 : numEq; // //numEq; //number of wellEq is only numEq for polymer
diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp
index 13829e2a6..88b816601 100644
--- a/opm/autodiff/StandardWell_impl.hpp
+++ b/opm/autodiff/StandardWell_impl.hpp
@@ -797,6 +797,50 @@ namespace Opm
                 OPM_THROW(std::runtime_error, "individual mobility for wells does not work in combination with solvent");
             }
         }
+
+        // modify the water mobility if polymer is present
+        if (has_polymer) {
+            // assume fully mixture for wells.
+            EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
+
+            if (wellType() == INJECTOR) {
+                const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(intQuants.pvtRegionIndex());
+                mob[ Water ] /= (extendEval(intQuants.waterViscosityCorrection()) * viscosityMultiplier.eval(polymerConcentration, /*extrapolate=*/true) );
+            }
+
+            /* if (PolymerModule::hasPlyshlog()) {
+                // compute the well water velocity with out shear effects.
+                const int numComp = numComponents();
+                const bool allow_cf = crossFlowAllowed(ebosSimulator);
+                const EvalWell& bhp = getBhp();
+                std::vector<EvalWell> cq_s(numComp,0.0);
+                computePerfRate(intQuants, mob, wellIndex()[perf], bhp, perfPressureDiffs()[perf], allow_cf, cq_s);
+                double area = 2 * M_PI * wells_rep_radius_[perf] * wells_perf_length_[perf];
+                const auto& materialLawManager = ebosSimulator.problem().materialLawManager();
+                const auto& scaledDrainageInfo =
+                        materialLawManager->oilWaterScaledEpsInfoDrainage(cell_idx);
+                const Scalar& Swcr = scaledDrainageInfo.Swcr;
+                const EvalWell poro = extendEval(intQuants.porosity());
+                const EvalWell Sw = extendEval(intQuants.fluidState().saturation(flowPhaseToEbosPhaseIdx(Water)));
+                // guard against zero porosity and no water
+                const EvalWell denom = Opm::max( (area * poro * (Sw - Swcr)), 1e-12);
+                EvalWell waterVelocity = cq_s[ Water ] / denom * extendEval(intQuants.fluidState().invB(flowPhaseToEbosPhaseIdx(Water)));
+
+                if (PolymerModule::hasShrate()) {
+                    // TODO Use the same conversion as for the reservoar equations.
+                    // Need the "permeability" of the well?
+                    // For now use the same formula as in legacy.
+                    waterVelocity *= PolymerModule::shrate( intQuants.pvtRegionIndex() ) / wells_bore_diameter_[perf];
+                }
+                EvalWell polymerConcentration = extendEval(intQuants.polymerConcentration());
+                EvalWell shearFactor = PolymerModule::computeShearFactor(polymerConcentration,
+                                                                         intQuants.pvtRegionIndex(),
+                                                                         waterVelocity);
+
+                // modify the mobility with the shear factor and recompute the well fluxes.
+                mob[ Water ] /= shearFactor;
+            } */
+        }
     }
 
 
diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp
index d494363a6..1287ebcae 100644
--- a/opm/autodiff/WellInterface.hpp
+++ b/opm/autodiff/WellInterface.hpp
@@ -78,6 +78,8 @@ namespace Opm
         typedef Dune::BlockVector<VectorBlockType> BVector;
         typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
 
+        typedef Ewoms::BlackOilPolymerModule<TypeTag> PolymerModule;
+
         static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
         static const bool has_polymer = GET_PROP_VALUE(TypeTag, EnablePolymer);