Cleaned up single cell Newton solver.

opm/polymer/TransportModelCompressiblePolymer.cpp

@@ -65,8 +65,6 @@ public:
     void computeGradientResS(const double* x, double* res, double* gradient) const;
     void computeGradientResC(const double* x, double* res, double* gradient) const;
     void computeJacobiRes(const double* x, double* dres_s_dsdc, double* dres_c_dsdc) const;
-    bool solveNewtonStepSC(const double*, const double*, double*) const;
-    bool solveNewtonStepC(const double* , const double*, double*) const;
 
     
 
@@ -635,12 +633,6 @@ namespace Opm
 	case Gradient:
 	    solveSingleCellGradient(cell);
 	    break;
-	case NewtonSimpleSC:
-	    solveSingleCellNewtonSimple(cell,true);
-	    break;
-	case NewtonSimpleC:
-	    solveSingleCellNewtonSimple(cell,false);
-	    break;	    
 	default:
 	    THROW("Unknown method " << method_);
 	}
@@ -879,21 +871,20 @@ namespace Opm
             double dres_s_dsdc[2];
             double dres_c_dsdc[2];
             scToc(x, x_c);
-            res_eq.computeJacobiRes(x_c, dres_s_dsdc, dres_c_dsdc);
+            double x_c_app[2];
-            double dFx_dx = (dres_s_dsdc[0]-x_c[1]*dres_s_dsdc[1]);
+            // The computation of the Jacobi fails for s=0 (we have an undetermined fraction 0/0).
-            double dFx_dy;
+            // When s is close to zero we replace x_c with x_c_app.
-            if (x[0] < 1e-2*tol_) {
+            x_c_app[1] = x_c[1];
-                dFx_dy = 0.0;
+            if (x_c[0] < 1e-2*tol_) {
+                x_c_app[0] = 1e-2*tol_;
             } else {
-                dFx_dy = (dres_s_dsdc[1]/x[0]);
+                x_c_app[0] = x_c[0];
-            }
-            double dFy_dx = (dres_c_dsdc[0]-x_c[1]*dres_c_dsdc[1]);
-            double dFy_dy;
-            if (x[0] < 1e-2*tol_) {
-                dFy_dy = 1.0 - res_eq.dps;
-            } else {
-                dFy_dy = (dres_c_dsdc[1]/x[0]);
             }
+            res_eq.computeJacobiRes(x_c_app, dres_s_dsdc, dres_c_dsdc);
+            double dFx_dx = (dres_s_dsdc[0]-x_c_app[1]*dres_s_dsdc[1]);
+            double dFx_dy = (dres_s_dsdc[1]/x_c_app[0]);
+            double dFy_dx = (dres_c_dsdc[0]-x_c_app[1]*dres_c_dsdc[1]);
+            double dFy_dy = (dres_c_dsdc[1]/x_c_app[0]);
             double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
             double alpha = 1.0;
             int max_lin_it = 100;
@@ -934,154 +925,6 @@ namespace Opm
 	}
     }
 
