From 71f6bac21e55cd2b0f121e314c5843a48fe779b1 Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Tue, 17 Jan 2012 14:39:09 +0100 Subject: [PATCH] Add bare-bones implementation of reordering algorithm. Lightly tested. --- Makefile.am | 17 +- opm/core/transport/reorder/nlsolvers.c | 263 ++++++++++++++++++ opm/core/transport/reorder/nlsolvers.h | 14 + opm/core/transport/reorder/reordersequence.c | 103 +++++++ opm/core/transport/reorder/reordersequence.h | 13 + opm/core/transport/reorder/tarjan.c | 180 ++++++++++++ opm/core/transport/reorder/tarjan.h | 19 ++ opm/core/transport/reorder/twophase.c | 123 ++++++++ opm/core/transport/reorder/twophase.h | 26 ++ .../transport/reorder/twophasetransport.c | 54 ++++ .../transport/reorder/twophasetransport.h | 16 ++ 11 files changed, 825 insertions(+), 3 deletions(-) create mode 100644 opm/core/transport/reorder/nlsolvers.c create mode 100644 opm/core/transport/reorder/nlsolvers.h create mode 100644 opm/core/transport/reorder/reordersequence.c create mode 100644 opm/core/transport/reorder/reordersequence.h create mode 100644 opm/core/transport/reorder/tarjan.c create mode 100644 opm/core/transport/reorder/tarjan.h create mode 100644 opm/core/transport/reorder/twophase.c create mode 100644 opm/core/transport/reorder/twophase.h create mode 100644 opm/core/transport/reorder/twophasetransport.c create mode 100644 opm/core/transport/reorder/twophasetransport.h diff --git a/Makefile.am b/Makefile.am index 4991ade9..3d1e5b2b 100644 --- a/Makefile.am +++ b/Makefile.am @@ -73,8 +73,14 @@ opm/core/pressure/mimetic/mimetic.c \ opm/core/pressure/mimetic/hybsys_global.c \ opm/core/pressure/mimetic/hybsys.c \ opm/core/transport/spu_explicit.c \ -opm/core/transport/spu_implicit.c \ -opm/core/transport/transport_source.c +opm/core/transport/spu_implicit.c \ +opm/core/transport/transport_source.c \ +opm/core/transport/reorder/reordersequence.c \ +opm/core/transport/reorder/twophase.c \ +opm/core/transport/reorder/nlsolvers.c \ +opm/core/transport/reorder/tarjan.c \ +opm/core/transport/reorder/twophasetransport.c + nobase_include_HEADERS = \ opm/core/eclipse/CornerpointChopper.hpp \ @@ -177,7 +183,12 @@ opm/core/transport/spu_explicit.h \ opm/core/transport/JacobianSystem.hpp \ opm/core/transport/spu_implicit.h \ opm/core/transport/ImplicitAssembly.hpp \ -opm/core/transport/transport_source.h +opm/core/transport/transport_source.h \ +opm/core/transport/reorder/nlsolvers.h \ +opm/core/transport/reorder/reordersequence.h \ +opm/core/transport/reorder/twophase.h \ +opm/core/transport/reorder/tarjan.h \ +opm/core/transport/reorder/twophasetransport.h if UMFPACK libopmcore_la_SOURCES += opm/core/linalg/call_umfpack.c diff --git a/opm/core/transport/reorder/nlsolvers.c b/opm/core/transport/reorder/nlsolvers.c new file mode 100644 index 00000000..5da76dfc --- /dev/null +++ b/opm/core/transport/reorder/nlsolvers.c @@ -0,0 +1,263 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#include +#include +#include +#include +#include + +#include "nlsolvers.h" + + +#ifdef MEXFILE +#include +extern int interrupted; +#define print mexPrintf +#define malloc mxMalloc +#define calloc mxCalloc +#define realloc mxRealloc +#define free mxFree +#else +#define print printf +#endif + + +static const char no_root_str[]= + " In %s:\n" + " With G(0) =% 5f, G(1) =% 5f, G(s) cannot have a zero in [0,1]!\n"; + +static const double EPS = 1e-14; + + +/* Start with bracket [0,1] with G(0)*G(1)<0. */ +/*---------------------------------------------------------------------------*/ +double +ridder (double (*G)(double, void*), void *data, double tol, int maxit, int *iterations) +/*---------------------------------------------------------------------------*/ +{ + *iterations = 0; + + double s2 = ((double *) data)[2]; + double G2 = G(s2, data); + if (fabs(G2) 0) + { + print(no_root_str, "ridder", G0, G1); + return -1.0; + } + + if (G0>G1) + { + double swap; + swap = G0; + G0 = G1; + G1 = swap; + } + + double s3=0; + double G3=10; + + int it = 0; + while ( (fabs(G3) > tol) && (it++ < maxit)) + { + /* find zero crossing of line segment [(s0,G0), (s1,G1)] */ + double root = sqrt(G2*G2-G0*G1); + if (fabs(root)G1 ? 1.0 : -1.0; + s3 = s2 + ( s2-s0 )*sgn*G2/root; + G3 = G(s3, data); + + + /* if G2*G3<0 */ + if (G2*G3 <= 0.0) + { + if (G2 > G3) + { + s0 = s3; + G0 = G3; + s1 = s2; + G1 = G2; + } + else + { + s0 = s2; + G0 = G2; + s1 = s3; + G1 = G3; + + } + + } + else if (G0*G3 <= 0.0) + { + s1 = s3; + G1 = G3; + } + else if (G1*G3 <= 0.0) + { + s0 = s3; + G0 = G3; + } + else + { + print( + "In ridder:\n" + "G0=%10.10f, G1=%10.10f, G3=%10.10f\n", + G0, G1, G3 + ); + getchar(); + } + s2 = 0.5*(s0+s1); + G2 = G(s2, data); + } + + *iterations = it; + return s3; +} + + + +/* Start with bracket [0,1] with G(0)*G(1)<0. Search by finding zero crossing + sN of line segment [(sL,GL), (sR,GR)]. Set SL=sN if G(sN<0), sR=sN + otherwise.*/ +/*---------------------------------------------------------------------------*/ +double +regulafalsi (double (*G)(double, void*), void *data, double tol, int maxit, + int *iterations) +/*---------------------------------------------------------------------------*/ +{ + *iterations = 0; + + double sn = 0.5;//((double *) data)[2]; /* "Undefined" value */ + double Gn = G(sn, data); + if (fabs(Gn) < tol) return sn; + + /* Initial guess is interval [0,1] of course */ + double s0=0.0; + double G0=G(s0, data); + if (fabs(G0) 0) { + print(no_root_str, "regulafalsi", G0, G1); + return -1.0; + } + + if (G0>G1) { + double swap; + swap = G0; + G0 = G1; + G1 = swap; + } + + int it = 0; + while ( (fabs(Gn) > tol) && (it++ < maxit)) { +#if 0 + /* Unmodified Regula-Falsi */ + /* maintain bracket with G1>G0 */ + if ( Gn>0 ) { + G1 = Gn; + s1 = sn; + } + else { + G0 = Gn; + s0 = sn; + } +#else + /* Modified Regula-Falsi*/ + if ((Gn>0.0)==(G0>0.0)) { + s0 = s1; + G0 = G1; + } + else { + // const double gamma_illinois = 0.5; + const double gamma_pegasus = G1/(G1+Gn); + G0 *= gamma_pegasus; + } + s1 = sn; + G1 = Gn; +#endif + + sn = s0 - (s1-s0)*G0/(G1-G0); + Gn = G(sn, data); + } + + *iterations = it; + return sn; +} + + + +/* Start with bracket [0,1] with G(0)*G(1)<0. Search by finding sN=0.5(sL+sR). + Set SL=sN if G(sN<0), sR=sN otherwise.*/ +/*---------------------------------------------------------------------------*/ +double +bisection (double (*G)(double, void*), void *data, double tol, int maxit, + int *iterations) +/*---------------------------------------------------------------------------*/ +{ + *iterations = 0; + + double sn = ((double *) data)[2]; + double Gn = G(sn, data); + if (fabs(Gn) 0.0) { + print(no_root_str, "bisection", G0, G1); + return -1.0; + } + + if (G0>G1) { + double swap; + swap = G0; + G0 = G1; + G1 = swap; + } + + int it=0; + while ( (fabs(Gn)>tol) && (it++ < maxit) ) { + if ( Gn>0 ) { + G1 = Gn; + s1 = sn; + } + else { + G0 = Gn; + s0 = sn; + } + + sn = 0.5*(s0+s1); + Gn = G(sn, data); + } + *iterations = it; + if (it >= maxit) print("Warning: convergence criterion not met\n"); + return sn; +} + + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/core/transport/reorder/nlsolvers.h b/opm/core/transport/reorder/nlsolvers.h new file mode 100644 index 00000000..1774218d --- /dev/null +++ b/opm/core/transport/reorder/nlsolvers.h @@ -0,0 +1,14 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#ifndef NLSOLVERS_H +#define NLSOLVERS_H + +double bisection (double (*)(double, void*), void*, double, int, int*); +double ridder (double (*)(double, void*), void*, double, int, int*); +double regulafalsi (double (*)(double, void*), void*, double, int, int*); + +#endif /* NLSOLVERS_H */ + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/core/transport/reorder/reordersequence.c b/opm/core/transport/reorder/reordersequence.c new file mode 100644 index 00000000..57029f42 --- /dev/null +++ b/opm/core/transport/reorder/reordersequence.c @@ -0,0 +1,103 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#include +#include +#include +#include +#include + + + +/* Construct adjacency matrix of upwind graph wrt flux. Column + indices are not sorted. */ +static void +make_upwind_graph(int nc, int *cellfaces, int *faceptr, int *face2cell, + const double *flux, int *ia, int *ja, int *work) +{ + /* Using topology (conn, cptr), and direction, construct adjacency + matrix of graph. */ + + int i, j, p, f, positive_sign; + double theflux; + + /* For each face, store upwind cell in work array */ + for (i=0; i 0 ) { + /* flux from cell i over face f */ + work[f] = i; + } + } + } + + /* Fill ia and ja */ + p = 0; + ia[0] = p; + for (i=0; inumber_of_cells, + grid->number_of_faces, + grid->cell_faces, + grid->cell_facepos, + grid->face_cells, + flux, + sequence, + components, + ncomponents); +} + + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ + + diff --git a/opm/core/transport/reorder/reordersequence.h b/opm/core/transport/reorder/reordersequence.h new file mode 100644 index 00000000..d4344c2a --- /dev/null +++ b/opm/core/transport/reorder/reordersequence.h @@ -0,0 +1,13 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ +#ifndef REORDERSEQUENCE_H_INCLUDED +#define REORDERSEQUENCE_H_INCLUDED + +struct UnstructuredGrid; +void compute_sequence(struct UnstructuredGrid *grid, const double *flux, + int *sequence, int *components, int *ncomponents); + +#endif + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/core/transport/reorder/tarjan.c b/opm/core/transport/reorder/tarjan.c new file mode 100644 index 00000000..f4f1a5a2 --- /dev/null +++ b/opm/core/transport/reorder/tarjan.c @@ -0,0 +1,180 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#include +#include +#include + + + + + +static int min(int a, int b){ return a= -2); + + if (status[c] == REMAINING) + { + status[c] = ia[c+1]-ia[c]; /* number of descendants of c */ + time[c] = link[c] = t++; + *cstack-- = c; /* push c on strongcomp stack */ + } + + + + /* if all descendants are processed */ + if (status[c] == 0) + { + + /* if c is root of strong component */ + if (link[c] == time[c]) + { + do + { + assert (cstack != cbottom); + + v = *++cstack; /* pop strong component stack */ + status[v] = DONE; + P[pos++] = v; /* store vertex in P */ + } + while ( v != c ); + + *Q++ = pos; /* store end point of component */ + ++*ncomp; + } + + ++stack; /* pop c */ + + if (stack != bottom) + { + link[*(stack+1)] = min(link[*(stack+1)], link[c]); + } + } + + + + /* if there are more descendants to consider */ + else + { + assert(status[c] > 0); + + child = ja[ia[c] + status[c]-1]; + --status[c]; /* Decrement descendant count of c*/ + + if (status[child] == REMAINING) + { + *stack-- = child; /* push child */ + + } + else if (status[child] >= 0) + { + link[c] = min(link[c], time[child]); + + } + else + { + assert(status[child] = DONE); + } + } + } + assert (cstack == cbottom); + } +} + + +#ifdef TEST_TARJAN +#include +#include +/* gcc -DTEST_TARJAN tarjan.c */ +int main() +{ + int size = 5; + int ia[] = {0, 1, 2, 4, 5, 5}; + int ja[] = {1, 2, 0, 3, 4}; + + int *work = malloc(3*size*sizeof *work); + + int number_of_components; + int *permutation = malloc(size * sizeof *permutation); + int *ptr = malloc(size * sizeof *ptr); + int i,j; + tarjan (size, ia, ja, permutation, ptr, &number_of_components, work); + + fprintf(stderr, "Number of components: %d\n", number_of_components); + for (i=0; i */ + +#ifndef TARJAN_H_INCLUDED +#define TARJAN_H_INCLUDED +#ifdef __cplusplus +extern "C" { +#endif + + void tarjan (int size, int *ia, int *ja, int *rowP, int *P, + int *ncomp, int *work); + +#ifdef __cplusplus +} +#endif +#endif /* TARJAN_H_INCLUDED */ + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/core/transport/reorder/twophase.c b/opm/core/transport/reorder/twophase.c new file mode 100644 index 00000000..3b172ae2 --- /dev/null +++ b/opm/core/transport/reorder/twophase.c @@ -0,0 +1,123 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#include +#include +#include + +static struct Parameters get_parameters(struct vdata *vd, + const struct cdata *cd, + int cell); +static double fluxfun(double s); + +static double +fluxfun(double s) +{ + return s*s/(s*s + (1-s)*(1-s)); +} + +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) */ +}; +#include +static double +G(double s, void *data) +{ + struct Parameters *p = data; + fprintf(stderr, "s:%f dtpv:%f s0:%f, influx:%f outflux:%f G:%f\n", + s, + p->dtpv, + p->s0, + p->influx, + p->outflux, + s - p->s0 + p->dtpv*(p->outflux*fluxfun(s) + p->influx)); + /* G(s) = s - s0 + dt/pv*(influx - outflux*f(s) ) */ + return s - p->s0 + p->dtpv*(p->outflux*fluxfun(s) + p->influx); +} + +static struct Parameters +get_parameters(struct vdata *vd, const struct cdata *cd, int cell) +{ + struct UnstructuredGrid *g = cd->grid; + struct Parameters p; + + p.s0 = vd->saturation[cell]; + p.influx = cd->source[cell] > 0 ? -cd->source[cell] : 0.0; + p.outflux = cd->source[cell] <= 0 ? -cd->source[cell] : 0.0; + p.dtpv = cd->dt/cd->porevolume[cell]; + + fprintf(stderr, "src:%f pv:%f\n", + cd->source[cell], + cd->porevolume[cell]); + + int i; + vd->saturation[cell] = 0; + for (i=g->cell_facepos[cell]; icell_facepos[cell+1]; ++i) { + int f = g->cell_faces[i]; + double flux; + int other; + + /* Compute cell flux*/ + if (cell == g->face_cells[2*f]) { + flux = cd->darcyflux[f]; + other = g->face_cells[2*f+1]; + } + else { + flux =-cd->darcyflux[f]; + other = g->face_cells[2*f]; + } + + fprintf(stderr, "flux:%f\n", flux); + if (other != -1) { + if (flux < 0.0) { + p.influx += flux*vd->fractionalflow[other]; + } + else { + p.outflux += flux; + } + } + } + return p; +} + + +static enum Method {RIDDER, REGULAFALSI, BISECTION} method = REGULAFALSI; + +void solve(void *vdata, const void *cdata, int cell) +{ + const struct cdata *cd = cdata; + struct vdata *vd = vdata; + + struct Parameters p = get_parameters(vdata, cdata, cell); + + int iterations; + int maxit = 20; + double TOL = 1e-9; + + switch (method) { + default: + case BISECTION: + vd->saturation[cell] = bisection ( &G, &p, TOL, maxit, &iterations); + break; + + case RIDDER: + vd->saturation[cell] = ridder ( &G, &p, TOL, maxit, &iterations); + break; + + case REGULAFALSI: + vd->saturation[cell] = regulafalsi( &G, &p, TOL, maxit, &iterations); + break; + } + + vd->fractionalflow[cell] = fluxfun(vd->saturation[cell]); +} + + + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ + diff --git a/opm/core/transport/reorder/twophase.h b/opm/core/transport/reorder/twophase.h new file mode 100644 index 00000000..82082dae --- /dev/null +++ b/opm/core/transport/reorder/twophase.h @@ -0,0 +1,26 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#ifndef TWOPHASE_H_INCLUDED +#define TWOPHASE_H_INCLUDED + + +struct vdata { + double *saturation; /* one per cell */ + double *fractionalflow; /* one per cell */ +}; + +struct UnstructuredGrid; +struct cdata { + struct UnstructuredGrid *grid; + const double *darcyflux; /* one flux per face in cdata::grid*/ + const double *source; /* one source per cell */ + const double *porevolume; /* one volume per cell */ + double dt; +}; + +void solve(void *vdata, const void *cdata, int cell); +#endif /* TWOPHASE_H_INCLUDED */ + +/* Local Variables: */ +/* c-basic-offset:4 */ +/* End: */ diff --git a/opm/core/transport/reorder/twophasetransport.c b/opm/core/transport/reorder/twophasetransport.c new file mode 100644 index 00000000..a9b0cc24 --- /dev/null +++ b/opm/core/transport/reorder/twophasetransport.c @@ -0,0 +1,54 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#include +#include +#include +#include +#include +#include +#include + +int twophasetransport( + const double *porevolume, + const double *source, + double dt, + struct UnstructuredGrid *grid, + const double *darcyflux, + double *saturation) +{ + int sz, i; + + int ncomponents; + int *sequence; + int *components; + + struct cdata cd = {grid, darcyflux, source, porevolume, dt}; + struct vdata vd; + + /* Compute sequence of single-cell problems */ + sequence = malloc( grid->number_of_cells * sizeof *sequence); + components = malloc(( grid->number_of_cells+1) * sizeof *components); + + compute_sequence(grid, darcyflux, sequence, components, &ncomponents); + assert(ncomponents == grid->number_of_cells); + + + vd.saturation = saturation; + vd.fractionalflow = malloc(grid->number_of_cells * + sizeof *vd.fractionalflow); + for(i=0; inumber_of_cells; ++i) + { + vd.fractionalflow[i] = 0.0; + } + + /* Assume all strong components are single-cell domains. */ + for(i=0; inumber_of_cells; ++i) + { + fprintf(stderr,"i:%d\n", i); + solve(&vd, &cd, sequence[i]); + } + + free(vd.fractionalflow); + free(sequence); + free(components); +} diff --git a/opm/core/transport/reorder/twophasetransport.h b/opm/core/transport/reorder/twophasetransport.h new file mode 100644 index 00000000..28caf23c --- /dev/null +++ b/opm/core/transport/reorder/twophasetransport.h @@ -0,0 +1,16 @@ +/* Copyright 2011 (c) Jostein R. Natvig */ + +#ifndef TWOPHASETRANSPORT_H_INCLUDED +#define TWOPHASETRANSPORT_H_INCLUDED + +struct UnstructuredGrid; +int twophasetransport( + const double *porevolume, + const double *source, + double dt, + struct UnstructuredGrid *grid, + const double *darcyflux, + double *saturation); + +#endif /* TWOPHASETRANSPORT_H_INCLUDED */ +