From ef053a188346eb21f093ec5639ff6712a6dd44d6 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 8 Oct 2012 10:29:45 +0200 Subject: [PATCH] Cleaned up single cell Newton solver. --- .../TransportModelCompressiblePolymer.cpp | 236 +----------------- 1 file changed, 12 insertions(+), 224 deletions(-) diff --git a/opm/polymer/TransportModelCompressiblePolymer.cpp b/opm/polymer/TransportModelCompressiblePolymer.cpp index 99d630938..b9abe0afd 100644 --- a/opm/polymer/TransportModelCompressiblePolymer.cpp +++ b/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 dFx_dx = (dres_s_dsdc[0]-x_c[1]*dres_s_dsdc[1]); - double dFx_dy; - if (x[0] < 1e-2*tol_) { - dFx_dy = 0.0; + double x_c_app[2]; + // The computation of the Jacobi fails for s=0 (we have an undetermined fraction 0/0). + // When s is close to zero we replace x_c with x_c_app. + x_c_app[1] = x_c[1]; + if (x_c[0] < 1e-2*tol_) { + x_c_app[0] = 1e-2*tol_; } else { - dFx_dy = (dres_s_dsdc[1]/x[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]); + x_c_app[0] = x_c[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_it0) ? 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