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.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-02-23 17:00:43 +01:00
commit f0fc7bf3c0
5 changed files with 206 additions and 475 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

296
examples/Utilities.cpp
View File

@ -104,286 +104,40 @@ namespace Opm
}
void writeVtkDataAllCartesian(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size,
const std::vector<double>& pressure,
const std::vector<double>& saturation,
std::ostream& vtk_file)
void
compute_totmob_omega(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega)
{
// Dimension is hardcoded in the prototype and the next two lines,
// but the rest is flexible (allows dimension == 2 or 3).
int dimension = 3;
int num_cells = dims[0]*dims[1]*dims[2];
ASSERT(dimension == 2 || dimension == 3);
ASSERT(num_cells = dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2]));
vtk_file << "# vtk DataFile Version 2.0\n";
vtk_file << "Structured Grid\n \n";
vtk_file << "ASCII \n";
vtk_file << "DATASET STRUCTURED_POINTS\n";
vtk_file << "DIMENSIONS "
<< dims[0] + 1 << " "
<< dims[1] + 1 << " ";
if (dimension == 3) {
vtk_file << dims[2] + 1;
} else {
vtk_file << 1;
int num_cells = props.numCells();
int num_phases = props.numPhases();
totmob.resize(num_cells);
omega.resize(num_cells);
ASSERT(int(s.size()) == num_cells*num_phases);
std::vector<int> cells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
cells[cell] = cell;
}
vtk_file << "\n";
vtk_file << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n";
vtk_file << "SPACING " << cell_size[0] << " " << cell_size[1];
if (dimension == 3) {
vtk_file << " " << cell_size[2];
} else {
vtk_file << " " << 0.0;
std::vector<double> kr(num_cells*num_phases);
props.relperm(num_cells, &s[0], &cells[0], &kr[0], 0);
const double* mu = props.viscosity();
for (int cell = 0; cell < num_cells; ++cell) {
totmob[cell] = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
totmob[cell] += kr[2*cell + phase]/mu[phase];
}
vtk_file << "\n";
vtk_file << "CELL_DATA " << num_cells << '\n';
vtk_file << "SCALARS pressure float" << '\n';
vtk_file << "LOOKUP_TABLE pressure_table " << '\n';
for (int i = 0; i < num_cells; ++i) {
vtk_file << pressure[i] << '\n';
}
vtk_file << "SCALARS saturation float" << '\n';
vtk_file << "LOOKUP_TABLE saturation_table " << '\n';
for (int i = 0; i < num_cells; ++i) {
double s = saturation[2*i];
if (s > 1e-10) {
vtk_file << s << '\n';
} else {
vtk_file << 0.0 << '\n';
const double* rho = props.density();
for (int cell = 0; cell < num_cells; ++cell) {
omega[cell] = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
omega[cell] += rho[phase]*(kr[2*cell + phase]/mu[phase])/totmob[cell];
}
}
}
typedef std::map<std::string, std::string> PMap;
struct Tag
{
Tag(const std::string& tag, const PMap& props, std::ostream& os)
: name_(tag), os_(os)
{
indent(os);
os << "<" << tag;
for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) {
os << " " << it->first << "=\"" << it->second << "\"";
}
os << ">\n";
++indent_;
}
Tag(const std::string& tag, std::ostream& os)
: name_(tag), os_(os)
{
indent(os);
os << "<" << tag << ">\n";
++indent_;
}
~Tag()
{
--indent_;
indent(os_);
os_ << "</" << name_ << ">\n";
}
static void indent(std::ostream& os)
{
for (int i = 0; i < indent_; ++i) {
os << " ";
}
}
private:
static int indent_;
std::string name_;
std::ostream& os_;
};
int Tag::indent_ = 0;
void writeVtkDataGeneralGrid(const UnstructuredGrid* grid,
const DataMap& data,
std::ostream& os)
{
if (grid->dimensions != 3) {
THROW("Vtk output for 3d grids only");
}
os.precision(12);
os << "<?xml version=\"1.0\"?>\n";
PMap pm;
pm["type"] = "UnstructuredGrid";
Tag vtkfiletag("VTKFile", pm, os);
Tag ugtag("UnstructuredGrid", os);
int num_pts = grid->number_of_nodes;
int num_cells = grid->number_of_cells;
pm.clear();
pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts);
pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells);
Tag piecetag("Piece", pm, os);
{
Tag pointstag("Points", os);
pm.clear();
pm["type"] = "Float64";
pm["Name"] = "Coordinates";
pm["NumberOfComponents"] = "3";
pm["format"] = "ascii";
Tag datag("DataArray", pm, os);
for (int i = 0; i < num_pts; ++i) {
Tag::indent(os);
os << grid->node_coordinates[3*i + 0] << ' '
<< grid->node_coordinates[3*i + 1] << ' '
<< grid->node_coordinates[3*i + 2] << '\n';
}
}
{
Tag cellstag("Cells", os);
pm.clear();
pm["type"] = "Int32";
pm["NumberOfComponents"] = "1";
pm["format"] = "ascii";
std::vector<int> cell_numpts;
cell_numpts.reserve(num_cells);
{
pm["Name"] = "connectivity";
Tag t("DataArray", pm, os);
int hf = 0;
for (int c = 0; c < num_cells; ++c) {
std::set<int> cell_pts;
for (; hf < grid->cell_facepos[c+1]; ++hf) {
int f = grid->cell_faces[hf];
const int* fnbeg = grid->face_nodes + grid->face_nodepos[f];
const int* fnend = grid->face_nodes + grid->face_nodepos[f+1];
cell_pts.insert(fnbeg, fnend);
}
cell_numpts.push_back(cell_pts.size());
Tag::indent(os);
std::copy(cell_pts.begin(), cell_pts.end(),
std::ostream_iterator<int>(os, " "));
os << '\n';
}
}
{
pm["Name"] = "offsets";
Tag t("DataArray", pm, os);
int offset = 0;
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
offset += cell_numpts[c];
os << offset << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
std::vector<int> cell_foffsets;
cell_foffsets.reserve(num_cells);
{
pm["Name"] = "faces";
Tag t("DataArray", pm, os);
const int* fp = grid->cell_facepos;
int offset = 0;
for (int c = 0; c < num_cells; ++c) {
Tag::indent(os);
os << fp[c+1] - fp[c] << '\n';
++offset;
for (int hf = fp[c]; hf < fp[c+1]; ++hf) {
int f = grid->cell_faces[hf];
const int* np = grid->face_nodepos;
int f_num_pts = np[f+1] - np[f];
Tag::indent(os);
os << f_num_pts << ' ';
++offset;
std::copy(grid->face_nodes + np[f],
grid->face_nodes + np[f+1],
std::ostream_iterator<int>(os, " "));
os << '\n';
offset += f_num_pts;
}
cell_foffsets.push_back(offset);
}
}
{
pm["Name"] = "faceoffsets";
Tag t("DataArray", pm, os);
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << cell_foffsets[c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
{
pm["type"] = "UInt8";
pm["Name"] = "types";
Tag t("DataArray", pm, os);
const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << "42 ";
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
}
{
pm.clear();
if (data.find("saturation") != data.end()) {
pm["Scalars"] = "saturation";
} else if (data.find("pressure") != data.end()) {
pm["Scalars"] = "pressure";
}
Tag celldatatag("CellData", pm, os);
pm.clear();
pm["type"] = "Int32";
pm["NumberOfComponents"] = "1";
pm["format"] = "ascii";
pm["type"] = "Float64";
for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) {
pm["Name"] = dit->first;
const std::vector<double>& field = *(dit->second);
// We always print only the first data item for every
// cell, using 'stride'.
// This is a hack to get water saturation nicely.
// \TODO: Extend to properly printing vector data.
const int stride = field.size()/grid->number_of_cells;
Tag ptag("DataArray", pm, os);
const int num_per_line = 5;
for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) {
Tag::indent(os);
}
os << field[stride*c] << ' ';
if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) {
os << '\n';
}
}
}
}
}
void toWaterSat(const std::vector<double>& sboth, std::vector<double>& sw)
{

149
examples/Utilities.hpp
View File

@ -59,133 +59,6 @@ namespace Opm
{
/// Concrete grid class constructing a
/// corner point grid from a deck,
/// or a cartesian grid.
class Grid
{
public:
Grid(const Opm::EclipseGridParser& deck)
{
// Extract data from deck.
const std::vector<double>& zcorn = deck.getFloatingPointValue("ZCORN");
const std::vector<double>& coord = deck.getFloatingPointValue("COORD");
const std::vector<int>& actnum = deck.getIntegerValue("ACTNUM");
std::vector<int> dims;
if (deck.hasField("DIMENS")) {
dims = deck.getIntegerValue("DIMENS");
} else if (deck.hasField("SPECGRID")) {
dims = deck.getSPECGRID().dimensions;
} else {
THROW("Deck must have either DIMENS or SPECGRID.");
}
// Collect in input struct for preprocessing.
struct grdecl grdecl;
grdecl.zcorn = &zcorn[0];
grdecl.coord = &coord[0];
grdecl.actnum = &actnum[0];
grdecl.dims[0] = dims[0];
grdecl.dims[1] = dims[1];
grdecl.dims[2] = dims[2];
// Process and compute.
ug_ = preprocess(&grdecl, 0.0);
compute_geometry(ug_);
}
Grid(int nx, int ny)
{
ug_ = create_cart_grid_2d(nx, ny);
}
Grid(int nx, int ny, int nz)
{
ug_ = create_cart_grid_3d(nx, ny, nz);
}
Grid(int nx, int ny, int nz,
double dx, double dy, double dz)
{
ug_ = create_hexa_grid_3d(nx, ny, nz, dx, dy, dz);
}
~Grid()
{
free_grid(ug_);
}
virtual const UnstructuredGrid* c_grid() const
{
return ug_;
}
private:
// Disable copying and assignment.
Grid(const Grid& other);
Grid& operator=(const Grid& other);
struct UnstructuredGrid* ug_;
};
class PressureSolver
{
public:
PressureSolver(const UnstructuredGrid* g,
const IncompPropertiesInterface& props)
: htrans_(g->cell_facepos[ g->number_of_cells ]),
trans_ (g->number_of_faces),
gpress_(g->cell_facepos[ g->number_of_cells ])
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]);
h_ = ifs_tpfa_construct(gg);
}
~PressureSolver()
{
ifs_tpfa_destroy(h_);
}
template <class State>
void
solve(const UnstructuredGrid* g ,
const ::std::vector<double>& totmob,
const ::std::vector<double>& src ,
State& state )
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
// No gravity
std::fill(gpress_.begin(), gpress_.end(), double(0.0));
ifs_tpfa_assemble(gg, &trans_[0], &src[0], &gpress_[0], h_);
using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
linsolve.solve(h_->A, h_->b, h_->x);
ifs_tpfa_press_flux(gg, &trans_[0], h_,
&state.pressure()[0],
&state.faceflux()[0]);
}
private:
::std::vector<double> htrans_;
::std::vector<double> trans_ ;
::std::vector<double> gpress_;
struct ifs_tpfa_data* h_;
};
void
compute_porevolume(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props,
@ -198,23 +71,11 @@ namespace Opm
std::vector<double>& totmob);
void writeVtkDataAllCartesian(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size,
const std::vector<double>& pressure,
const std::vector<double>& saturation,
std::ostream& vtk_file);
typedef std::map<std::string, const std::vector<double>*> DataMap;
void writeVtkDataGeneralGrid(const UnstructuredGrid* grid,
const DataMap& data,
std::ostream& os);
void
compute_totmob_omega(const Opm::IncompPropertiesInterface& props,
const std::vector<double>& s,
std::vector<double>& totmob,
std::vector<double>& omega);
void toWaterSat(const std::vector<double>& sboth, std::vector<double>& sw);

