OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Build fluid.c and make it into a quadratic Corey fluid.

Browse Source
This commit is contained in:
Atgeirr Flø Rasmussen 2012-01-20 14:07:23 +01:00
parent 68c6936cbd
commit 29cdeedde3
2 changed files with 19 additions and 88 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile.am
View File

@ -80,6 +80,7 @@ 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/reorder/fluid.c \
opm/core/transport/reorder/reordersequence.c \
opm/core/transport/reorder/twophase.c \
opm/core/transport/reorder/nlsolvers.c \
@ -195,6 +196,7 @@ 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/reorder/fluid.h \
opm/core/transport/reorder/nlsolvers.h \
opm/core/transport/reorder/reordersequence.h \
opm/core/transport/reorder/twophase.h \

105
opm/core/transport/reorder/fluid.c
View File

@ -1,46 +1,18 @@
#include <math.h>
#ifdef MATLAB_MEX_FILE
#include <mex.h>
#include "fluid.h"
#else
#include <stdio.h>
#include <fluid.h>
#endif
#include <opm/core/transport/reorder/fluid.h>
static int nr;
static double *viscw;
static double *visco;
static double *nw;
static double *no;
static double *srw;
static double *sro;
/*---------------------------------------------------------------------------*/
double fluxfun(double sw, const int reg)
/*---------------------------------------------------------------------------*/
{
double den, so, mw, mo;
#ifdef MATLAB_MEX_FILE
mxAssert ((0 <= reg) && (reg < nr),"Region number out of bounds.");
#else
if(0 <= reg)
{
fprintf(stderr, "Region number out of bounds.");
exit(EXIT_FAILURE);
}
#endif
den = 1 - srw[reg] - sro[reg];
so = (1 - sw - sro[reg])/den;
sw = (sw - srw[reg])/den;
sw = sw > 1.0 ? 1.0 : sw;
sw = sw < 0.0 ? 0.0 : sw;
mw = pow(sw, nw[reg]) / viscw[reg];
mo = pow(so, no[reg]) / visco[reg];
double so, mw, mo, vw, vo;
/* Hardcoding behaviour to make test program work */
so = 1.0 - sw;
vw = 0.001;
vo = 0.003;
mw = sw*sw / vw;
mo = so*so / vo;
return mw / (mw + mo);
}
@ -48,65 +20,22 @@ double fluxfun(double sw, const int reg)
double dfluxfun(double sw, const int reg)
/*---------------------------------------------------------------------------*/
{
double den, so, mw, mo, dmw, dmo, fw;
double so, mw, mo, vw, vo, dmw, dmo, fw;
/* Hardcoding behaviour to make test program work */
so = 1.0 - sw;
vw = 0.001;
vo = 0.003;
mw = sw*sw / vw;
mo = so*so / vo;
den = 1 - srw[reg] - sro[reg];
so = (1 - sw - sro[reg])/den;
sw = (sw - srw[reg])/den;
sw = sw > 1.0 ? 1.0 : sw;
sw = sw < 0.0 ? 0.0 : sw;
mw = pow(sw, nw[reg]) / viscw[reg];
mo = pow(so, no[reg]) / visco[reg];
dmw = nw[reg] * pow(sw, nw[reg]-1) / viscw[reg];
dmo = no[reg] * pow(so, no[reg]-1) / visco[reg];
dmw = 2.0*sw / vw;
dmo = 2.0*so / vo;
fw = mw / (mw+mo);
return (dmw * (1-fw) - dmo * fw) / (mw+mo);
}
#ifdef MATLAB_MEX_FILE
/* ------------------------------------------------------------------ */
static const mxArray*
getField(const mxArray *a, const char *field)
/* ------------------------------------------------------------------ */
{
int fld_no = mxGetFieldNumber(a, field);
mxAssert(fld_no >= 0, "Missing field in fluid opts.");
return mxGetFieldByNumber(a , 0, fld_no);
}
/*---------------------------------------------------------------------------*/
void init_fluid(const mxArray *arr)
/*---------------------------------------------------------------------------*/
{
viscw = mxGetPr(getField(arr, "viscw"));
visco = mxGetPr(getField(arr, "visco"));
srw = mxGetPr(getField(arr, "srw"));
sro = mxGetPr(getField(arr, "sro"));
nw = mxGetPr(getField(arr, "nw"));
no = mxGetPr(getField(arr, "no"));
nr = mxGetNumberOfElements(getField(arr, "no"));
}
#else
void init_fluid(void)
{
fprintf(stderr, "Not implemented");
}
#endif
/*---------------------------------------------------------------------------*/
int getNumSatRegions(void)
/*---------------------------------------------------------------------------*/
{
return nr;
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */