From 29cdeedde33a6bb9528d54b91b3d288075fcf78c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Fri, 20 Jan 2012 14:07:23 +0100
Subject: [PATCH] Build fluid.c and make it into a quadratic Corey fluid.

---
 Makefile.am                        |   2 +
 opm/core/transport/reorder/fluid.c | 105 +++++------------------------
 2 files changed, 19 insertions(+), 88 deletions(-)

diff --git a/Makefile.am b/Makefile.am
index 6cf45cd0..b21a6ed5 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -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		\
diff --git a/opm/core/transport/reorder/fluid.c b/opm/core/transport/reorder/fluid.c
index 44d11474..d777ff16 100644
--- a/opm/core/transport/reorder/fluid.c
+++ b/opm/core/transport/reorder/fluid.c
@@ -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:                */