diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index 105542b66..ed3d528d7 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -495,19 +495,14 @@ main(int argc, char** argv) c_vals_visc[1] = 7.0; std::vector visc_mult_vals(2, -1e100); visc_mult_vals[0] = 1.0; - // polyprop.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0); - visc_mult_vals[1] = 1.0; - std::vector c_vals_ads(3, -1e100); + visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0); + std::vector c_vals_ads(2, -1e100); c_vals_ads[0] = 0.0; - c_vals_ads[1] = 2.0; - c_vals_ads[2] = 8.0; + c_vals_ads[1] = 8.0; // Here we set up adsorption equal to zero. - std::vector ads_vals(3, -1e100); + std::vector ads_vals(2, -1e100); ads_vals[0] = 0.0; ads_vals[1] = 0.0; - ads_vals[2] = 0.0; - // ads_vals[1] = 0.0; - // ads_vals[2] = 0.0; polyprop.set(c_max, mix_param, rock_density, dead_pore_vol, res_factor, c_max_ads, static_cast(ads_index), c_vals_visc, visc_mult_vals, c_vals_ads, ads_vals); @@ -609,6 +604,10 @@ main(int argc, char** argv) method = Opm::TransportModelPolymer::Newton; } else if (method_string == "Gradient") { method = Opm::TransportModelPolymer::Gradient; + } else if (method_string == "NewtonSimpleSC") { + method = Opm::TransportModelPolymer::NewtonSimpleSC; + } else if (method_string == "NewtonSimpleC") { + method = Opm::TransportModelPolymer::NewtonSimpleC; } else { THROW("Unknown method: " << method_string); } @@ -767,6 +766,11 @@ main(int argc, char** argv) } } while (!well_control_passed); + // Update pore volumes if rock is compressible. + if (rock_comp->isActive()) { + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + } + // Process transport sources (to include bdy terms and well flows). Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, wells->c_wells(), well_state.perfRates(), reorder_src); diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp index e7289ce39..d48067b52 100644 --- a/opm/polymer/SimulatorPolymer.cpp +++ b/opm/polymer/SimulatorPolymer.cpp @@ -296,6 +296,11 @@ namespace Opm ptime += pt; } while (false); + // Update pore volumes if rock is compressible. + if (rock_comp_props_ && rock_comp_props_->isActive()) { + computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol); + } + // Process transport sources (to include bdy terms and well flows). Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0, wells_, well_state.perfRates(), transport_src); diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp index a1cef6273..544c96f34 100644 --- a/opm/polymer/TransportModelPolymer.cpp +++ b/opm/polymer/TransportModelPolymer.cpp @@ -26,7 +26,7 @@ #include #include #include - +#include // Choose error policy for scalar solves here. typedef Opm::RegulaFalsi RootFinder; @@ -114,8 +114,11 @@ namespace return std::max(std::abs(res[0]), std::abs(res[1])); } - bool solveNewtonStep(const double* , const Opm::TransportModelPolymer::ResidualEquation&, - const double*, double*); + bool solveNewtonStepSC(const double* , const Opm::TransportModelPolymer::ResidualEquation&, + const double*, double*); + bool solveNewtonStepC(const double* , const Opm::TransportModelPolymer::ResidualEquation&, + const double*, double*); + // Define a piecewise linear curve along which we will look for zero of the "s" or "r" residual. @@ -333,7 +336,7 @@ namespace Opm comp_term = tm.source_[cell]; // Note: this assumes that all source flux is water. dtpv = tm.dt_/tm.porevolume_[cell]; porosity = tm.porosity_[cell]; - s = -1e100; + s = s0; for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { int f = tm.grid_.cell_faces[i]; @@ -388,7 +391,7 @@ namespace Opm // Solve for s first. // s = modifiedRegulaFalsi(res_s, std::max(tm.smin_[2*cell], dps), tm.smax_[2*cell], // tm.maxit_, tm.tol_, iters_used); - s = RootFinder::solve(res_s, s0, 0.0, 1.0, + s = RootFinder::solve(res_s, s, 0.0, 1.0, tm.maxit_, tm.tol_, iters_used); double ff; tm.fracFlow(s, c, cmax0, cell, ff); @@ -641,6 +644,12 @@ namespace Opm case Gradient: solveSingleCellGradient(cell); break; + case NewtonSimpleSC: + solveSingleCellNewtonSimple(cell,true); + break; + case NewtonSimpleC: + solveSingleCellNewtonSimple(cell,false); + break; default: THROW("Unknown method " << method_); } @@ -857,8 +866,8 @@ namespace Opm // res_eq.computeResidual(x, res, mc, ff); } - double x_min[2] = { 0.0, 0.0 }; - double x_max[2] = { 1.0, polyprops_.cMax()*1.1 }; + const double x_min[2] = { 0.0, 0.0 }; + const double x_max[2] = { 1.0, polyprops_.cMax()*1.1 }; bool successfull_newton_step = true; // initialize x_new to avoid warning double x_new[2] = {0.0, 0.0}; @@ -932,6 +941,147 @@ namespace Opm } } + void TransportModelPolymer::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]; + }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()*1.1 }; + 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=1e-6; + 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("Newton 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 TransportModelPolymer::solveMultiCell(const int num_cells, const int* cells) { @@ -1471,12 +1621,12 @@ namespace return res_eq_.computeResidualC(x_c); } - bool solveNewtonStep(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq, + bool solveNewtonStepSC(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq, const double* res, double* x_new) { double dres_s_dsdc[2]; double dres_c_dsdc[2]; - double dx=1e-11; + double dx=1e-8; double x[2]; if(!(x[0]>0)){ x[0] = dx; @@ -1502,6 +1652,30 @@ namespace return true; } } + bool solveNewtonStepC(const double* xx, const Opm::TransportModelPolymer::ResidualEquation& res_eq, + const double* res, double* x_new) { + + 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]; + } + res_eq.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; + } + } } // Anonymous namespace diff --git a/opm/polymer/TransportModelPolymer.hpp b/opm/polymer/TransportModelPolymer.hpp index c8e6adf44..198e72f7f 100644 --- a/opm/polymer/TransportModelPolymer.hpp +++ b/opm/polymer/TransportModelPolymer.hpp @@ -40,7 +40,7 @@ namespace Opm { public: - enum SingleCellMethod { Bracketing, Newton, Gradient }; + enum SingleCellMethod { Bracketing, Newton, Gradient, NewtonSimpleSC, NewtonSimpleC}; enum GradientMethod { Analytic, FinDif }; // Analytic is chosen (hard-coded) /// Construct solver. @@ -106,6 +106,7 @@ namespace Opm void solveSingleCellBracketing(int cell); void solveSingleCellNewton(int cell); void solveSingleCellGradient(int cell); + void solveSingleCellNewtonSimple(int cell,bool use_sc); class ResidualEquation; void initGravity(const double* grav);