Merge from upstream.

examples/Utilities.cpp

@@ -211,8 +211,7 @@ namespace Opm
 
 
     void writeVtkDataGeneralGrid(const UnstructuredGrid* grid,
-				 const std::vector<double>& pressure,
+				 const DataMap& data,
-				 const std::vector<double>& saturation,
 				 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";
-	    {
+	    for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) {
-		pm["Name"] = "pressure";
+		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 << pressure[c] << ' ';
+		    os << field[stride*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] << ' ';
 		    if (c % num_per_line == num_per_line - 1
 			|| c == num_cells - 1) {
 			os << '\n';

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<std::string, const std::vector<double>*> DataMap;
+
     void writeVtkDataGeneralGrid(const UnstructuredGrid* grid,
-				 const std::vector<double>& pressure,
+				 const DataMap& data,
-				 const std::vector<double>& saturation,
 				 std::ostream& os);
 
 

examples/polymer_reorder.cpp

@@ -26,6 +26,7 @@
 #include <opm/core/pressure/tpfa/trans_tpfa.h>
 
 #include <opm/core/utility/cart_grid.h>
+#include <opm/core/utility/linearInterpolation.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <opm/core/utility/StopWatch.hpp>
 #include <opm/core/utility/Units.hpp>
@@ -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<double> sw_;
+    std::vector<double> krw_;
+    std::vector<double> so_;
+    std::vector<double> 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<double>& d = *(it->second);
+	std::copy(d.begin(), d.end(), std::ostream_iterator<double>(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));
+	// 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 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.visc_mult_vals[1] = param.getDefault("c_max_viscmult", 30.0);
-	polydata.c_vals_ads = polydata.c_vals_visc;
+	polydata.visc_mult_vals[1] = 20.0;
-	polydata.ads_vals.resize(2);
+	polydata.c_vals_ads.resize(3);
-	polydata.ads_vals[0] = 1.0;
+	polydata.c_vals_ads[0] = 0.0;
-	polydata.ads_vals[1] = param.getDefault("c_max_ads", 0.0025);
+	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<double> 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<double> 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<double>("injection_rate_per_day")/Opm::unit::day;
+    }
     std::vector<double> 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;

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,
+static struct ParametersCRes get_parameters_c(struct PolymerSolverData *d, int cell);
-					      NonlinearSolverCtrl* ctrl);
 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]; i<g->cell_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
+    std::cout << "Cell: " << cell
-	      << "    out: " << p.outflux << std::endl;
+	      << "    in: " << p.influx
+	      << "    in(polymer): " << p.influx_polymer
+	      << "    out: " << p.outflux << '\n'
+	      << "           c0: " << p.c0
+	      << "    s0: " << p.s0 << std::endl;
 #endif
     return p;
 }

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);

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; i<grid->number_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);*/
     }