From c63d81733292717c698b378ed1b5f26de735e2d6 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Fri, 24 Feb 2012 17:35:47 +0100 Subject: [PATCH] Added varying bounded box for allowable values of c and c in the splitting residual solver. --- opm/polymer/TransportModelPolymer.cpp | 268 +++++++++++++++++++++++--- opm/polymer/TransportModelPolymer.hpp | 1 + 2 files changed, 245 insertions(+), 24 deletions(-) diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp index 4e5ed521c..e6482399b 100644 --- a/opm/polymer/TransportModelPolymer.cpp +++ b/opm/polymer/TransportModelPolymer.cpp @@ -424,7 +424,7 @@ namespace Opm // parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding // rectangle, x_out=(s_out, c_out) denotes the values of s and c at that point. - void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, double& t_max_arg) + void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, double& t_max_out, double& t_out_out) { double t0, t1; x[0] = x_arg[0]; @@ -455,7 +455,8 @@ namespace Opm x_out[0] = x[0] + t_out*direction[0]; x_out[1] = x[1] + t_out*direction[1]; t_max = t_out + 1; - t_max_arg = t_max; + t_max_out = t_max; + t_out_out = t_out; } @@ -465,8 +466,8 @@ namespace Opm x_new[0] = x[0] + t*direction[0]; x_new[1] = x[1] + t*direction[1]; } else { - x_new[0] = (t_max - t)*x_out[0] + end_point[0]*(t - t_out); - x_new[1] = (t_max - t)*x_out[1] + end_point[1]*(t - t_out); + x_new[0] = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out)); + x_new[1] = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out)); } } @@ -480,10 +481,10 @@ namespace Opm s = x[0] + t*direction[0]; c = x[1] + t*direction[1]; } else { - s = (t_max - t)*x_out[0] + end_point[0]*(t - t_out); - c = (t_max - t)*x_out[1] + end_point[1]*(t - t_out); + s = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out)); + c = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out)); } - return s - s0 + dtpv*(outflux*tm.fracFlow(s, c, cell) + influx); + return s - s0 + dtpv*(outflux*tm.fracFlow(s, c, cell) + influx); } @@ -555,12 +556,12 @@ namespace Opm x_new[0] = x[0] + t*direction[0]; x_new[1] = x[1] + t*direction[1]; } else { - x_new[0] = (t_max - t)*x_out[0] + end_point[0]*(t - t_out); - x_new[1] = (t_max - t)*x_out[1] + end_point[1]*(t - t_out); + x_new[0] = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out)); + x_new[1] = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out)); } } - void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, double& t_max_arg) + void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, double& t_max_out, double& t_out_out) { bool if_t0 = true; bool if_t1 = true; @@ -609,7 +610,8 @@ namespace Opm x_out[0] = x[0] + t_out*direction[0]; x_out[1] = x[1] + t_out*direction[1]; t_max = t_out + 1; - t_max_arg = t_max; + t_max_out = t_max; + t_out_out = t_out; } double operator()(double t) const @@ -620,8 +622,8 @@ namespace Opm s = x[0] + t*direction[0]; c = x[1] + t*direction[1]; } else { - s = (t_max - t)*x_out[0] + end_point[0]*(t - t_out); - c = (t_max - t)*x_out[1] + end_point[1]*(t - t_out); + s = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out)); + c = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out)); } double ff = tm.fracFlow(s, c, cell); double mc = tm.computeMc(c); @@ -635,6 +637,171 @@ namespace Opm } }; + struct TransportModelPolymer::ResidualDir + { + int cell; + int s_or_c; // s_or_c = 0 if s direction, s_or_c = 1 if c direction, + double s0; + double c0; + double cmax0; + double influx; // sum_j min(v_ij, 0)*f(s_j) + double influx_polymer; // sum_j min(v_ij, 0)*f(s_j)*mc(c_j) + double outflux; // sum_j max(v_ij, 0) + double porosity; + double dtpv; // dt/pv(i) + double direction[2]; + double end_point[2]; + double x_max[2]; + double x_min[2]; + double t_out; + double t_max; // t_max = t_out + 1 + double x_out[2]; + double x[2]; + const TransportModelPolymer& tm; + + ResidualDir(const TransportModelPolymer& tmodel, int cell_index) + : tm(tmodel) + { + cell = cell_index; + s0 = tm.saturation_[cell]; + c0 = tm.concentration_[cell]; + cmax0 = tm.cmax_[cell]; + double dflux = -tm.source_[cell]; + bool src_is_inflow = dflux < 0.0; + influx = src_is_inflow ? dflux : 0.0; + influx_polymer = src_is_inflow ? dflux*tm.computeMc(tm.inflow_c_) : 0.0; + outflux = !src_is_inflow ? dflux : 0.0; + dtpv = tm.dt_/tm.porevolume_[cell]; + porosity = tm.porosity_[cell]; + + for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) { + int f = tm.grid_.cell_faces[i]; + double flux; + int other; + // Compute cell flux + if (cell == tm.grid_.face_cells[2*f]) { + flux = tm.darcyflux_[f]; + other = tm.grid_.face_cells[2*f+1]; + } else { + flux =-tm.darcyflux_[f]; + other = tm.grid_.face_cells[2*f]; + } + // Add flux to influx or outflux, if interior. + if (other != -1) { + if (flux < 0.0) { + influx += flux*tm.fractionalflow_[other]; + influx_polymer += flux*tm.fractionalflow_[other]*tm.mc_[other]; + } else { + outflux += flux; + } + } + } + } + + + // For a given point x=(s,c) in the s,c plane, set up a piecewise linear curve wich starts + // from "x" with slope "direction", hits the bound of the rectangle + // [s_min,s_max]x[c_min,c_max] and continue in a straight line to "end_point". The curve is + // parametrized by t in [0, t_max], t_out is equal to t when the curve hits the bounding + // rectangle, x_out=(s_out, c_out) denotes the values of s and c at that point. + + void setup(const double* x_arg, const double* direction_arg, const double* end_point_arg, const double* x_min_arg, const double* x_max_arg, const int& s_or_c_arg, double& t_max_out, double& t_out_out) + { + bool if_t0 = true; + bool if_t1 = true; + double t0, t1; + + s_or_c = s_or_c_arg; + x[0] = x_arg[0]; + x[1] = x_arg[1]; + x_max[0] = x_max_arg[0]; + x_max[1] = x_max_arg[1]; + x_min[0] = x_min_arg[0]; + x_min[1] = x_min_arg[1]; + direction[0] = direction_arg[0]; + direction[1] = direction_arg[1]; + end_point[0] = end_point_arg[0]; + end_point[1] = end_point_arg[1]; + if ((end_point[0]-x[0])*direction[0] + (end_point[1]-x[1])*direction[1] < 0) { + direction[0] *= -1; + direction[1] *= -1; + } + if (direction[0] == 0) { + if_t0 = false; + } else { + if (direction[0] > 0) { + t0 = (x_max[0] - x[0])/direction[0]; + } else { + t0 = (x_min[0] - x[0])/direction[0]; + } + } + if (direction[1] == 0) { + if_t1 = false; + } else { + if (direction[1] > 0) { + t1 = (x_max[1] - x[1])/direction[1]; + } else { + t1 = (x_min[1] - x[1])/direction[1]; + } + } + if (if_t0 && if_t1) { + t_out = std::min(t0, t1); + } else { + if (if_t0) { + t_out = t0; + } else { + t_out = t1; + } + } + x_out[0] = x[0] + t_out*direction[0]; + x_out[1] = x[1] + t_out*direction[1]; + t_max = t_out + 1; + t_max_out = t_max; + t_out_out = t_out; + } + + + // Compute x=(s,c) for a given t (t is the parameter for the piecewise linear curve) + void compute_new_x(double* x_new, const double t) { + if (t <= t_out) { + x_new[0] = x[0] + t*direction[0]; + x_new[1] = x[1] + t*direction[1]; + } else { + x_new[0] = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out)); + x_new[1] = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out)); + } + } + + + + double operator()(double t) const + { + double s; + double c; + if (t <= t_out) { + s = x[0] + t*direction[0]; + c = x[1] + t*direction[1]; + } else { + s = 1/(t_max-t_out)*((t_max - t)*x_out[0] + end_point[0]*(t - t_out)); + c = 1/(t_max-t_out)*((t_max - t)*x_out[1] + end_point[1]*(t - t_out)); + } + if (s_or_c == 0) { + return s - s0 + dtpv*(outflux*tm.fracFlow(s, c, cell) + influx); + } else if (s_or_c == 1) { + double ff = tm.fracFlow(s, c, cell); + double mc = tm.computeMc(c); + double dps = tm.polyprops_.dps; + double rhor = tm.polyprops_.rhor; + double ads0 = tm.polyprops_.adsorbtion(std::max(c0, cmax0)); + double ads = tm.polyprops_.adsorbtion(std::max(c, cmax0)); + return (s - dps)*c - (s0 - dps)*c0 + + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + + dtpv*(outflux*ff*mc + influx_polymer); + } else { + std::cout << "problem!" << std::endl; + } + } + }; void TransportModelPolymer::solveSingleCell(const int cell) { @@ -667,7 +834,7 @@ namespace Opm void TransportModelPolymer::solveSingleCellSplitting(int cell) { const int max_iters_falsi = 20; - const double tol = 1e-9; + const double tol = 1e-7; int iters_used_falsi = 0; const int max_iters_split = 20; int iters_used_split = 0; @@ -675,6 +842,9 @@ namespace Opm Residual residual(*this, cell); ResidualSDir residual_s_dir(*this, cell); ResidualCDir residual_c_dir(*this, cell); + // const int sdir = 0; + // const int cdir = 1; + // ResidualDir residual_dir(*this, cell); double x[2] = {saturation_[cell], concentration_[cell]}; double res[2]; residual.computeResidual(x, res); @@ -688,6 +858,7 @@ namespace Opm double x_max[2] = {1.0, polyprops_.c_max_limit}; double t; double t_max; + double t_out; double direction[2]; double end_point[2]; double gradient[2]; @@ -699,15 +870,24 @@ namespace Opm end_point[1] = x_min[1]; direction[0] = end_point[0] - x[0]; direction[1] = end_point[1] - x[1]; - residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_s_dir(t_out) >= 0) { + t_max = t_out; + } } else { end_point[0] = x_min[0]; end_point[1] = x_max[1]; direction[0] = end_point[0] - x[0]; direction[1] = end_point[1] - x[1]; - residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_s_dir(t_out) <= 0) { + t_max = t_out; + } } t = modifiedRegulaFalsi(residual_s_dir, 0., t_max, max_iters_falsi, tol, iters_used_falsi); + if (std::abs(residual_s_dir(t)) > tol) { + std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl; + } residual_s_dir.compute_new_x(x, t); } res_s_done = true; @@ -719,15 +899,24 @@ namespace Opm end_point[1] = x_max[1]; direction[0] = end_point[0] - x[0]; direction[1] = end_point[1] - x[1]; - residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_c_dir(t_out) >= 0) { + t_max = t_out; + } } else { end_point[0] = x_min[0]; end_point[1] = x_min[1]; direction[0] = end_point[0] - x[0]; direction[1] = end_point[1] - x[1]; - residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_c_dir(t_out) <= 0) { + t_max = t_out; + } } t = modifiedRegulaFalsi(residual_c_dir, 0., t_max, max_iters_falsi, tol, iters_used_falsi); + if (std::abs(residual_c_dir(t)) > tol) { + std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl; + } residual_c_dir.compute_new_x(x, t); } res_s_done = false; @@ -739,15 +928,37 @@ namespace Opm direction[0] = -gradient[1]; direction[1] = gradient[0]; if (res[1] < 0) { + // We update the bounding box (Here we assume that the curve res_s(s,c)=0 is + // increasing). We do it only for a significantly large res[1] + if (res[1] < -tol ) { + x_min[0] = x[0]; + x_min[1] = x[1]; + } + // end_point[0] = x_max[0]; end_point[1] = x_max[1]; - residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_c_dir(t_out) >= 0) { + t_max = t_out; + } } else { + // We update the bounding box (Here we assume that the curve res_s(s,c)=0 is increasing) + if (res[1] > tol) { + x_max[0] = x[0]; + x_max[1] = x[1]; + } + // end_point[0] = x_min[0]; end_point[1] = x_min[1]; - residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_c_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_c_dir(t_out) <= 0) { + t_max = t_out; + } } t = modifiedRegulaFalsi(residual_c_dir, 0., t_max, max_iters_falsi, tol, iters_used_falsi); + if (std::abs(residual_c_dir(t)) > tol) { + std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl; + } residual_c_dir.compute_new_x(x, t); residual.computeGradient(x, res, gradient, false, 1); res_s_done = false; @@ -757,13 +968,22 @@ namespace Opm if (res[0] < 0) { end_point[0] = x_max[0]; end_point[1] = x_min[1]; - residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_s_dir(t_out) >= 0) { + t_max = t_out; + } } else { end_point[0] = x_min[0]; end_point[1] = x_max[1]; - residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max); + residual_s_dir.setup(x, direction, end_point, x_min, x_max, t_max, t_out); + if (residual_s_dir(t_out) <= 0) { + t_max = t_out; + } } t = modifiedRegulaFalsi(residual_s_dir, 0., t_max, max_iters_falsi, tol, iters_used_falsi); + if (std::abs(residual_s_dir(t)) > tol) { + std::cout << "modifiedRegulaFalsi did not produce result under tolerance." << std::endl; + } residual_s_dir.compute_new_x(x, t); res_s_done = true; residual.computeGradient(x, res, gradient, true, 1); @@ -830,11 +1050,11 @@ namespace Opm } while (((max_s_change > tol_) || (max_c_change > tol_)) && ++num_iters < maxit_); if (max_s_change > tol_) { THROW("In solveMultiCell(), we did not converge after " - << num_iters << " iterations. Delta s = " << max_s_change); + << num_iters << " iterations. Delta s = " << max_s_change); } if (max_c_change > tol_) { THROW("In solveMultiCell(), we did not converge after " - << num_iters << " iterations. Delta c = " << max_c_change); + << num_iters << " iterations. Delta c = " << max_c_change); } std::cout << "Solved " << num_cells << " cell multicell problem in " << num_iters << " iterations." << std::endl; diff --git a/opm/polymer/TransportModelPolymer.hpp b/opm/polymer/TransportModelPolymer.hpp index 71306ede4..432ba5ebd 100644 --- a/opm/polymer/TransportModelPolymer.hpp +++ b/opm/polymer/TransportModelPolymer.hpp @@ -131,6 +131,7 @@ namespace Opm // Residual functions which are used in splitting method struct ResidualCDir; struct ResidualSDir; + struct ResidualDir; struct Residual; double fracFlow(double s, double c, int cell) const;