diff --git a/Makefile.am b/Makefile.am index 31bb00fd7..4942ceb8b 100644 --- a/Makefile.am +++ b/Makefile.am @@ -2,5 +2,14 @@ ACLOCAL_AMFLAGS = -I m4 SUBDIRS = . examples + +# Declare libraries +lib_LTLIBRARIES = libopmpolymer.la + +libopmpolymer_la_SOURCES = \ +opm/polymer/polymermodel.cpp \ +opm/polymer/polymertransport.cpp + nobase_include_HEADERS = \ -opm/polymer/HelloPolymer.hpp +opm/polymer/polymermodel.hpp \ +opm/polymer/polymertransport.hpp diff --git a/examples/Makefile.am b/examples/Makefile.am index f17a3ce17..d8324aa03 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,6 +1,8 @@ AM_CPPFLAGS = \ -I$(top_srcdir) +LDADD = $(top_builddir)/libopmpolymer.la + noinst_PROGRAMS = \ hello \ polymer_reorder diff --git a/examples/polymer_reorder.cpp b/examples/polymer_reorder.cpp index 685de3a96..17cd84a62 100644 --- a/examples/polymer_reorder.cpp +++ b/examples/polymer_reorder.cpp @@ -37,7 +37,7 @@ #include -#include +#include #include #include @@ -223,13 +223,13 @@ main(int argc, char** argv) // Also, for anything but noflow boundaries, // boundary flows must be accumulated into // source term following the same convention. - twophasetransport(&porevol[0], - &reorder_src[0], - stepsize, - const_cast(grid->c_grid()), - props.get(), - &state.faceflux()[0], - &reorder_sat[0]); + polymertransport(&porevol[0], + &reorder_src[0], + stepsize, + const_cast(grid->c_grid()), + props.get(), + &state.faceflux()[0], + &reorder_sat[0]); Opm::toBothSat(reorder_sat, state.saturation()); current_time += stepsize; diff --git a/opm/polymer/polymermodel.cpp b/opm/polymer/polymermodel.cpp new file mode 100644 index 000000000..8827150e1 --- /dev/null +++ b/opm/polymer/polymermodel.cpp @@ -0,0 +1,161 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ +/* Copyright 2012 (c) SINTEF */ + +#include +#include + +#include +#include +#include +#include + + +/* Parameters used in solution of single-cell boundary-value problem */ +struct Parameters +{ + double s0; + double influx; /* sum_j min(v_ij, 0)*f(s_j) */ + double outflux; /* sum_j max(v_ij, 0) */ + double dtpv; /* dt/pv(i) */ + int cell; + const Opm::IncompPropertiesInterface* props; +}; + + +static struct Parameters get_parameters(struct PolymerSolverData *d, int cell); +static double residual(double s, void *data); +static double fluxfun_props(double s, int cell, const Opm::IncompPropertiesInterface* props); + +void +destroy_solverdata(struct PolymerSolverData *d) +{ + if (d!=NULL) + { + free(d->fractionalflow); + } + free(d); +} + +struct PolymerSolverData * +init_solverdata(struct UnstructuredGrid *grid, + const Opm::IncompPropertiesInterface* props, + const double *darcyflux, + const double *porevolume, + const double *source, + const double dt, + double *saturation) +{ + int i; + struct PolymerSolverData *d = (struct PolymerSolverData*) malloc(sizeof *d); + + if(d!=NULL) + { + d->grid = grid; + d->props = props; + d->darcyflux = darcyflux; + d->porevolume = porevolume; + d->source = source; + d->dt = dt; + + d->saturation = saturation; + d->fractionalflow = (double*) malloc(grid->number_of_cells * + sizeof *d->fractionalflow); + if (d->fractionalflow == NULL) + { + destroy_solverdata(d); + d = NULL; + } + for(i=0; inumber_of_cells; ++i) + { + d->fractionalflow[i] = 0.0; + } + } + return d; +} + +/* Solver for single-cell bvp calls root-finder in nlsolvers.c */ +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); + // 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); +} + + +/* ====================== Internals =================================*/ + + +/* Residual function r(s) for a single-cell bvp */ +/* + * r(s) = s - s0 + dt/pv*(influx - outflux*f(s) ) + */ +/* influx is water influx, outflux is total outflux */ +static double +residual(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); +} + +static struct Parameters +get_parameters(struct PolymerSolverData *d, int cell) +{ + int i; + struct UnstructuredGrid *g = d->grid; + struct Parameters p; + double flux; + int f, other; + + 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; + p.dtpv = d->dt/d->porevolume[cell]; + p.cell = cell; + p.props = d->props; + + d->saturation[cell] = 0; + for (i=g->cell_facepos[cell]; icell_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]; + } + else { + p.outflux += flux; + } + } + } + return p; +} + +static double fluxfun_props(double s, int cell, const Opm::IncompPropertiesInterface* props) +{ + const double* visc = props->viscosity(); + double sat[2] = { s, 1.0 - s }; + double mob[2]; + props->relperm(1, sat, &cell, mob, NULL); + mob[0] /= visc[0]; + mob[1] /= visc[1]; + return mob[0]/(mob[0] + mob[1]); +} + + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/polymer/polymermodel.hpp b/opm/polymer/polymermodel.hpp new file mode 100644 index 000000000..6b1e1b425 --- /dev/null +++ b/opm/polymer/polymermodel.hpp @@ -0,0 +1,49 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ +/* Copyright 2012 (c) SINTEF */ + +#ifndef POLYMER_HPP_INCLUDED +#define POLYMER_HPP_INCLUDED + + + +struct UnstructuredGrid; +namespace Opm +{ + class IncompPropertiesInterface; +} + +struct PolymerSolverData { + struct UnstructuredGrid *grid; + const Opm::IncompPropertiesInterface* props; + 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 *fractionalflow; /* one per cell */ +}; + +struct NonlinearSolverCtrl; + + +void +polymer_solvecell(void *data, struct NonlinearSolverCtrl *ctrl, int cell); + +void +destroy_solverdata(struct PolymerSolverData *d); + +struct PolymerSolverData * +init_solverdata(struct UnstructuredGrid *grid, + const Opm::IncompPropertiesInterface* props, + const double *darcyflux, + const double *porevolume, + const double *source, + const double dt, + double *saturation); + + +#endif /* POLYMER_H_INCLUDED */ + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/polymer/polymertransport.cpp b/opm/polymer/polymertransport.cpp new file mode 100644 index 000000000..dd9d0cb4b --- /dev/null +++ b/opm/polymer/polymertransport.cpp @@ -0,0 +1,77 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ +/* Copyright 2012 (c) SINTEF */ + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + + +void polymertransport( + const double *porevolume, + const double *source, + double dt, + struct UnstructuredGrid *grid, + const Opm::IncompPropertiesInterface* props, + const double *darcyflux, + double *saturation) +{ + int i; + + int ncomponents; + int *sequence; + int *components; + + struct PolymerSolverData *data = init_solverdata(grid, props, darcyflux, + porevolume, source, dt, saturation); + + + 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); + + 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) + { +#ifdef MATLAB_MEX_FILE + if (interrupt_signal) + { + mexPrintf("Reorder loop interrupted by user: %d of %d " + "cells finished.\n", i, grid->number_of_cells); + break; + } +#endif + polymer_solvecell(data, &ctrl, sequence[i]); + /*print("iterations:%d residual:%20.16f\n", + * ctrl.iterations, ctrl.residual);*/ + } + + destroy_solverdata(data); + + free(sequence); + free(components); +} + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/polymer/polymertransport.hpp b/opm/polymer/polymertransport.hpp new file mode 100644 index 000000000..c642abd1a --- /dev/null +++ b/opm/polymer/polymertransport.hpp @@ -0,0 +1,24 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ +/* Copyright 2012 (c) SINTEF */ + +#ifndef POLYMERTRANSPORT_HPP_INCLUDED +#define POLYMERTRANSPORT_HPP_INCLUDED + +namespace Opm +{ + class IncompPropertiesInterface; +} + +struct UnstructuredGrid; + +void polymertransport( + const double *porevolume, + const double *source, + double dt, + struct UnstructuredGrid *grid, + const Opm::IncompPropertiesInterface* props, + const double *darcyflux, + double *saturation); + + +#endif /* POLYMERTRANSPORT_HPP_INCLUDED */