95
examples/polymer_reorder.cpp
View File

@ -22,22 +22,22 @@
#include "Utilities.hpp"
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/utility/cart_grid.h>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/utility/writeVtkData.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/cpgpreprocess/cgridinterface.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/linalg/LinearSolverUmfpack.hpp>
// #include <opm/core/linalg/LinearSolverIstl.hpp>
#include <opm/polymer/polymertransport.hpp>
#include <opm/polymer/polymermodel.hpp>
@ -177,15 +177,28 @@ private:
double polymerInflowAtTime(double time)
class PolymerInflow
{
if (time >= 300.0*Opm::unit::day && time < 800.0*Opm::unit::day) {
return 5.0;
public:
PolymerInflow(const double starttime,
const double endtime,
const double amount)
: stime_(starttime), etime_(endtime), amount_(amount)
{
}
double operator()(double time)
{
if (time >= stime_ && time < etime_) {
return amount_;
} else {
return 0.0;
}
}
private:
double stime_;
double etime_;
double amount_;
};
@ -207,7 +220,7 @@ void outputState(const UnstructuredGrid* grid,
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
dm["concentration"] = &state.concentration();
Opm::writeVtkDataGeneralGrid(grid, dm, vtkfile);
Opm::writeVtkData(grid, dm, vtkfile);
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
@ -250,15 +263,16 @@ main(int argc, char** argv)
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
boost::scoped_ptr<Opm::Grid> grid;
boost::scoped_ptr<Opm::GridManager> grid;
boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
Opm::PolymerData polydata;
if (use_deck) {
THROW("We do not yet read polymer keywords from deck.");
std::cout << "**** Warning: We do not yet read polymer keywords from deck. "
"Polymer behaviour is hardcoded." << std::endl;
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename);
// Grid init
grid.reset(new Opm::Grid(deck));
grid.reset(new Opm::GridManager(deck));
// Rock and fluid init
const int* gc = grid->c_grid()->global_cell;
std::vector<int> global_cell(gc, gc + grid->c_grid()->number_of_cells);
@ -271,10 +285,12 @@ main(int argc, char** argv)
const int dx = param.getDefault("dx", 1.0);
const int dy = param.getDefault("dy", 1.0);
const int dz = param.getDefault("dz", 1.0);
grid.reset(new Opm::Grid(nx, ny, nz, dx, dy, dz));
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 AdHocProps(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
}
// Setting polydata defaults to mimic a simple example case.
polydata.c_max_limit = param.getDefault("c_max_limit", 5.0);
polydata.omega = param.getDefault("omega", 1.0);
@ -296,22 +312,48 @@ main(int argc, char** argv)
// polydata.ads_vals[1] = param.getDefault("c_max_ads", 0.0025);
polydata.ads_vals[1] = 0.0015;
polydata.ads_vals[2] = 0.0025;
}
double poly_start = param.getDefault("poly_start_days", 300.0)*Opm::unit::day;
double poly_end = param.getDefault("poly_end_days", 800.0)*Opm::unit::day;
double poly_amount = param.getDefault("poly_amount", 5.0);
PolymerInflow poly_inflow(poly_start, poly_end, poly_amount);
bool new_code = param.getDefault("new_code", false);
int method = param.getDefault("method", 1);
// Extra rock init.
std::vector<double> porevol;
compute_porevolume(grid->c_grid(), *props, porevol);
double tot_porevol = std::accumulate(porevol.begin(), porevol.end(), 0.0);
// Solvers init.
Opm::PressureSolver psolver(grid->c_grid(), *props);
// Gravity init.
double gravity[3] = { 0.0 };
double g = param.getDefault("gravity", 0.0);
bool use_gravity = g != 0.0;
if (use_gravity) {
gravity[grid->c_grid()->dimensions - 1] = g;
if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
}
Opm::TransportModelPolymer tmodel(*grid->c_grid(), props->porosity(), &porevol[0], *props, polydata, method);
// Solvers init.
Opm::LinearSolverUmfpack linsolver;
// Opm::LinearSolverIstl linsolver(param);
Opm::IncompTpfa psolver(*grid->c_grid(),
props->permeability(),
use_gravity ? gravity : 0,
linsolver);
const int method = param.getDefault("method", 1);
const double nltol = param.getDefault("nl_tolerance", 1e-9);
const int maxit = param.getDefault("nl_maxiter", 30);
Opm::TransportModelPolymer tmodel(*grid->c_grid(), props->porosity(), &porevol[0], *props, polydata,
method, nltol, maxit);
// State-related and source-related variables init.
std::vector<double> totmob;
std::vector<double> omega; // Empty dummy unless/until we include gravity here.
double init_sat = param.getDefault("init_sat", 0.0);
ReservoirState state(grid->c_grid(), props->numPhases(), init_sat);
// We need a separate reorder_sat, because the reorder
@ -362,15 +404,24 @@ main(int argc, char** argv)
outputState(grid->c_grid(), state, pstep, output_dir);
}
if (use_gravity) {
compute_totmob_omega(*props, state.saturation(), totmob, omega);
} else {
compute_totmob(*props, state.saturation(), totmob);
}
pressure_timer.start();
psolver.solve(grid->c_grid(), totmob, src, state);
psolver.solve(totmob, omega, src, state.pressure(), state.faceflux());
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
double inflow_c = polymerInflowAtTime(current_time);
const double inflowc0 = poly_inflow(current_time + 1e-5*stepsize);
const double inflowc1 = poly_inflow(current_time + (1.0 - 1e-5)*stepsize);
if (inflowc0 != inflowc1) {
std::cout << "**** Warning: polymer inflow rate changes during timestep. Using rate near start of step.";
}
const double inflow_c = inflowc0;
Opm::toWaterSat(state.saturation(), reorder_sat);
// We must treat reorder_src here,
// if we are to handle anything but simple water

