OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added profiling branch. Compute number of residual evaluations.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-06-11 14:44:21 +02:00
parent 9eb9ba372d
commit 02c5935865
3 changed files with 90 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

39
examples/polymer_reorder.cpp
View File

@ -28,7 +28,7 @@
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/utility/SimulatorTimer.hpp>
@ -38,7 +38,7 @@
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/polymer/IncompPropertiesDefaultPolymer.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
@ -53,7 +53,7 @@
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/ColumnExtract.hpp>
#include <opm/polymer/PolymerState.hpp>
#include <opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp>
@ -78,14 +78,15 @@
#include <iterator>
#include <vector>
#include <numeric>
#include <list>
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<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(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<Opm::TransportModelPolymer::Newton_Iter> 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<double> 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);

81
opm/polymer/TransportModelPolymer.cpp
View File

@ -24,7 +24,7 @@
#include <opm/core/utility/RootFinders.hpp>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <cmath>
#include <list>
// Choose error policy for scalar solves here.
typedef Opm::RegulaFalsi<Opm::WarnAndContinueOnError> 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];
}

20
opm/polymer/TransportModelPolymer.hpp
View File

@ -24,6 +24,7 @@
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <vector>
#include <list>
class UnstructuredGrid;
@ -82,6 +83,24 @@ namespace Opm
std::vector<double>& concentration,
std::vector<double>& 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<Newton_Iter> 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