Initial work on supporting polymer transport. Work in progress.

examples/polymer_reorder.cpp

@@ -38,6 +38,7 @@
 #include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
 
 #include <opm/polymer/polymertransport.hpp>
+#include <opm/polymer/polymermodel.hpp>
 
 #include <boost/filesystem/convenience.hpp>
 #include <boost/scoped_ptr.hpp>
@@ -67,7 +68,8 @@ public:
         : press_ (g->number_of_cells, 0.0),
           fpress_(g->number_of_faces, 0.0),
           flux_  (g->number_of_faces, 0.0),
-          sat_   (num_phases * g->number_of_cells, 0.0)
+          sat_   (num_phases * g->number_of_cells, 0.0),
+	  concentration_(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;
@@ -80,17 +82,20 @@ public:
     ::std::vector<double>& facepressure() { return fpress_; }
     ::std::vector<double>& faceflux    () { return flux_  ; }
     ::std::vector<double>& saturation  () { return sat_   ; }
+    ::std::vector<double>& concentration() { return concentration_; }
 
     const ::std::vector<double>& pressure    () const { return press_ ; }
     const ::std::vector<double>& facepressure() const { return fpress_; }
     const ::std::vector<double>& faceflux    () const { return flux_  ; }
     const ::std::vector<double>& saturation  () const { return sat_   ; }
+    const ::std::vector<double>& concentration() const { return concentration_; }
 
 private:
     ::std::vector<double> press_ ;
     ::std::vector<double> fpress_;
     ::std::vector<double> flux_  ;
     ::std::vector<double> sat_   ;
+    ::std::vector<double> concentration_   ;
 };
 
 
@@ -142,7 +147,9 @@ main(int argc, char** argv)
     bool use_deck = param.has("deck_filename");
     boost::scoped_ptr<Opm::Grid> grid;
     boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
+    PolymerData polydata;
     if (use_deck) {
+	THROW("We do not yet read polymer keywords from deck.");
 	std::string deck_filename = param.get<std::string>("deck_filename");
 	Opm::EclipseGridParser deck(deck_filename);
 	// Grid init
@@ -159,6 +166,14 @@ main(int argc, char** argv)
 	grid.reset(new Opm::Grid(nx, ny, nz));
 	// 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);
+	polydata.omega = param.getDefault("omega", 1.0);
+	polydata.c_vals.resize(2);
+	polydata.c_vals[0] = 0.0;
+	polydata.c_vals[0] = polydata.c_max_limit;
+	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);
     }
 
     // Extra rock init.
@@ -228,8 +243,10 @@ main(int argc, char** argv)
 			 stepsize,
 			 const_cast<UnstructuredGrid*>(grid->c_grid()),
 			 props.get(),
+			 &polydata,
 			 &state.faceflux()[0],
-			 &reorder_sat[0]);
+			 &reorder_sat[0],
+			 &state.concentration()[0]);
 	Opm::toBothSat(reorder_sat, state.saturation());
 
 	current_time += stepsize;

opm/polymer/polymermodel.cpp

@@ -1,30 +1,46 @@
 /* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
 /* Copyright 2012 (c) SINTEF */
 
-#include <stdlib.h>
-#include <stdio.h>
 
 #include <opm/polymer/polymermodel.hpp>
 #include <opm/core/grid.h>
 #include <opm/core/transport/reorder/nlsolvers.h>
 #include <opm/core/fluid/IncompPropertiesInterface.hpp>
 
+#include <cstdlib>
+#include <cstdio>
+#include <cmath>
+
 
 /* Parameters used in solution of single-cell boundary-value problem */
+
+
 struct Parameters
 {
+    double c;
     double s0;
-    double influx;  /* sum_j min(v_ij, 0)*f(s_j) */
+    // double c0;
-    double outflux; /* sum_j max(v_ij, 0)        */
     double dtpv;    /* dt/pv(i)                  */
+    double influx;  /* sum_j min(v_ij, 0)*f(s_j) */
+    // double influx_polymer;
+    double outflux; /* sum_j max(v_ij, 0)        */
+    // double dps;
+    // double rhor;
+    // double phi;
     int cell;
     const Opm::IncompPropertiesInterface* props;
+    const PolymerData* polydata;
 };
 
 
 static struct Parameters get_parameters(struct PolymerSolverData *d, int cell);
