update scaling factors for mass balance equations.

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  ---------

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();
+            }
+        }
     }
 
 

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;
     }
 
 

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 );
 

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;
     };

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

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