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

merge. (included changes made in profiling branch in default branch).

Browse Source
This commit is contained in:
Xavier Raynaud
2012-06-14 09:38:45 +02:00
parent 502215832d 9df87db47f
commit 9ce282b069
5 changed files with 70 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
examples/Makefile.am
View File

@@ -6,7 +6,8 @@ LDFLAGS = $(BOOST_LDFLAGS)
LDADD = $(top_builddir)/libopmpolymer.la \
$(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB)
$(BOOST_SYSTEM_LIB) \
/usr/local/lib/libopmcore.la
noinst_PROGRAMS = \

35
examples/polymer_reorder.cpp
View File

@@ -79,14 +79,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;
@@ -116,6 +117,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();
}
@@ -478,15 +496,16 @@ main(int argc, char** argv)
std::vector<double> 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] = 20.0;
visc_mult_vals[1] = 1.0;
std::vector<double> c_vals_ads(3, -1e100);
c_vals_ads[0] = 0.0;
c_vals_ads[1] = 2.0;
c_vals_ads[2] = 8.0;
// Here we set up adsorption equal to zero.
std::vector<double> ads_vals(3, -1e100);
ads_vals[0] = 0.0;
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,
@@ -588,6 +607,8 @@ main(int argc, char** argv)
method = Opm::TransportModelPolymer::Bracketing;
} else if (method_string == "Newton") {
method = Opm::TransportModelPolymer::Newton;
} else if (method_string == "Gradient") {
method = Opm::TransportModelPolymer::Gradient;
} else {
THROW("Unknown method: " << method_string);
}
@@ -684,7 +705,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.
@@ -862,7 +883,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);

2
opm/polymer/PolymerProperties.cpp
View File

@@ -211,7 +211,7 @@ namespace Opm
double rk = 1 + (res_factor_ - 1)*c_ads/c_max_ads_;
eff_relperm_wat = relperm[0]/rk;
if (if_with_der) {
deff_relperm_wat_ds = (drelperm_ds[0]-drelperm_ds[3])/rk; //derivative with respect to sw
deff_relperm_wat_ds = (drelperm_ds[0]-drelperm_ds[2])/rk; //derivative with respect to sw
deff_relperm_wat_dc = dc_ads_dc*relperm[0]/(rk*rk*c_max_ads_);
} else {
deff_relperm_wat_ds = -1.0;

29
opm/polymer/TransportModelPolymer.cpp
View File

@@ -25,7 +25,7 @@
#include <opm/core/utility/miscUtilities.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;
@@ -51,9 +51,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;
@@ -200,6 +200,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);
@@ -239,6 +241,7 @@ namespace Opm
toWaterSat(saturation, saturation_);
concentration_ = &concentration[0];
cmax_ = &cmax[0];
res_counts.clear();
reorderAndTransport(grid_, darcyflux);
toBothSat(saturation_, saturation);
}
@@ -254,7 +257,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_;
@@ -263,7 +266,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,
@@ -285,6 +288,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_);
@@ -311,8 +315,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;
@@ -371,6 +375,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));
}
@@ -393,6 +399,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)
@@ -413,7 +420,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;
@@ -552,12 +559,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);
@@ -577,6 +586,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);
@@ -586,6 +596,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]));
}
}
@@ -678,7 +689,7 @@ namespace Opm
double res[2];
double mc;
double ff;
res_eq.computeResidual(x, res, mc, ff);
res_eq.computeResidual(x, res, mc, ff);
if (norm(res) <= tol_) {
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
fractionalflow_[cell] = ff;

19
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;
@@ -113,6 +114,24 @@ namespace Opm
const double* gravflux);
int solveGravityColumn(const std::vector<int>& cells);
// 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_;