Rewriting the detection function with for loop.

It is more automatic while it remains to change to some more flexible
form.

examples/sim_fibo_ad.cpp

@@ -222,6 +222,7 @@ try
         SimulatorReport episodeReport = simulator.run(simtimer, state, well_state);
 
         ++simtimer;
+        std::cin.ignore();
 
         outputWriter.writeTimeStep(simtimer, state, well_state.basicWellState());
         fullReport += episodeReport;

opm/autodiff/FullyImplicitBlackoilSolver_impl.hpp

@@ -263,7 +263,7 @@ namespace {
             computeWellConnectionPressures(state, xw);
         }
 
-        const int    maxit = 15;
+        const int    maxit = 25;
 
         std::vector<std::vector<double>> residual_history;
 
@@ -285,22 +285,22 @@ namespace {
                   << std::setw(9) << it << std::setprecision(9)
                   << std::setw(18) << r0 << std::endl;
 
-        // compute the size of the linear system
+        // compute the size of the non-linear system
-        int linearSize = 0;
+        int sizeNonLinear = 0;
 
         std::vector<ADB>::const_iterator massBalanceIt = residual_.material_balance_eq.begin();
         const std::vector<ADB>::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
 
         for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) {
-            linearSize += (*massBalanceIt).size();
+            sizeNonLinear += (*massBalanceIt).size();
         }
-        linearSize += residual_.well_flux_eq.size();
+        sizeNonLinear += residual_.well_flux_eq.size();
-        linearSize += residual_.well_eq.size();
+        sizeNonLinear += residual_.well_eq.size();
 
-        std::cout << " the size of the linear system is " << linearSize << std::endl;
+        std::cout << " the size of the linear system is " << sizeNonLinear<< std::endl;
         // std::cin.ignore();
 
-        V dxOld = V::Zero(linearSize);
+        V dxOld = V::Zero(sizeNonLinear);
 
         bool isOscillate = false;
         bool isStagnate = false;
@@ -1704,66 +1704,53 @@ namespace {
         // The detection of oscillation in two primary variable results in the report of the detection
         // of oscillation for the solver 
         // Stagnate is not used for any treatment here.
-        const Opm::PhaseUsage& pu = fluid_.phaseUsage();
 
         oscillate = false;
         stagnate = false;
 
-        if ( it < 3 ) {
+        if ( it < 2 ) {
             return;
         }
 
-        bool oscillateWater = false;
+        int oscillatePhase = 0;
-        bool oscillateOil  = false;
-        bool oscillateGas = false;
 
-        bool stagnateWater = true;
+        stagnate = true;
-        bool stagnateOil = true;
-        bool stagnateGas = true;
 
+        // bool oscillateWater = false;
+        // bool oscillateOil  = false;
+        // bool oscillateGas = false;
 
-        if ( active_[Water]) {
+        // bool stagnateWater = true;
-            const int pos = pu.phase_pos[Water];
+        // bool stagnateOil = true;
-            double relChange1 = std::fabs((residual_history[it][Water] - residual_history[it - 2][Water]) /
+        // bool stagnateGas = true;
-                                           residual_history[it][Water] );
-            double relChange2 = std::fabs((residual_history[it][Water] - residual_history[it - 1][Water]) /
-                                           residual_history[it][Water] );
-            oscillateWater = (relChange1 < relaxRelTol) && (relChange2 > relaxRelTol);
 
-            double relChange3 = std::fabs((residual_history[it - 1][Water] - residual_history[it - 2][Water]) /
+        std::cout << " residual_history " << std::endl;
-                                           residual_history[it - 2][Water] );
+        for( int i = it-2; i <=it; i++ ){
-            stagnateWater = relChange3 < 1.e-3;
+            for ( int j = 0; j<5; j++ ){
+                std::cout << " " << residual_history[i][j];
+            }
+            std::cout << std::endl;
+        }
+        for (int phaseIdx= 0; phaseIdx < fluid_.numPhases(); ++ phaseIdx){
+            if (active_[phaseIdx]) {
+                double relChange1 = std::fabs((residual_history[it][phaseIdx] - residual_history[it - 2][phaseIdx]) /
+                                               residual_history[it][phaseIdx]);
+                double relChange2 = std::fabs((residual_history[it][phaseIdx] - residual_history[it - 1][phaseIdx]) /
+                                               residual_history[it][phaseIdx]);
+                oscillatePhase += (relChange1 < relaxRelTol) && (relChange2 > relaxRelTol);
+
+                double relChange3 = std::fabs((residual_history[it - 1][phaseIdx] - residual_history[it - 2][phaseIdx]) /
+                                               residual_history[it - 2][phaseIdx]);
+                stagnate = stagnate && (relChange3 > 1.e-3);
+
+                std::cout << " relChange1 " << relChange1 << " relChange2 " << relChange2
+                          << " relChange3 " << relChange3 << std::endl;
+            }
         }
 
-        if ( active_[Oil]) {
+        oscillate = (oscillatePhase > 1);
-            const int pos = pu.phase_pos[Oil];
-            double relChange1 = std::fabs((residual_history[it][Oil] - residual_history[it - 2][Oil]) /
-                                           residual_history[it][Oil] );
-            double relChange2 = std::fabs((residual_history[it][Oil] - residual_history[it - 1][Oil]) /
-                                           residual_history[it][Oil] );
-            oscillateOil = (relChange1 < relaxRelTol) && (relChange2 > relaxRelTol);
-
-            double relChange3 = std::fabs((residual_history[it - 1][Oil] - residual_history[it - 2][Oil]) /
-                                           residual_history[it - 2][Oil] );
-            stagnateOil = relChange3 < 1.e-3;
-        }
-
-        if ( active_[Gas]) {
-            const int pos = pu.phase_pos[Gas];
-            double relChange1 = std::fabs((residual_history[it][Gas] - residual_history[it - 2][Gas]) /
-                                           residual_history[it][Gas] );
-            double relChange2 = std::fabs((residual_history[it][Gas] - residual_history[it - 1][Gas]) /
-                                           residual_history[it][Gas] );
-            oscillateGas = (relChange1 < relaxRelTol) && (relChange2 > relaxRelTol);
-
-            double relChange3 = std::fabs((residual_history[it - 1][Gas] - residual_history[it - 2][Gas]) /
-                                           residual_history[it - 2][Gas] );
-            stagnateGas = relChange3 < 1.e-3;
-        }
-
-        stagnate = stagnateWater && stagnateOil && stagnateGas;
-        oscillate = ( oscillateWater + oscillateOil + oscillateGas ) > 1;
 
+        std::cin.ignore();
     }