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

First compiling version of full polymer residual. Not yet working.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-02-06 14:01:56 +01:00
parent 0ad4fab860
commit d17cfd1372
5 changed files with 160 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
examples/polymer_reorder.cpp
View File

@ -63,14 +63,16 @@
class ReservoirState {
class ReservoirState
{
public:
ReservoirState(const UnstructuredGrid* g, const int num_phases = 2)
: 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),
concentration_(g->number_of_cells, 0.0)
concentration_(g->number_of_cells, 0.0),
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;
@ -79,24 +81,27 @@ public:
int numPhases() const { return sat_.size()/press_.size(); }
::std::vector<double>& pressure () { return press_ ; }
::std::vector<double>& facepressure() { return fpress_; }
::std::vector<double>& faceflux () { return flux_ ; }
::std::vector<double>& saturation () { return sat_ ; }
::std::vector<double>& concentration() { return concentration_; }
std::vector<double>& pressure () { return press_ ; }
std::vector<double>& facepressure() { return fpress_; }
std::vector<double>& faceflux () { return flux_ ; }
std::vector<double>& saturation () { return sat_ ; }
std::vector<double>& concentration() { return concentration_; }
std::vector<double>& cmax() { return cmax_; }
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_; }
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_; }
const std::vector<double>& cmax() const { return cmax_; }
private:
::std::vector<double> press_ ;
::std::vector<double> fpress_;
::std::vector<double> flux_ ;
::std::vector<double> sat_ ;
::std::vector<double> concentration_ ;
std::vector<double> press_ ;
std::vector<double> fpress_;
std::vector<double> flux_ ;
std::vector<double> sat_ ;
std::vector<double> concentration_;
std::vector<double> cmax_;
};
@ -169,12 +174,18 @@ main(int argc, char** argv)
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.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.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", 30.0);
}
// Extra rock init.
@ -252,6 +263,7 @@ main(int argc, char** argv)
// source term following the same convention.
transport_timer.start();
polymertransport(&porevol[0],
props->porosity(),
&reorder_src[0],
stepsize,
const_cast<UnstructuredGrid*>(grid->c_grid()),
@ -259,7 +271,8 @@ main(int argc, char** argv)
&polydata,
&state.faceflux()[0],
&reorder_sat[0],
&state.concentration()[0]);
&state.concentration()[0],
&state.cmax()[0]);
Opm::toBothSat(reorder_sat, state.saturation());
transport_timer.stop();
double tt = transport_timer.secsSinceStart();

121
opm/polymer/polymermodel.cpp
View File

@ -19,14 +19,9 @@ struct ParametersSRes
{
double c;
double s0;
// double c0;
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;
@ -34,10 +29,19 @@ struct ParametersSRes
struct ParametersCRes
{
double s0;
double c0;
double cmax0;
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 porosity;
PolymerSolverData* psdata;
int cell;
NonlinearSolverCtrl* ctrl;
double s;
const PolymerData* polydata;
};
@ -51,6 +55,9 @@ static double fluxfun_props(double s,
int cell,
const Opm::IncompPropertiesInterface* props,
const PolymerData* polydata);
static double compute_mc(double c,
const Opm::IncompPropertiesInterface* props,
const PolymerData* polydata);
void
destroy_solverdata(struct PolymerSolverData *d)
@ -58,6 +65,7 @@ destroy_solverdata(struct PolymerSolverData *d)
if (d!=NULL)
{
free(d->fractionalflow);
free(d->mc);
}
free(d);
}
@ -68,10 +76,12 @@ init_solverdata(struct UnstructuredGrid *grid,
const PolymerData* polydata,
const double *darcyflux,
const double *porevolume,
const double *porosity,
const double *source,
const double dt,
double *saturation,
double *concentration)
double *concentration,
double *cmax)
{
int i;
struct PolymerSolverData *d = (struct PolymerSolverData*) malloc(sizeof *d);
@ -83,14 +93,18 @@ init_solverdata(struct UnstructuredGrid *grid,
d->polydata = polydata;
d->darcyflux = darcyflux;
d->porevolume = porevolume;
d->porosity = porosity;
d->source = source;
d->dt = dt;
d->saturation = saturation;
d->concentration = concentration;
d->cmax = cmax;
d->fractionalflow = (double*) malloc(grid->number_of_cells *
sizeof *d->fractionalflow);
if (d->fractionalflow == NULL)
d->mc = (double*) malloc(grid->number_of_cells *
sizeof *d->mc);
if (d->fractionalflow == NULL || d->mc == NULL)
{
destroy_solverdata(d);
d = NULL;
@ -110,8 +124,10 @@ void polymer_solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell)
struct ParametersCRes prm = get_parameters_c(d, cell, ctrl);
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);
d->mc[cell] = compute_mc(d->concentration[cell], d->props, d->polydata);
}
@ -134,13 +150,25 @@ residual_s(double s, void *data)
static double
residual_c(double c, void *data)
{
struct ParametersCRes *prm_c = (struct ParametersCRes*) data;
int cell = prm_c->cell;
struct ParametersSRes prm_s = get_parameters_s(prm_c->psdata, cell);
struct ParametersCRes *p = (struct ParametersCRes*) 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, prm_c->ctrl);
prm_c->s = s;
return c - 0.0;
double s = find_zero(residual_s, &prm_s, p->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);
double dps = p->polydata->dps;
double s0 = p->s0;
double c0 = p->c0;
double rhor = p->polydata->rhor;
double porosity = p->porosity;
double cmax0 = p->cmax0;
double ads0 = p->polydata->adsorbtion(std::max(c0, cmax0));
double ads = p->polydata->adsorbtion(std::max(c, cmax0));
return (s - dps)*c - (s0 - dps)*c0
+ rhor*((1.0 - porosity)/porosity)*(ads - ads0)
+ p->dtpv*(p->outflux*ff*mc + p->influx_polymer);
}
static struct ParametersSRes
@ -190,12 +218,49 @@ get_parameters_s(struct PolymerSolverData *d, int cell)
static struct ParametersCRes
get_parameters_c(struct PolymerSolverData *d, int cell, NonlinearSolverCtrl* ctrl)
{
ParametersCRes prm;
prm.psdata = d;
prm.cell = cell;
prm.ctrl = ctrl;
prm.s = -1e100;
return prm;
int i;
struct UnstructuredGrid *g = d->grid;
double flux;
int f, other;
ParametersCRes p;
p.c0 = d->concentration[cell];
p.cmax0 = d->cmax[cell];
p.s0 = d->saturation[cell];
p.dtpv = d->dt/d->porevolume[cell];
p.influx = d->source[cell] > 0 ? -d->source[cell] : 0.0;
p.influx_polymer = d->source[cell] > 0 ? -d->source[cell] : 0.0; // TODO. Wrong if nonzero, mult by mc.
for (i=g->cell_facepos[cell]; i<g->cell_facepos[cell+1]; ++i) {
f = g->cell_faces[i];
/* Compute cell flux*/
if (cell == g->face_cells[2*f]) {
flux = d->darcyflux[f];
other = g->face_cells[2*f+1];
}
else {
flux =-d->darcyflux[f];
other = g->face_cells[2*f];
}
if (other != -1) {
if (flux < 0.0) {
p.influx += flux*d->fractionalflow[other];
p.influx_polymer += flux*d->fractionalflow[other]*d->mc[other];
}
else {
p.outflux += flux;
}
}
}
p.outflux = d->source[cell] <= 0 ? -d->source[cell] : 0.0;
p.porosity = d->porosity[cell];
p.psdata = d;
p.cell = cell;
p.ctrl = ctrl;
p.s = -1e100;
p.polydata = d->polydata;
return p;
}
@ -223,6 +288,24 @@ static double fluxfun_props(double s, double c, int cell,
return mob[0]/(mob[0] + mob[1]);
}
static double compute_mc(double c,
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;
return c/(inv_mu_w_eff*mu_p_eff);
}
/* Local Variables: */
/* c-basic-offset:4 */

