diff --git a/examples/Utilities.cpp b/examples/Utilities.cpp index 02671d4b0..10e5fa7d8 100644 --- a/examples/Utilities.cpp +++ b/examples/Utilities.cpp @@ -104,286 +104,40 @@ namespace Opm } - - void writeVtkDataAllCartesian(const std::tr1::array& dims, - const std::tr1::array& cell_size, - const std::vector& pressure, - const std::vector& saturation, - std::ostream& vtk_file) + void + compute_totmob_omega(const Opm::IncompPropertiesInterface& props, + const std::vector& s, + std::vector& totmob, + std::vector& 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 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 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 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_ << "\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 << "\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(num_pts); - pm["NumberOfCells"] = boost::lexical_cast(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 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 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(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 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(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& 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& sboth, std::vector& sw) { diff --git a/examples/Utilities.hpp b/examples/Utilities.hpp index 4df0abcc7..f550a4237 100644 --- a/examples/Utilities.hpp +++ b/examples/Utilities.hpp @@ -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& zcorn = deck.getFloatingPointValue("ZCORN"); - const std::vector& coord = deck.getFloatingPointValue("COORD"); - const std::vector& actnum = deck.getIntegerValue("ACTNUM"); - std::vector 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(g); - tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]); - - h_ = ifs_tpfa_construct(gg); - } - - ~PressureSolver() - { - ifs_tpfa_destroy(h_); - } - - template - void - solve(const UnstructuredGrid* g , - const ::std::vector& totmob, - const ::std::vector& src , - State& state ) - { - UnstructuredGrid* gg = const_cast(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 htrans_; - ::std::vector trans_ ; - ::std::vector gpress_; - - struct ifs_tpfa_data* h_; - }; - - - - void compute_porevolume(const UnstructuredGrid* g, const Opm::IncompPropertiesInterface& props, @@ -198,23 +71,11 @@ namespace Opm std::vector& totmob); - - - void writeVtkDataAllCartesian(const std::tr1::array& dims, - const std::tr1::array& cell_size, - const std::vector& pressure, - const std::vector& saturation, - std::ostream& vtk_file); - - - - typedef std::map*> DataMap; - - void writeVtkDataGeneralGrid(const UnstructuredGrid* grid, - const DataMap& data, - std::ostream& os); - - + void + compute_totmob_omega(const Opm::IncompPropertiesInterface& props, + const std::vector& s, + std::vector& totmob, + std::vector& omega); void toWaterSat(const std::vector& sboth, std::vector& sw); diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index b0ee466f0..eb546e6d7 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -22,22 +22,22 @@ #include "Utilities.hpp" -#include -#include +#include -#include +#include +#include +#include #include #include #include #include -#include #include -#include #include #include -#include +#include +// #include #include #include @@ -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; - } else { - return 0.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 grid; + boost::scoped_ptr grid; boost::scoped_ptr 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("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 global_cell(gc, gc + grid->c_grid()->number_of_cells); @@ -271,47 +285,75 @@ 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); - polydata.rhor = param.getDefault("rock_density", 1000.0); - polydata.dps = param.getDefault("dead_pore_space", 0.15); - polydata.c_vals_visc.resize(2); - polydata.c_vals_visc[0] = 0.0; - polydata.c_vals_visc[1] = 7.0; - polydata.visc_mult_vals.resize(2); - polydata.visc_mult_vals[0] = 1.0; - // polydata.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0); - polydata.visc_mult_vals[1] = 20.0; - polydata.c_vals_ads.resize(3); - polydata.c_vals_ads[0] = 0.0; - polydata.c_vals_ads[1] = 2.0; - polydata.c_vals_ads[2] = 8.0; - polydata.ads_vals.resize(3); - polydata.ads_vals[0] = 0.0; - // polydata.ads_vals[1] = param.getDefault("c_max_ads", 0.0025); - polydata.ads_vals[1] = 0.0015; - polydata.ads_vals[2] = 0.0025; } + + // 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); + polydata.rhor = param.getDefault("rock_density", 1000.0); + polydata.dps = param.getDefault("dead_pore_space", 0.15); + polydata.c_vals_visc.resize(2); + polydata.c_vals_visc[0] = 0.0; + polydata.c_vals_visc[1] = 7.0; + polydata.visc_mult_vals.resize(2); + polydata.visc_mult_vals[0] = 1.0; + // polydata.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0); + polydata.visc_mult_vals[1] = 20.0; + polydata.c_vals_ads.resize(3); + polydata.c_vals_ads[0] = 0.0; + polydata.c_vals_ads[1] = 2.0; + polydata.c_vals_ads[2] = 8.0; + polydata.ads_vals.resize(3); + polydata.ads_vals[0] = 0.0; + // 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 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 totmob; + std::vector 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); } - compute_totmob(*props, state.saturation(), totmob); + 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 diff --git a/opm/polymer/TransportModelPolymer.cpp b/opm/polymer/TransportModelPolymer.cpp index d861fa2ee..4e5ed521c 100644 --- a/opm/polymer/TransportModelPolymer.cpp +++ b/opm/polymer/TransportModelPolymer.cpp @@ -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 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 s0(num_cells); + std::vector c0(num_cells); + std::vector 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; } diff --git a/opm/polymer/TransportModelPolymer.hpp b/opm/polymer/TransportModelPolymer.hpp index 3cf88d843..71306ede4 100644 --- a/opm/polymer/TransportModelPolymer.hpp +++ b/opm/polymer/TransportModelPolymer.hpp @@ -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 smin_; + std::vector smax_; + double tol_; + double maxit_; const double* darcyflux_; // one flux per grid face const double* source_; // one source per cell