From b39b9a85a2854b64f4ed34b9b17e0edab3761ff5 Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Wed, 30 Mar 2016 15:41:03 +0200
Subject: [PATCH] making the interleaved solver works for blackoil polymer
 simulator.

The CPR solver does not work yet. An error will be thrown if people
specify to use CPR linear solver.
---
 opm/autodiff/FlowMainPolymer.hpp              | 32 +++++++++++++++++--
 .../fullyimplicit/BlackoilPolymerModel.hpp    |  2 ++
 .../BlackoilPolymerModel_impl.hpp             | 21 ++++++++++++
 3 files changed, 53 insertions(+), 2 deletions(-)

diff --git a/opm/autodiff/FlowMainPolymer.hpp b/opm/autodiff/FlowMainPolymer.hpp
index b594b3c0b..1a5e647f8 100644
--- a/opm/autodiff/FlowMainPolymer.hpp
+++ b/opm/autodiff/FlowMainPolymer.hpp
@@ -38,6 +38,10 @@ namespace Opm
     {
     protected:
         using Base = FlowMainBase<FlowMainPolymer<Grid, Simulator>, Grid, Simulator>;
+        using Base::eclipse_state_;
+        using Base::param_;
+        using Base::fis_solver_;
+        using Base::parallel_information_;
         friend Base;
 
         // Set in setupGridAndProps()
@@ -90,10 +94,34 @@ namespace Opm
         // Setup linear solver.
         // Writes to:
         //   fis_solver_
+        // Currently, the CPR solver is not ready for polymer solver yet
         void setupLinearSolver()
         {
-            OPM_MESSAGE("Caution: polymer solver currently only works with direct linear solver.");
-            Base::fis_solver_.reset(new NewtonIterationBlackoilSimple(Base::param_, Base::parallel_information_));
+            const std::string cprSolver = "cpr";
+            const std::string interleavedSolver = "interleaved";
+            const std::string directSolver = "direct";
+            const std::string flowDefaultSolver = interleavedSolver;
+
+            std::shared_ptr<const Opm::SimulationConfig> simCfg = eclipse_state_->getSimulationConfig();
+            std::string solver_approach = flowDefaultSolver;
+
+            if (param_.has("solver_approach")) {
+                solver_approach = param_.template get<std::string>("solver_approach");
+            }  else {
+                if (simCfg->useCPR()) {
+                    solver_approach = cprSolver;
+                }
+            }
+
+            if (solver_approach == cprSolver) {
+                OPM_THROW( std::runtime_error , "CPR solver is not ready for use with polymer solver yet.");
+            } else if (solver_approach == interleavedSolver) {
+                fis_solver_.reset(new NewtonIterationBlackoilInterleaved(param_, parallel_information_));
+            } else if (solver_approach == directSolver) {
+                fis_solver_.reset(new NewtonIterationBlackoilSimple(param_, parallel_information_));
+            } else {
+                OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized.");
+            }
         }
 
 
diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
index 53d452e21..da06207cf 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
@@ -234,6 +234,8 @@ namespace Opm {
                                            const SolutionState& state,
                                            WellState& xw);
 
+        void updateEquationsScaling();
+
         void
         computeMassFlux(const int               actph ,
                         const V&                transi,
diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
index 83dead022..fe2a98056 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
@@ -107,6 +107,7 @@ namespace Opm {
             if (!active_[Water]) {
                 OPM_THROW(std::logic_error, "Polymer must solved in water!\n");
             }
+            residual_.matbalscale.resize(fluid_.numPhases() + 1, 1.1169); // use the same as the water phase
             // If deck has polymer, residual_ should contain polymer equation.
             rq_.resize(fluid_.numPhases() + 1);
             residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null());
@@ -339,6 +340,26 @@ namespace Opm {
             residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0])
                                                + ops_.div*rq_[poly_pos_].mflux;
         }
+
+
+        if (param_.update_equations_scaling_) {
+            updateEquationsScaling();
+        }
+
+    }
+
+
+
+
+
+    template <class Grid>
+    void BlackoilPolymerModel<Grid>::updateEquationsScaling()
+    {
+        Base::updateEquationsScaling();
+        if (has_polymer_) {
+            const int water_pos = fluid_.phaseUsage().phase_pos[Water];
+            residual_.matbalscale[poly_pos_] = residual_.matbalscale[water_pos];
+        }
     }