19
opm/polymer/polymermodel.hpp
View File

@ -18,11 +18,19 @@ struct PolymerData
double omega;
double viscMult(double c) const
{
return Opm::linearInterpolation(c_vals, visc_mult_vals, c);
return Opm::linearInterpolation(c_vals_visc, visc_mult_vals, c);
}
double rhor;
double dps;
double adsorbtion(double c) const
{
return Opm::linearInterpolation(c_vals_ads, ads_vals, c);
}
std::vector<double> c_vals;
std::vector<double> c_vals_visc;
std::vector<double> visc_mult_vals;
std::vector<double> c_vals_ads;
std::vector<double> ads_vals;
};
@ -32,11 +40,14 @@ struct PolymerSolverData {
const PolymerData* polydata;
const double *darcyflux; /* one flux per face in cdata::grid*/
const double *porevolume; /* one volume per cell */
const double *porosity;
const double *source; /* one source per cell */
double dt;
double *saturation; /* one per cell */
double *concentration; /* one per cell */
double *cmax;
double *fractionalflow; /* one per cell */
double *mc;
};
struct NonlinearSolverCtrl;
@ -54,10 +65,12 @@ init_solverdata(struct UnstructuredGrid *grid,
const PolymerData* polydata,
const double *darcyflux,
const double *porevolume,
const double *porosity,
const double *source,
const double dt,
double *saturation,
double *concentration);
double *concentration,
double *cmax);
#endif /* POLYMER_H_INCLUDED */

8
opm/polymer/polymertransport.cpp
View File

@ -16,6 +16,7 @@
void polymertransport(
const double *porevolume,
const double *porosity,
const double *source,
double dt,
struct UnstructuredGrid *grid,
@ -23,7 +24,8 @@ void polymertransport(
const PolymerData* polydata,
const double *darcyflux,
double *saturation,
double *concentration)
double *concentration,
double *cmax)
{
int i;
@ -32,8 +34,8 @@ void polymertransport(
int *components;
PolymerSolverData *data = init_solverdata(grid, props, polydata, darcyflux,
porevolume, source, dt, saturation,
concentration);
porevolume, porosity, source, dt, saturation,
concentration, cmax);
struct NonlinearSolverCtrl ctrl;

4
opm/polymer/polymertransport.hpp
View File

@ -14,6 +14,7 @@ struct PolymerData;
void polymertransport(
const double *porevolume,
const double *porosity,
const double *source,
double dt,
struct UnstructuredGrid *grid,
@ -21,7 +22,8 @@ void polymertransport(
const PolymerData* polydata,
const double *darcyflux,
double *saturation,
double *concentration);
double *concentration,
double *cmax);
#endif /* POLYMERTRANSPORT_HPP_INCLUDED */