diff --git a/examples/Utilities.cpp b/examples/Utilities.cpp index 3b5366458..11a965258 100644 --- a/examples/Utilities.cpp +++ b/examples/Utilities.cpp @@ -211,8 +211,7 @@ namespace Opm void writeVtkDataGeneralGrid(const UnstructuredGrid* grid, - const std::vector& pressure, - const std::vector& saturation, + const DataMap& data, std::ostream& os) { if (grid->dimensions != 3) { @@ -350,37 +349,32 @@ namespace Opm } { pm.clear(); - pm["Scalars"] = "saturation"; + 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"; - { - pm["Name"] = "pressure"; + 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 << pressure[c] << ' '; - if (c % num_per_line == num_per_line - 1 - || c == num_cells - 1) { - os << '\n'; - } - } - } - { - pm["Name"] = "saturation"; - 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 << saturation[2*c] << ' '; + os << field[stride*c] << ' '; if (c % num_per_line == num_per_line - 1 || c == num_cells - 1) { os << '\n'; diff --git a/examples/Utilities.hpp b/examples/Utilities.hpp index 5b69443be..3ff7c6976 100644 --- a/examples/Utilities.hpp +++ b/examples/Utilities.hpp @@ -106,6 +106,12 @@ namespace Opm 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_); @@ -204,9 +210,10 @@ namespace Opm + typedef std::map*> DataMap; + void writeVtkDataGeneralGrid(const UnstructuredGrid* grid, - const std::vector& pressure, - const std::vector& saturation, + const DataMap& data, std::ostream& os); diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index 7932b27e4..80f4df660 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -26,6 +26,7 @@ #include #include +#include #include #include #include @@ -60,13 +61,64 @@ +class AdHocProps : public Opm::IncompPropertiesBasic +{ +public: + AdHocProps(const Opm::parameter::ParameterGroup& param, int dim, int num_cells) + : Opm::IncompPropertiesBasic(param, dim, num_cells) + { + ASSERT(numPhases() == 2); + sw_.resize(3); + sw_[0] = 0.2; + sw_[1] = 0.7; + sw_[2] = 1.0; + krw_.resize(3); + krw_[0] = 0.0; + krw_[1] = 0.7; + krw_[2] = 1.0; + so_.resize(2); + so_[0] = 0.3; + so_[1] = 0.8; + kro_.resize(2); + kro_[0] = 0.0; + kro_[1] = 1.0; + } + virtual void relperm(const int n, + const double* s, + const int* /*cells*/, + double* kr, + double* dkrds) const + { + ASSERT(dkrds == 0); + for (int i = 0; i < n; ++i) { + kr[2*i] = krw(s[2*i]); + kr[2*i+1] = kro(s[2*i+1]); + } + } + + +private: + double krw(double s) const + { + return Opm::linearInterpolation(sw_, krw_, s); + } + + double kro(double s) const + { + return Opm::linearInterpolation(so_, kro_, s); + } + std::vector sw_; + std::vector krw_; + std::vector so_; + std::vector kro_; +}; class ReservoirState { public: - ReservoirState(const UnstructuredGrid* g, const int num_phases = 2) + ReservoirState(const UnstructuredGrid* g, const int num_phases = 2, const double init_sat = 0.0) : press_ (g->number_of_cells, 0.0), fpress_(g->number_of_faces, 0.0), flux_ (g->number_of_faces, 0.0), @@ -75,7 +127,8 @@ public: cmax_(g->number_of_cells, 0.0) { for (int cell = 0; cell < g->number_of_cells; ++cell) { - sat_[num_phases*cell + num_phases - 1] = 1.0; + sat_[num_phases*cell] = init_sat; + sat_[num_phases*cell + num_phases - 1] = 1.0 - init_sat; } } @@ -109,7 +162,11 @@ private: double polymerInflowAtTime(double time) { - return time >= 4.0*Opm::unit::day ? 1.0 : 0.0; + if (time >= 300.0*Opm::unit::day && time < 800.0*Opm::unit::day) { + return 5.0; + } else { + return 0.0; + } } @@ -121,13 +178,30 @@ void outputState(const UnstructuredGrid* grid, const int step, const std::string& output_dir) { + // Write data in VTK format. std::ostringstream vtkfilename; vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); if (!vtkfile) { THROW("Failed to open " << vtkfilename.str()); } - Opm::writeVtkDataGeneralGrid(grid, state.pressure(), state.saturation(), vtkfile); + Opm::DataMap dm; + dm["saturation"] = &state.saturation(); + dm["pressure"] = &state.pressure(); + dm["concentration"] = &state.concentration(); + Opm::writeVtkDataGeneralGrid(grid, dm, vtkfile); + + // Write data (not grid) in Matlab format + for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { + std::ostringstream fname; + fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat"; + std::ofstream file(fname.str().c_str()); + if (!file) { + THROW("Failed to open " << fname.str()); + } + const std::vector& d = *(it->second); + std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); + } } @@ -176,23 +250,34 @@ main(int argc, char** argv) const int nx = param.getDefault("nx", 100); const int ny = param.getDefault("ny", 100); const int nz = param.getDefault("nz", 1); - grid.reset(new Opm::Grid(nx, ny, nz)); + 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)); // Rock and fluid init. - props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - polydata.c_max_limit = param.getDefault("c_max_limit", 1.0); + // 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[0] = polydata.c_max_limit; + 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.c_vals_ads = polydata.c_vals_visc; - polydata.ads_vals.resize(2); - polydata.ads_vals[0] = 1.0; - polydata.ads_vals[1] = param.getDefault("c_max_ads", 0.0025); + // 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; } // Extra rock init. @@ -205,11 +290,15 @@ main(int argc, char** argv) // State-related and source-related variables init. std::vector totmob; - ReservoirState state(grid->c_grid(), props->numPhases()); + 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 // code expects a scalar sw, not both sw and so. std::vector reorder_sat(grid->c_grid()->number_of_cells); double flow_per_sec = 0.1*tot_porevol/Opm::unit::day; + if (param.has("injection_rate_per_day")) { + flow_per_sec = param.get("injection_rate_per_day")/Opm::unit::day; + } std::vector src (grid->c_grid()->number_of_cells, 0.0); src[0] = flow_per_sec; src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec; diff --git a/opm/polymer/polymermodel.cpp b/opm/polymer/polymermodel.cpp index e7f254f9a..180e24623 100644 --- a/opm/polymer/polymermodel.cpp +++ b/opm/polymer/polymermodel.cpp @@ -43,15 +43,13 @@ struct ParametersCRes double porosity; PolymerSolverData* psdata; int cell; - NonlinearSolverCtrl* ctrl; double s; const PolymerData* polydata; }; static struct ParametersSRes get_parameters_s(struct PolymerSolverData *d, int cell); -static struct ParametersCRes get_parameters_c(struct PolymerSolverData *d, int cell, - NonlinearSolverCtrl* ctrl); +static struct ParametersCRes get_parameters_c(struct PolymerSolverData *d, int cell); static double residual_s(double s, void *data); static double residual_c(double c, void *data); static double fluxfun_props(double s, @@ -124,12 +122,21 @@ init_solverdata(struct UnstructuredGrid *grid, } /* Solver for single-cell bvp calls root-finder in nlsolvers.c */ -void polymer_solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell) +void polymer_solvecell(void *data, int cell) { struct PolymerSolverData *d = (struct PolymerSolverData*) data; - struct ParametersCRes prm = get_parameters_c(d, cell, ctrl); - d->concentration[cell] = find_zero(residual_c, &prm, ctrl); + struct NonlinearSolverCtrl ctrl; + ctrl.maxiterations = 20; + ctrl.nltolerance = 1e-9; + ctrl.method = NonlinearSolverCtrl::REGULAFALSI; + ctrl.initialguess = 0.3*d->polydata->c_max_limit; + ctrl.min_bracket = 0.0; + ctrl.max_bracket = d->polydata->c_max_limit; + + struct ParametersCRes prm = get_parameters_c(d, cell); + + d->concentration[cell] = find_zero(residual_c, &prm, &ctrl); d->cmax[cell] = std::max(d->cmax[cell], d->concentration[cell]); d->saturation[cell] = prm.s; d->fractionalflow[cell] = fluxfun_props(d->saturation[cell], d->concentration[cell], cell, d->props, d->polydata); @@ -160,7 +167,14 @@ residual_c(double c, void *data) int cell = p->cell; struct ParametersSRes prm_s = get_parameters_s(p->psdata, cell); prm_s.c = c; - double s = find_zero(residual_s, &prm_s, p->ctrl); + struct NonlinearSolverCtrl ctrl; + ctrl.maxiterations = 20; + ctrl.nltolerance = 1e-9; + ctrl.method = NonlinearSolverCtrl::REGULAFALSI; + ctrl.initialguess = 0.5; + ctrl.min_bracket = 0.2; // TODO: Make this a proper s_min value. + ctrl.max_bracket = 1.0; + double s = find_zero(residual_s, &prm_s, &ctrl); p->s = s; double ff = fluxfun_props(s, c, p->cell, prm_s.props, p->polydata); double mc = compute_mc(c, prm_s.props, p->polydata); @@ -176,7 +190,7 @@ residual_c(double c, void *data) + rhor*((1.0 - porosity)/porosity)*(ads - ads0) + p->dtpv*(p->outflux*ff*mc + p->influx_polymer); #ifdef EXTRA_DEBUG_OUTPUT - std::cout << "res(c) = " << res << std::endl; + std::cout << "c = " << c << " s = " << s << " c-residual = " << res << std::endl; #endif return res; } @@ -226,7 +240,7 @@ get_parameters_s(struct PolymerSolverData *d, int cell) static struct ParametersCRes -get_parameters_c(struct PolymerSolverData *d, int cell, NonlinearSolverCtrl* ctrl) +get_parameters_c(struct PolymerSolverData *d, int cell) { int i; struct UnstructuredGrid *g = d->grid; @@ -241,11 +255,10 @@ get_parameters_c(struct PolymerSolverData *d, int cell, NonlinearSolverCtrl* ctr double src = d->source[cell]; p.influx = src > 0 ? -src : 0.0; p.influx_polymer = src > 0 ? -src*compute_mc(d->inflow_c, d->props, d->polydata) : 0.0; - p.outflux = d->source[cell] <= 0 ? -d->source[cell] : 0.0; + p.outflux = src <= 0 ? -src : 0.0; p.porosity = d->porosity[cell]; p.psdata = d; p.cell = cell; - p.ctrl = ctrl; p.s = -1e100; p.polydata = d->polydata; for (i=g->cell_facepos[cell]; icell_facepos[cell+1]; ++i) { @@ -272,8 +285,12 @@ get_parameters_c(struct PolymerSolverData *d, int cell, NonlinearSolverCtrl* ctr } } #ifdef EXTRA_DEBUG_OUTPUT - std::cout << "in: " << p.influx << " in(polymer): " << p.influx_polymer - << " out: " << p.outflux << std::endl; + std::cout << "Cell: " << cell + << " in: " << p.influx + << " in(polymer): " << p.influx_polymer + << " out: " << p.outflux << '\n' + << " c0: " << p.c0 + << " s0: " << p.s0 << std::endl; #endif return p; } diff --git a/opm/polymer/polymermodel.hpp b/opm/polymer/polymermodel.hpp index 82a165a65..9a1d1b6c0 100644 --- a/opm/polymer/polymermodel.hpp +++ b/opm/polymer/polymermodel.hpp @@ -51,11 +51,9 @@ struct PolymerSolverData { double *mc; }; -struct NonlinearSolverCtrl; - void -polymer_solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell); +polymer_solvecell(void *data, int cell); void destroy_solverdata(struct PolymerSolverData *d); diff --git a/opm/polymer/polymertransport.cpp b/opm/polymer/polymertransport.cpp index b7f2a67e4..db98d83d8 100644 --- a/opm/polymer/polymertransport.cpp +++ b/opm/polymer/polymertransport.cpp @@ -39,10 +39,6 @@ void polymertransport( inflow_c, saturation, concentration, cmax); - - struct NonlinearSolverCtrl ctrl; - - /* Compute sequence of single-cell problems */ sequence = (int*) malloc( grid->number_of_cells * sizeof *sequence); components = (int*) malloc(( grid->number_of_cells+1) * sizeof *components); @@ -50,13 +46,6 @@ void polymertransport( compute_sequence(grid, darcyflux, sequence, components, &ncomponents); assert(ncomponents == grid->number_of_cells); - - - ctrl.maxiterations = 20; - ctrl.nltolerance = 1e-9; - ctrl.method = NonlinearSolverCtrl::REGULAFALSI; - ctrl.initialguess = 0.5; - /* Assume all strong components are single-cell domains. */ for(i=0; inumber_of_cells; ++i) { @@ -68,7 +57,7 @@ void polymertransport( break; } #endif - polymer_solvecell(data, &ctrl, sequence[i]); + polymer_solvecell(data, sequence[i]); /*print("iterations:%d residual:%20.16f\n", * ctrl.iterations, ctrl.residual);*/ }