-static double residual(double s, void *data);
+static double residual_s(double s, void *data);
-static double fluxfun_props(double s, int cell, const Opm::IncompPropertiesInterface* props);
+static double residual_c(double c, void *data);
+static double fluxfun_props(double s,
+			    double c,
+			    int cell,
+			    const Opm::IncompPropertiesInterface* props, 
+			    const PolymerData* polydata);
 
 void
 destroy_solverdata(struct PolymerSolverData *d)
@@ -39,11 +55,13 @@ destroy_solverdata(struct PolymerSolverData *d)
 struct PolymerSolverData *
 init_solverdata(struct UnstructuredGrid *grid,
 		const Opm::IncompPropertiesInterface* props,
+		const PolymerData* polydata,
 		const double *darcyflux,
                 const double *porevolume,
 		const double *source,
                 const double dt,
-		double *saturation)
+		double *saturation,
+		double *concentration)
 {
     int i;
     struct PolymerSolverData *d = (struct PolymerSolverData*) malloc(sizeof *d);
@@ -52,12 +70,14 @@ init_solverdata(struct UnstructuredGrid *grid,
     {
         d->grid       = grid;
 	d->props      = props;
+	d->polydata   = polydata;
         d->darcyflux  = darcyflux;
         d->porevolume = porevolume;
         d->source     = source;
         d->dt         = dt;
 
         d->saturation     = saturation;
+	d->concentration  = concentration;
         d->fractionalflow = (double*) malloc(grid->number_of_cells *
 					     sizeof *d->fractionalflow);
         if (d->fractionalflow == NULL)
@@ -79,11 +99,11 @@ void polymer_solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell)
     struct PolymerSolverData   *d   = (struct PolymerSolverData*) data;
     struct Parameters   prm = get_parameters(d, cell);
 
-    d->saturation[cell] = find_zero(residual, &prm, ctrl);
+    d->saturation[cell] = find_zero(residual_s, &prm, ctrl);
     // double ff1 = fluxfun_props(d->saturation[cell], cell, d->props);
     // double ff2 = fluxfun(d->saturation[cell], -999);
     // printf("New = %f   old = %f\n", ff1, ff2);
-    d->fractionalflow[cell] = fluxfun_props(d->saturation[cell], cell, d->props);
+    d->fractionalflow[cell] = fluxfun_props(d->saturation[cell], d->concentration[cell], cell, d->props, d->polydata);
 }
 
 
@@ -96,10 +116,21 @@ void polymer_solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell)
  */
 /* influx is water influx, outflux is total outflux */
 static double
-residual(double s, void *data)
+residual_s(double s, void *data)
 {
     struct Parameters *p = (struct Parameters*) data;
-    return s - p->s0 +  p->dtpv*(p->outflux*fluxfun_props(s, p->cell, p->props) + p->influx);
+    double c = p->c;
+    return s - p->s0 +  p->dtpv*(p->outflux*fluxfun_props(s, c, p->cell, p->props, p->polydata) + p->influx);
+}
+
+static double
+residual_c(double c, void *data)
+{
+    (void) c;
+    (void) data;
+    // struct Parameters *p = (struct Parameters*) data;
+    //    return s - p->s0 +  p->dtpv*(p->outflux*fluxfun_props(s, p->cell, p->props) + p->influx);
+    return 0.0;
 }
 
 static struct Parameters
@@ -111,6 +142,7 @@ get_parameters(struct PolymerSolverData *d, int cell)
     double flux;
     int f, other;
 
+    p.c       = d->concentration[cell];
     p.s0      = d->saturation[cell];
     p.influx  = d->source[cell] >  0 ? -d->source[cell] : 0.0;
     p.outflux = d->source[cell] <= 0 ? -d->source[cell] : 0.0;
@@ -141,17 +173,31 @@ get_parameters(struct PolymerSolverData *d, int cell)
             }
         }
     }
+    p.polydata = d->polydata;
     return p;
 }
 
