From 02c5935865d9df6d9a2c348a534570cc6d4598ec Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Mon, 11 Jun 2012 14:44:21 +0200 Subject: [PATCH] Added profiling branch. Compute number of residual evaluations. --- examples/polymer_reorder.cpp | 39 +++++++++---- opm/polymer/TransportModelPolymer.cpp | 81 ++++++++++++++------------- opm/polymer/TransportModelPolymer.hpp | 20 +++++++ 3 files changed, 90 insertions(+), 50 deletions(-) diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index b6cb3e48e..2b59791a3 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -28,7 +28,7 @@ #include #include #include -#include +#include #include #include #include @@ -38,7 +38,7 @@ #include #include -#include +#include #include #include @@ -53,7 +53,7 @@ #include #include -#include +#include #include #include @@ -78,14 +78,15 @@ #include #include #include - +#include static void outputState(const UnstructuredGrid& grid, const Opm::PolymerState& state, const int step, - const std::string& output_dir) + const std::string& output_dir, + const Opm::TransportModelPolymer& reorder_model) { // Write data in VTK format. std::ostringstream vtkfilename; @@ -115,6 +116,23 @@ static void outputState(const UnstructuredGrid& grid, const std::vector& d = *(it->second); std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); } + + std::ostringstream fname; + fname << output_dir << "/" << "residualcounts" << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + + typedef std::list ListRes; + + const ListRes& res_counts = reorder_model.res_counts; + for (ListRes::const_iterator it = res_counts.begin(); it != res_counts.end(); ++it) { + file << it->res_s << "," << it->cell << "," << std::setprecision(15) << it->s << "," << std::setprecision(15) << it->c << "\n"; + } + file.close(); + + } @@ -392,7 +410,7 @@ main(int argc, char** argv) grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz)); // Rock and fluid init. // props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - props.reset(new Opm::IncompPropertiesDefaultPolymer(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); + props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); // Wells init. wells.reset(new Opm::WellsManager()); // Timer init. @@ -442,9 +460,8 @@ main(int argc, char** argv) c_vals_ads[2] = 8.0; std::vector ads_vals(3, -1e100); ads_vals[0] = 0.0; - // polyprop.ads_vals[1] = param.getDefault("c_max_ads", 0.0025); - ads_vals[1] = 0.0015; - ads_vals[2] = 0.0025; + 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, @@ -646,7 +663,7 @@ main(int argc, char** argv) // Report timestep and (optionally) write state to disk. simtimer.report(std::cout); if (output && (simtimer.currentStepNum() % output_interval == 0)) { - outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir); + outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir, reorder_model); } // Solve pressure. @@ -860,7 +877,7 @@ main(int argc, char** argv) << "\n Transport time: " << ttime << std::endl; if (output) { - outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir); + outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir, reorder_model); outputWaterCut(watercut, output_dir); if (wells->c_wells()) { outputWellReport(wellreport, output_dir); diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp index f28cf8439..8413ddd23 100644 --- a/opm/polymer/TransportModelPolymer.cpp +++ b/opm/polymer/TransportModelPolymer.cpp @@ -24,7 +24,7 @@ #include #include #include - +#include // Choose error policy for scalar solves here. typedef Opm::RegulaFalsi RootFinder; @@ -50,9 +50,9 @@ public: double ads0; GradientMethod gradient_method; - const TransportModelPolymer& tm; + TransportModelPolymer& tm; - ResidualEquation(const TransportModelPolymer& tmodel, int cell_index); + ResidualEquation(TransportModelPolymer& tmodel, int cell_index); void computeResidual(const double* x, double* res) const; void computeResidual(const double* x, double* res, double& mc, double& ff) const; double computeResidualS(const double* x) const; @@ -126,7 +126,7 @@ namespace CurveInSCPlane(); void setup(const double* x, const double* direction, const double* end_point, const double* x_min, - const double* x_max, const bool adjust_dir, + const double* x_max, const bool adjust_dir, const double tol, double& t_max_out, double& t_out_out); void computeXOfT(double*, const double) const; @@ -202,6 +202,8 @@ namespace Opm } visc_ = props.viscosity(); + res_counts.clear(); + // Set up smin_ and smax_ int num_cells = props.numCells(); smin_.resize(props.numPhases()*num_cells); @@ -233,6 +235,7 @@ namespace Opm saturation_ = saturation; concentration_ = concentration; cmax_ = cmax; + res_counts.clear(); reorderAndTransport(grid_, darcyflux); } @@ -247,7 +250,7 @@ namespace Opm // Influxes are negative, outfluxes positive. struct TransportModelPolymer::ResidualS { - const TransportModelPolymer& tm_; + TransportModelPolymer& tm_; const int cell_; const double s0_; const double cmax0_; @@ -256,7 +259,7 @@ namespace Opm const double comp_term_; // q - sum_j v_ij const double dtpv_; // dt/pv(i) const double c_; - explicit ResidualS(const TransportModelPolymer& tmodel, + explicit ResidualS(TransportModelPolymer& tmodel, const int cell, const double s0, const double cmax0, @@ -278,6 +281,7 @@ namespace Opm } double operator()(double s) const { + tm_.res_counts.push_back(Newton_Iter(true, cell_, s, c_)); double ff; tm_.fracFlow(s, c_, cmax0_, cell_, ff); return s - s0_ + dtpv_*(outflux_*ff + influx_ + s*comp_term_); @@ -304,8 +308,8 @@ namespace Opm double porosity; double dtpv; // dt/pv(i) mutable double s; // Mutable in order to change it with every operator() call to be the last computed s value. - const TransportModelPolymer& tm; - explicit ResidualC(const TransportModelPolymer& tmodel, int cell_index) + TransportModelPolymer& tm; + explicit ResidualC(TransportModelPolymer& tmodel, int cell_index) : tm(tmodel) { cell = cell_index; @@ -364,6 +368,8 @@ namespace Opm + rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0) + dtpv*(outflux*ff*mc + influx_polymer) + dtpv*(s_arg*c_arg*(1.0 - dps) - rhor*c_ads)*comp_term; + tm.res_counts.push_back(Newton_Iter(true, cell, s_arg, c_arg)); + tm.res_counts.push_back(Newton_Iter(false, cell, s_arg, c_arg)); } @@ -386,6 +392,7 @@ namespace Opm tm.polyprops_.adsorption(c0, cmax0, c_ads0); double c_ads; tm.polyprops_.adsorption(c, cmax0, c_ads); + tm.res_counts.push_back(Newton_Iter(false, cell, s, c)); double res = (1 - dps)*s*c - (1 - dps)*s0*c0 + rhor*((1.0 - porosity)/porosity)*(c_ads - c_ads0) + dtpv*(outflux*ff*mc + influx_polymer) @@ -406,7 +413,7 @@ namespace Opm // ResidualEquation gathers parameters to construct the residual, computes its // value and the values of its derivatives. - TransportModelPolymer::ResidualEquation::ResidualEquation(const TransportModelPolymer& tmodel, int cell_index) + TransportModelPolymer::ResidualEquation::ResidualEquation(TransportModelPolymer& tmodel, int cell_index) : tm(tmodel) { gradient_method = Analytic; @@ -545,12 +552,14 @@ namespace Opm } if (if_res_s) { res[0] = s - s0 + dtpv*(outflux*ff + influx + s*comp_term); + tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1])); } if (if_res_c) { res[1] = (1 - dps)*s*c - (1 - dps)*s0*c0 + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + dtpv*(outflux*ff*mc + influx_polymer) + dtpv*(s*c*(1.0 - dps) - rhor*ads)*comp_term; + tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1])); } if (if_dres_s_dsdc) { dres_s_dsdc[0] = 1 + dtpv*(outflux*dff_dsdc[0] + comp_term); @@ -570,6 +579,7 @@ namespace Opm tm.fracFlow(s, c, cmax0, cell, ff); if (if_res_s) { res[0] = s - s0 + dtpv*(outflux*ff + influx + s*comp_term); + tm.res_counts.push_back(Newton_Iter(true, cell, x[0], x[1])); } if (if_res_c) { tm.computeMc(c, mc); @@ -579,6 +589,7 @@ namespace Opm + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + dtpv*(outflux*ff*mc + influx_polymer) + dtpv*(s*c*(1.0 - dps) - rhor*ads)*comp_term; + tm.res_counts.push_back(Newton_Iter(false, cell, x[0], x[1])); } } @@ -658,13 +669,6 @@ namespace Opm // curve. In these cases, we can use a robust 1d solver. void TransportModelPolymer::solveSingleCellNewton(int cell) { - // the tolerance for 1d solver is set as a function of the residual, because if we are far - // from the solution we do not need a very accurate 1d solver (recall that the 1d solver - // solves for only one of the two residuals) - // The tolerance falsi_tol is improved by - // (reduc_factor_falsi_tol * "previous residual") at each step - double falsi_tol; - const double reduc_factor_falsi_tol = 1e-2; int iters_used_falsi = 0; const int max_iters_split = maxit_; int iters_used_split = 0; @@ -683,7 +687,6 @@ namespace Opm return; } - falsi_tol = std::max(reduc_factor_falsi_tol*norm(res), tol_); // 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 }; @@ -706,7 +709,13 @@ namespace Opm while ((norm(res) > tol_) && (iters_used_split < max_iters_split)) { // We first try a Newton step if (counter_drop_newton == 0 && solveNewtonStep(x, res_eq, res, x_new)) { - res_eq.computeResidual(x_new, res_new, mc, ff); + if (check_interval(x_new, x_min, x_max)) { + // for testing + res_eq.computeResidual(x_new, res_new, mc, ff); + } else { + res_eq.computeResidual(x_new, res_new, mc, ff); + } + unsuccessfull_newton_step = false; not_so_successfull_newton_step = false; if (norm(res_new) > norm(res) || x_new[0] < x_min[0] || x_new[1] < x_min[1] || x_new[0] > x_max[0] || x_new[1] > x_max[1]) { @@ -714,7 +723,7 @@ namespace Opm } else { x[0] = x_new[0]; x[1] = x_new[1]; - if (norm(res_new) > 1e-1*norm(res) && norm(res_new) < 1e1*tol_) { + if (norm(res_new) > 1e-1*norm(res)) { // We are close to the solution and Newton does not perform well. // Then, we drop Newton for a given number of iterations. not_so_successfull_newton_step = true; @@ -722,9 +731,6 @@ namespace Opm } res[0] = res_new[0]; res[1] = res_new[1]; - if (check_interval(x, x_min, x_max)) { - res_eq.computeResidual(x, res, mc, ff); - } iters_used_split += 1; } } else { @@ -748,13 +754,12 @@ namespace Opm // We start with the zero curve of the s and r residual we are closest to. bool adjust_dir = true; if (std::abs(res[0]) < std::abs(res[1])) { - falsi_tol = std::max(reduc_factor_falsi_tol*std::abs(res[0]), tol_); - if (res[0] < -falsi_tol) { + if (res[0] < -tol_) { direction[0] = x_max[0] - x[0]; direction[1] = x_min[1] - x[1]; adjust_dir = true; if_res_s = true; - } else if (res[0] > falsi_tol) { + } else if (res[0] > tol_) { direction[0] = x_min[0] - x[0]; direction[1] = x_max[1] - x[1]; adjust_dir = true; @@ -763,17 +768,16 @@ namespace Opm res_eq.computeGradientResS(x, res, gradient); direction[0] = -gradient[1]; direction[1] = gradient[0]; - adjust_dir = false; + adjust_dir = true; if_res_s = false; } } else { - falsi_tol = std::max(reduc_factor_falsi_tol*std::abs(res[1]), tol_); - if (res[1] < -falsi_tol) { + if (res[1] < -tol_) { direction[0] = x_max[0] - x[0]; direction[1] = x_max[1] - x[1]; adjust_dir = true; if_res_s = false; - } else if (res[1] > falsi_tol) { + } else if (res[1] > tol_) { direction[0] = x_min[0] - x[0]; direction[1] = x_min[1] - x[1]; adjust_dir = true; @@ -782,7 +786,7 @@ namespace Opm res_eq.computeGradientResC(x, res, gradient); direction[0] = -gradient[1]; direction[1] = gradient[0]; - adjust_dir = false; + adjust_dir = true; if_res_s = true; } } @@ -790,42 +794,41 @@ namespace Opm if (res[0] < 0) { end_point[0] = x_max[0]; end_point[1] = x_min[1]; - res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, t_max, t_out); + res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out); if (res_s_on_curve(t_out) >= 0) { t_max = t_out; } } else { end_point[0] = x_min[0]; end_point[1] = x_max[1]; - res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, t_max, t_out); + res_s_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out); if (res_s_on_curve(t_out) <= 0) { t_max = t_out; } } // Note: In some experiments modifiedRegularFalsi does not yield a result under the given tolerance. - t = RootFinder::solve(res_s_on_curve, 0., t_max, maxit_, falsi_tol, iters_used_falsi); + t = RootFinder::solve(res_s_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi); res_s_on_curve.curve.computeXOfT(x, t); } else { if (res[1] < 0) { end_point[0] = x_max[0]; end_point[1] = x_max[1]; - res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, t_max, t_out); + res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out); if (res_c_on_curve(t_out) >= 0) { t_max = t_out; } } else { end_point[0] = x_min[0]; end_point[1] = x_min[1]; - res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, t_max, t_out); + res_c_on_curve.curve.setup(x, direction, end_point, x_min, x_max, adjust_dir, tol_, t_max, t_out); if (res_c_on_curve(t_out) <= 0) { t_max = t_out; } } - t = RootFinder::solve(res_c_on_curve, 0., t_max, maxit_, falsi_tol, iters_used_falsi); + t = RootFinder::solve(res_c_on_curve, 0., t_max, maxit_, tol_, iters_used_falsi); res_c_on_curve.curve.computeXOfT(x, t); } - check_interval(x, x_min, x_max); res_eq.computeResidual(x, res, mc, ff); iters_used_split += 1; } @@ -1272,7 +1275,7 @@ namespace // rectangle. x_out=(s_out, c_out) denotes the values of s and c at that point. void CurveInSCPlane::setup(const double* x, const double* direction, const double* end_point, const double* x_min, - const double* x_max, const bool adjust_dir, + const double* x_max, const bool adjust_dir, const double tol, double& t_max_out, double& t_out_out) { x_[0] = x[0]; @@ -1286,7 +1289,7 @@ namespace end_point_[0] = end_point[0]; end_point_[1] = end_point[1]; const double size_direction = std::abs(direction_[0]) + std::abs(direction_[1]); - if (size_direction == 0) { + if (size_direction < tol) { direction_[0] = end_point_[0]-x_[0]; direction_[1] = end_point_[1]-x_[1]; } diff --git a/opm/polymer/TransportModelPolymer.hpp b/opm/polymer/TransportModelPolymer.hpp index 4386c4413..93e61dc48 100644 --- a/opm/polymer/TransportModelPolymer.hpp +++ b/opm/polymer/TransportModelPolymer.hpp @@ -24,6 +24,7 @@ #include #include #include +#include class UnstructuredGrid; @@ -82,6 +83,24 @@ namespace Opm std::vector& concentration, std::vector& cmax); + // for testing + class Newton_Iter { + public: + bool res_s; + int cell; + double s; + double c; + + Newton_Iter(bool res_s_val, int cell_val, double s_val, double c_val) { + res_s = res_s_val; + cell = cell_val; + s = s_val; + c = c_val; + } + }; + + std::list res_counts; + private: const UnstructuredGrid& grid_; @@ -129,6 +148,7 @@ namespace Opm void computeMc(double c, double& mc) const; void computeMcWithDer(double c, double& mc, double& dmc_dc) const; void mobility(double s, double c, int cell, double* mob) const; + }; } // namespace Opm