Copied and renamed files and functions from opm-core as base for polymer reorder solver.

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

examples/Makefile.am

@@ -1,6 +1,8 @@
 AM_CPPFLAGS =					\
 -I$(top_srcdir)
 
+LDADD = $(top_builddir)/libopmpolymer.la
+
 noinst_PROGRAMS =				\
 hello						\
 polymer_reorder

examples/polymer_reorder.cpp

@@ -37,7 +37,7 @@
 
 #include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
 
-#include <opm/core/transport/reorder/twophasetransport.hpp>
+#include <opm/polymer/polymertransport.hpp>
 
 #include <boost/filesystem/convenience.hpp>
 #include <boost/scoped_ptr.hpp>
@@ -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<UnstructuredGrid*>(grid->c_grid()),
-			  props.get(),
-			  &state.faceflux()[0],
-			  &reorder_sat[0]);
+	polymertransport(&porevol[0],
+			 &reorder_src[0],
+			 stepsize,
+			 const_cast<UnstructuredGrid*>(grid->c_grid()),
+			 props.get(),
+			 &state.faceflux()[0],
+			 &reorder_sat[0]);
 	Opm::toBothSat(reorder_sat, state.saturation());
 
 	current_time += stepsize;

opm/polymer/polymermodel.cpp

@@ -0,0 +1,161 @@
+/* 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>
+
+
+/* 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; i<grid->number_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]; 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];
+            }
+            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:                */

opm/polymer/polymermodel.hpp

@@ -0,0 +1,49 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+/* 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:                */

opm/polymer/polymertransport.cpp

@@ -0,0 +1,77 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+/* Copyright 2012 (c) SINTEF */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <assert.h>
+
+#include <opm/polymer/polymermodel.hpp>
+#include <opm/polymer/polymertransport.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/transport/reorder/reordersequence.h>
+#include <opm/core/transport/reorder/tarjan.h>
+#include <opm/core/transport/reorder/nlsolvers.h>
+
+
+
+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; i<grid->number_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:                */

opm/polymer/polymertransport.hpp

@@ -0,0 +1,24 @@
+/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
+/* 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 */