-    void TransportModelCompressiblePolymer::solveSingleCellNewtonSimple(int cell,bool use_sc)
-    {
-	const int max_iters_split = maxit_;
-	int iters_used_split = 0;
-
-	// Check if current state is an acceptable solution.
-	ResidualEquation res_eq(*this, cell);
-	double x[2] = {saturation_[cell], concentration_[cell]};
-	double res[2];
-	double mc;
-	double ff;
-	res_eq.computeResidual(x, res, mc, ff);
-	if (norm(res) <= tol_) {
-	    cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
- 	    fractionalflow_[cell] = ff;
-	    mc_[cell] = mc;
-	    return;
-	}else{
-	    //*
-	    x[0] = saturation_[cell]-res[0];
-	    if((x[0]>1) || (x[0]<0)){
-		x[0] = 0.5;
-	        x[1] = x[1];
-	    }
-	    if(x[0]>0){
-                x[1] =  concentration_[cell]*saturation_[cell]-res[1];
-                x[1] = x[1]/x[0];
-                if(x[1]> polyprops_.cMax()){
-                    x[1]= polyprops_.cMax()/2.0;		
-                }
-                if(x[1]<0){
-                    x[1]=0;
-                }
-	    }else{
-		x[1]=0;
-	    }	    	    
-	    //x[0]=0.5;x[1]=polyprops_.cMax()/2.0;
-	    res_eq.computeResidual(x, res, mc, ff);
-	    //*/
-	}
-
-
-	// double x_min[2] = { std::max(polyprops_.deadPoreVol(), smin_[2*cell]), 0.0 };
-        double x_min[2] = { 0.0, 0.0 };
-	double x_max[2] = { 1.0, polyprops_.cMax()*adhoc_safety_ };
-	bool successfull_newton_step = true;
-	double x_new[2];
-	double res_new[2];
-	ResSOnCurve res_s_on_curve(res_eq);
-	ResCOnCurve res_c_on_curve(res_eq);
-	
- 	while ((norm(res) > tol_) &&
-	       (iters_used_split < max_iters_split)  &&
-	       successfull_newton_step) {
-	    // We first try a Newton step
-	    double dres_s_dsdc[2];
-	    double dres_c_dsdc[2];
-	    double dx=tol_;
-	    double tmp_x[2];
-	    if(!(x[0]>0)){
-		tmp_x[0]=dx;
-		tmp_x[1]=0;
-	    }else{
-		tmp_x[0]=x[0];
-		tmp_x[1]=x[1];
-	    }
-	    res_eq.computeJacobiRes(tmp_x, dres_s_dsdc, dres_c_dsdc);
-	    double dFx_dx,dFx_dy,dFy_dx,dFy_dy;
-	    double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
-	    if(use_sc){
-		dFx_dx=(dres_s_dsdc[0]-tmp_x[1]*dres_s_dsdc[1]/tmp_x[0]);
-		dFx_dy= (dres_s_dsdc[1]/tmp_x[0]);
-		dFy_dx=(dres_c_dsdc[0]-tmp_x[1]*dres_c_dsdc[1]/tmp_x[0]);
-		dFy_dy= (dres_c_dsdc[1]/tmp_x[0]);
-	    }else{
-		dFx_dx= dres_s_dsdc[0];
-		dFx_dy= dres_s_dsdc[1];
-		dFy_dx= dres_c_dsdc[0];
-		dFy_dy= dres_c_dsdc[1];
-	    }
-	    det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
-	    if(det==0){
-		THROW("det is zero");
-	    }
-
-	    double alpha=1;
-	    int max_lin_it=10;
-	    int lin_it=0;
-	    /*
-	      std::cout << "Nonlinear " << iters_used_split
-	      << "  " << norm(res_new)
-	      << "  " << norm(res) << std::endl;*/
-	    /*
-	      std::cout  << "x" << std::endl;
-	      std::cout  << " " <<  x[0] << " " <<  x[1] << std::endl;
-	      std::cout  << "dF" << std::endl;
-	      std::cout  << " " <<  dFx_dx << " " <<  dFx_dy << std::endl;
-	      std::cout  << " " <<  dFy_dx << " " <<  dFy_dy << std::endl;
-	      std::cout  << "dres" << std::endl;
-	      std::cout  << " " <<  dres_s_dsdc[0] << " " <<  dres_s_dsdc[1] << std::endl;
-	      std::cout  << " " <<  dres_c_dsdc[0] << " " <<  dres_c_dsdc[1] << std::endl;
-	      std::cout  << std::endl;
-	    */
-	    res_new[0]=res[0]*2;
-	    res_new[1]=res[1]*2;
-	    double update[2]={(res[0]*dFy_dy - res[1]*dFx_dy)/det,
-			      (res[1]*dFx_dx - res[0]*dFy_dx)/det};
-	    while((norm(res_new)>norm(res)) && (lin_it<max_lin_it)) {
-		x_new[0] = x[0] - alpha*update[0];
-		if(use_sc){
-		    x_new[1] = x[1]*x[0] - alpha*update[1];
-		    x_new[1] = (x_new[0]>0) ?  x_new[1]/x_new[0] : 0.0;
-		}else{
-		    x_new[1] = x[1] - alpha*update[1];
-		}
-		check_interval(x_min, x_max, x_new);
-		res_eq.computeResidual(x_new, res_new, mc, ff);
-		//  std::cout << " " << res_new[0] << "  " << res_new[1] << std::endl;
-		alpha=alpha/2.0;
-		lin_it=lin_it+1;
-		//  std::cout << "Linear iterations" << lin_it << "  " << norm(res_new) << std::endl;
-	    }
-	    if (lin_it>=max_lin_it) {
-		successfull_newton_step = false;
-	    } else  {
-		x[0] = x_new[0];
-		x[1] = x_new[1];
-		res[0] = res_new[0];
-		res[1] = res_new[1];
-		iters_used_split += 1;
-		successfull_newton_step = true;;		    
-	    }
-	    //	    std::cout << "Nonlinear " << iters_used_split << "  " << norm(res) << std::endl;
-	}
-		
-	if ((iters_used_split >=  max_iters_split) || (norm(res) > tol_)) {
-	    MESSAGE("NewtonSimple for single cell did not work in cell number " << cell);
-	    solveSingleCellBracketing(cell);
-	} else {
-	    concentration_[cell] = x[1];
-	    cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
-	    saturation_[cell] = x[0];
-	    fractionalflow_[cell] = ff;
-	    mc_[cell] = mc;
-	}
-    }
-
-
 
     void TransportModelCompressiblePolymer::solveMultiCell(const int num_cells, const int* cells)
     {
@@ -1526,61 +1369,6 @@ namespace Opm
         }
     }
 