-static double fluxfun_props(double s, int cell, const Opm::IncompPropertiesInterface* props)
+static double fluxfun_props(double s, double c, int cell,
+			    const Opm::IncompPropertiesInterface* props,
+			    const PolymerData* pd)
 {
     const double* visc = props->viscosity();
+    double c_max_limit = pd->c_max_limit;
+    double cbar = c/c_max_limit;
+    double mu_w = visc[0];
+    double mu_m = pd->viscMult(c)*mu_w;
+    double mu_p = pd->viscMult(pd->c_max_limit)*mu_w;
+    double omega = pd->omega;
+    double mu_m_omega = std::pow(mu_m, omega);
+    double mu_w_e   = mu_m_omega*std::pow(mu_w, 1.0 - omega);
+    double mu_p_eff = mu_m_omega*std::pow(mu_p, 1.0 - omega);
+    double inv_mu_w_eff = (1.0 - cbar)/mu_w_e + cbar/mu_p_eff;
+    double inv_visc_eff[2] = { inv_mu_w_eff, 1.0/visc[1] };
     double sat[2] = { s, 1.0 - s };
     double mob[2];
     props->relperm(1, sat, &cell, mob, NULL);
-    mob[0] /= visc[0];
+    mob[0] *= inv_visc_eff[0];
-    mob[1] /= visc[1];
+    mob[1] *= inv_visc_eff[1];
     return mob[0]/(mob[0] + mob[1]);
 }
 

opm/polymer/polymermodel.hpp

@@ -4,7 +4,7 @@
 #ifndef POLYMER_HPP_INCLUDED
 #define POLYMER_HPP_INCLUDED
 
-
+#include <opm/core/utility/linearInterpolation.hpp>
 
 struct UnstructuredGrid;
 namespace Opm
@@ -12,14 +12,30 @@ namespace Opm
     class IncompPropertiesInterface;
 }
 
+struct PolymerData
+{
+    double c_max_limit;
+    double omega;
+    double viscMult(double c) const
+    {
+	return Opm::linearInterpolation(c_vals, visc_mult_vals, c);
+    }
+
+    std::vector<double> c_vals;
+    std::vector<double> visc_mult_vals;
+};
+
+
 struct PolymerSolverData  {
     struct UnstructuredGrid *grid;
     const Opm::IncompPropertiesInterface* props;
+    const PolymerData* polydata;
     const double            *darcyflux;   /* one flux per face  in cdata::grid*/
     const double            *porevolume;  /* one volume per cell */
     const double            *source;      /* one source per cell */
     double                   dt;
     double                  *saturation;      /* one per cell */
+    double                  *concentration;   /* one per cell */
     double                  *fractionalflow;  /* one per cell */
 };
 
@@ -35,11 +51,13 @@ destroy_solverdata(struct PolymerSolverData *d);
 struct PolymerSolverData *
 init_solverdata(struct UnstructuredGrid *grid,
 		const Opm::IncompPropertiesInterface* props,
+		const PolymerData* polydata,
 		const double *darcyflux,
                 const double *porevolume,
 		const double *source,
                 const double dt,
-		double *saturation);
+		double *saturation,
+		double *concentration);
 
 
 #endif /* POLYMER_H_INCLUDED */

opm/polymer/polymertransport.cpp

@@ -20,8 +20,10 @@ void polymertransport(
     double                   dt,
     struct UnstructuredGrid *grid,
     const Opm::IncompPropertiesInterface* props,
+    const PolymerData* polydata,
     const double            *darcyflux,
-    double                  *saturation)
+    double                  *saturation,
+    double                  *concentration)
 {
     int    i;
 
@@ -29,8 +31,9 @@ void polymertransport(
     int    *sequence;
     int    *components;
 
-    struct PolymerSolverData *data = init_solverdata(grid, props, darcyflux,
+    PolymerSolverData *data = init_solverdata(grid, props, polydata, darcyflux,
-                                              porevolume, source, dt, saturation);
+					      porevolume, source, dt, saturation,
+					      concentration);
 
 
     struct NonlinearSolverCtrl ctrl;

opm/polymer/polymertransport.hpp

@@ -10,6 +10,7 @@ namespace Opm
 }
 
 struct UnstructuredGrid;
+struct PolymerData;
 
 void polymertransport(
     const double *porevolume,
@@ -17,8 +18,10 @@ void polymertransport(
     double dt,
     struct UnstructuredGrid *grid,
     const Opm::IncompPropertiesInterface* props,
+    const PolymerData* polydata,
     const double *darcyflux,
-    double *saturation);
+    double *saturation,
+    double *concentration);
 
 
 #endif /* POLYMERTRANSPORT_HPP_INCLUDED */