From b1b9c171b3ea23d44b289bc4bdba5f2f585daf0f Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Thu, 27 Aug 2015 16:58:44 +0200
Subject: [PATCH] update scaling factors for mass balance equations.

---
 opm/autodiff/BlackoilModelBase.hpp            |  3 ++
 opm/autodiff/BlackoilModelBase_impl.hpp       | 29 ++++++++++++++++++-
 opm/autodiff/BlackoilModelParameters.cpp      |  2 ++
 opm/autodiff/BlackoilModelParameters.hpp      |  3 ++
 opm/autodiff/LinearisedBlackoilResidual.hpp   |  2 ++
 opm/autodiff/NewtonIterationBlackoilCPR.cpp   |  3 +-
 .../NewtonIterationBlackoilInterleaved.cpp    |  3 +-
 7 files changed, 40 insertions(+), 5 deletions(-)

diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp
index 0f6229d3d..970f0f8f3 100644
--- a/opm/autodiff/BlackoilModelBase.hpp
+++ b/opm/autodiff/BlackoilModelBase.hpp
@@ -208,6 +208,9 @@ namespace Opm {
         /// The number of active phases in the model.
         int numPhases() const;
 
+        /// Update the scaling factors for mass balance equations
+        void updateEquationsScaling();
+
     protected:
 
         // ---------  Types and enums  ---------
diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp
index cdea87705..969437075 100644
--- a/opm/autodiff/BlackoilModelBase_impl.hpp
+++ b/opm/autodiff/BlackoilModelBase_impl.hpp
@@ -176,7 +176,8 @@ namespace detail {
         , isSg_(V::Zero(AutoDiffGrid::numCells(grid)))
         , residual_ ( { std::vector<ADB>(fluid.numPhases(), ADB::null()),
                         ADB::null(),
-                        ADB::null() } )
+                        ADB::null(),
+                        { 1.1169, 1.0031, 0.0031 }} ) // the default magic numbers
         , terminal_output_ (terminal_output)
     {
 #if HAVE_MPI
@@ -910,6 +911,32 @@ namespace detail {
             // OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]);
 
         }
+
+
+        if (param_.update_equations_scaling_) {
+            updateEquationsScaling();
+        }
+
+    }
+
+
+
+
+
+    template <class Grid, class Implementation>
+    void
+    BlackoilModelBase<Grid, Implementation>::updateEquationsScaling() {
+        ADB::V B;
+        const Opm::PhaseUsage& pu = fluid_.phaseUsage();
+        for ( int idx=0; idx<MaxNumPhases; ++idx )
+        {
+            if (active_[idx]) {
+                const int pos    = pu.phase_pos[idx];
+                const ADB& tempB = rq_[pos].b;
+                B = 1./tempB.value();
+                residual_.matbalscale[idx] = B.mean();
+            }
+        }
     }
 
 
diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp
index eb5927455..0cefa6139 100644
--- a/opm/autodiff/BlackoilModelParameters.cpp
+++ b/opm/autodiff/BlackoilModelParameters.cpp
@@ -47,6 +47,7 @@ namespace Opm
         tolerance_cnv_   = param.getDefault("tolerance_cnv", tolerance_cnv_);
         tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
         solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_);
+        update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
     }
 
 
@@ -63,6 +64,7 @@ namespace Opm
         tolerance_cnv_   = 1.0e-2;
         tolerance_wells_ = 1.0e-3;
         solve_welleq_initially_ = true;
+        update_equations_scaling_ = false;
     }
 
 
diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp
index 15e7042db..1d1e25525 100644
--- a/opm/autodiff/BlackoilModelParameters.hpp
+++ b/opm/autodiff/BlackoilModelParameters.hpp
@@ -46,6 +46,9 @@ namespace Opm
         /// Solve well equation initially
         bool solve_welleq_initially_;
 
+        /// Update scaling factors for mass balance equations
+        bool update_equations_scaling_;
+
         /// Construct from user parameters or defaults.
         explicit BlackoilModelParameters( const parameter::ParameterGroup& param );
 
diff --git a/opm/autodiff/LinearisedBlackoilResidual.hpp b/opm/autodiff/LinearisedBlackoilResidual.hpp
index d105e3026..acc3e8c83 100644
--- a/opm/autodiff/LinearisedBlackoilResidual.hpp
+++ b/opm/autodiff/LinearisedBlackoilResidual.hpp
@@ -63,6 +63,8 @@ namespace Opm
         /// pressure specification.
         ADB well_eq;
 
+        std::array<double, 3> matbalscale;
+
         /// The size of the non-linear system.
         int sizeNonLinear() const;
     };
diff --git a/opm/autodiff/NewtonIterationBlackoilCPR.cpp b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
index 27d1902dd..a0ed42034 100644
--- a/opm/autodiff/NewtonIterationBlackoilCPR.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilCPR.cpp
@@ -104,9 +104,8 @@ namespace Opm
         }
 
         // Scale material balance equations.
-        const double matbalscale[3] = { 1.1169, 1.0031, 0.0031 }; // HACK hardcoded instead of computed.
         for (int phase = 0; phase < np; ++phase) {
-            eqs[phase] = eqs[phase] * matbalscale[phase];
+            eqs[phase] = eqs[phase] * residual.matbalscale[phase];
         }
 
         // Add material balance equations (or other manipulations) to
diff --git a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
index 8b1451ef2..298136521 100644
--- a/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
+++ b/opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
@@ -102,9 +102,8 @@ namespace Opm
         }
 
         // Scale material balance equations.
-        const double matbalscale[3] = { 1.1169, 1.0031, 0.0031 }; // HACK hardcoded instead of computed.
         for (int phase = 0; phase < np; ++phase) {
-            eqs[phase] = eqs[phase] * matbalscale[phase];
+            eqs[phase] = eqs[phase] * residual.matbalscale[phase];
         }
 
         // calculating the size for b