82
opm/polymer/TransportModelPolymer.cpp
View File

@ -34,12 +34,16 @@ namespace Opm
const double* porevolume,
const IncompPropertiesInterface& props,
const PolymerData& polyprops,
int method)
const int method,
const double tol,
const int maxit)
: grid_(grid),
porosity_(porosity),
porevolume_(porevolume),
props_(props),
polyprops_(polyprops),
tol_(tol),
maxit_(maxit),
darcyflux_(0),
source_(0),
dt_(0.0),
@ -55,6 +59,16 @@ namespace Opm
THROW("Property object must have 2 phases");
}
visc_ = props.viscosity();
// Set up smin_ and smax_
int num_cells = props.numCells();
smin_.resize(props.numPhases()*num_cells);
smax_.resize(props.numPhases()*num_cells);
std::vector<int> cells(num_cells);
for (int i = 0; i < num_cells; ++i) {
cells[i] = i;
}
props.satRange(props.numCells(), &cells[0], &smin_[0], &smax_[0]);
}
@ -180,13 +194,9 @@ namespace Opm
double operator()(double c) const
{
ResidualS res_s(tm, cell, s0, influx, outflux, dtpv, c);
const double a = 0.2; // TODO: Make this a proper s_min value.
const double b = 1.0;
const int maxit = 20;
const double tol = 1e-9;
int iters_used;
// Solve for s first.
s = modifiedRegulaFalsi(res_s, a, b, maxit, tol, iters_used);
s = modifiedRegulaFalsi(res_s, tm.smin_[2*cell], tm.smax_[2*cell], tm.maxit_, tm.tol_, iters_used);
double ff = tm.fracFlow(s, c, cell);
double mc = tm.computeMc(c);
double dps = tm.polyprops_.dps;
@ -642,10 +652,8 @@ namespace Opm
ResidualC res(*this, cell);
const double a = 0.0;
const double b = polyprops_.c_max_limit;
const int maxit = 20;
const double tol = 1e-9;
int iters_used;
concentration_[cell] = modifiedRegulaFalsi(res, a, b, maxit, tol, iters_used);
concentration_[cell] = modifiedRegulaFalsi(res, a, b, maxit_, tol_, iters_used);
cmax_[cell] = std::max(cmax_[cell], concentration_[cell]);
saturation_[cell] = res.lastSaturation();
fractionalflow_[cell] = fracFlow(saturation_[cell], concentration_[cell], cell);
@ -776,10 +784,60 @@ namespace Opm
}
}
void TransportModelPolymer::solveMultiCell(const int num_cells, const int* /*cells*/)
void TransportModelPolymer::solveMultiCell(const int num_cells, const int* cells)
{
THROW("TransportModelPolymer::solveMultiCell() not yet implemented, "
"got a component of size " << num_cells);
double max_s_change = 0.0;
double max_c_change = 0.0;
int num_iters = 0;
// Must store state variables before we start.
std::vector<double> s0(num_cells);
std::vector<double> c0(num_cells);
std::vector<double> cmax0(num_cells);
// Must set initial fractional flows etc. before we start.
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
fractionalflow_[cell] = fracFlow(saturation_[cell], concentration_[cell], cell);
mc_[cell] = computeMc(concentration_[cell]);
s0[i] = saturation_[cell];
c0[i] = concentration_[cell];
cmax0[i] = cmax_[i];
}
do {
int max_s_change_cell = -1;
int max_c_change_cell = -1;
max_s_change = 0.0;
max_c_change = 0.0;
for (int i = 0; i < num_cells; ++i) {
const int cell = cells[i];
const double old_s = saturation_[cell];
const double old_c = concentration_[cell];
saturation_[cell] = s0[i];
concentration_[cell] = c0[i];
cmax_[cell] = cmax0[i];
solveSingleCell(cell);
// std::cout << "cell = " << cell << " delta s = " << saturation_[cell] - old_s << std::endl;
if (max_s_change < std::fabs(saturation_[cell] - old_s)) {
max_s_change_cell = cell;
}
if (max_c_change < std::fabs(concentration_[cell] - old_c)) {
max_c_change_cell = cell;
}
max_s_change = std::max(max_s_change, std::fabs(saturation_[cell] - old_s));
max_c_change = std::max(max_c_change, std::fabs(concentration_[cell] - old_c));
}
// std::cout << "Iter = " << num_iters << " max_s_change = " << max_s_change
// << " in cell " << max_change_cell << std::endl;
} 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);
}
if (max_c_change > tol_) {
THROW("In solveMultiCell(), we did not converge after "
<< num_iters << " iterations. Delta c = " << max_c_change);
}
std::cout << "Solved " << num_cells << " cell multicell problem in "
<< num_iters << " iterations." << std::endl;
}

9
opm/polymer/TransportModelPolymer.hpp
View File

@ -75,12 +75,15 @@ namespace Opm
class TransportModelPolymer : public TransportModelInterface
{
public:
/// \TODO document me, especially method.
TransportModelPolymer(const UnstructuredGrid& grid,
const double* porosity,
const double* porevolume,
const IncompPropertiesInterface& props,
const PolymerData& polyprops,
int method);
const int method,
const double tol,
const int maxit);
/// Solve transport eqn with implicit Euler scheme, reordered.
/// \TODO Now saturation is expected to be one sw value per cell,
@ -105,6 +108,10 @@ namespace Opm
const double* porevolume_; // one volume per cell
const IncompPropertiesInterface& props_;
const PolymerData& polyprops_;
std::vector<double> smin_;
std::vector<double> smax_;
double tol_;
double maxit_;
const double* darcyflux_; // one flux per grid face
const double* source_; // one source per cell