-    bool TransportModelCompressiblePolymer::ResidualEquation::solveNewtonStepSC(const double* xx, 
-                           const double* res, double* x_new) const {
-
-    	double dres_s_dsdc[2];
-    	double dres_c_dsdc[2];
-	double dx=1e-8;
-	double x[2];
-	if(!(x[0]>0)){
-	    x[0] = dx;
-            x[1] = 0;
-	} else {
-	    x[0] = xx[0];
-            x[1] = xx[1];
-	}
-	computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
-	
-	double dFx_dx=(dres_s_dsdc[0]-x[1]*dres_s_dsdc[1]);
-	double dFx_dy= (dres_s_dsdc[1]/x[0]);
-	double dFy_dx=(dres_c_dsdc[0]-x[1]*dres_c_dsdc[1]);
-	double dFy_dy= (dres_c_dsdc[1]/x[0]);
-
-	double det = dFx_dx*dFy_dy - dFy_dx*dFx_dy;
-    	if (std::abs(det) < 1e-8) {
-    	    return false;
-    	} else {
-	    x_new[0] = x[0] - (res[0]*dFy_dy - res[1]*dFx_dy)/det;
-	    x_new[1] = x[1]*x[0] - (res[1]*dFx_dx - res[0]*dFy_dx)/det;
-	    x_new[1] = (x_new[0]>0) ?  x_new[1]/x_new[0] : 0.0;
-	    return true;
-	}
-    }
-    bool TransportModelCompressiblePolymer::ResidualEquation::solveNewtonStepC(const double* xx, const double* res, double* x_new) const {
-	
-	double dres_s_dsdc[2];
-	double dres_c_dsdc[2];
-        double dx=1e-8;
-        double x[2];
-        if(!(x[0]>0)){
-            x[0] = dx;
-            x[1] = 0;
-        } else {
-	    x[0] = xx[0];
-            x[1] = xx[1];
-        }
-        computeJacobiRes(x, dres_s_dsdc, dres_c_dsdc);
-        double det = dres_s_dsdc[0]*dres_c_dsdc[1] - dres_c_dsdc[0]*dres_s_dsdc[1];
-	if (std::abs(det) < 1e-8) {
-    	    return false;
-    	} else {
-    	    x_new[0] = x[0] - (res[0]*dres_c_dsdc[1] - res[1]*dres_s_dsdc[1])/det;
-	    x_new[1] = x[1] - (res[1]*dres_s_dsdc[0] - res[0]*dres_c_dsdc[0])/det;
-	    return true;
-	}
-    }
-
 
 } // namespace Opm