Merge remote-tracking branch 'hnil/master'

CMakeLists.txt

@@ -2,8 +2,10 @@ cmake_minimum_required (VERSION 2.6)
 project (opm-core)
 
 enable_language(Fortran)
-
+set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR})
-
+SET(CMAKE_BUILD_TYPE "debug")
+SET(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -D_GLIBCXX_DEBUG -pg -Wall -fopenmp -ggdb")
+SET(CMAKE_C_FLAGS_DEBUG "${CMAKE_C_FLAGS_DEBUG} -pg -Wall -fopenmp -ggdb")
 
 set(Boost_USE_STATIC_LIBS        ON)
 set(Boost_USE_MULTITHREADED      ON)
@@ -16,94 +18,50 @@ include_directories(${PROJECT_SOURCE_DIR} ${Boost_INCLUDE_DIRS})
 
 
 # The opmcore library
-add_library(opmcore
+FILE(GLOB_RECURSE C_FILES_CORE "opm/core/*.c")
-opm/core/eclipse/EclipseGridInspector.cpp
+FILE(GLOB_RECURSE CPP_FILES_CORE "opm/core/*.cpp")
-opm/core/eclipse/EclipseGridParser.cpp
+FILE(GLOB_RECURSE REMOVE_FILES   "processgrid.c")
-opm/core/fluid/blackoil/BlackoilPvtProperties.cpp
+FILE(GLOB_RECURSE REMOVE_FILESMX   "mx*.c")
-opm/core/fluid/blackoil/SinglePvtDead.cpp
+FILE(GLOB_RECURSE REMOVE_FILESAGMG   "*AGMG.cpp" "*test*")
-opm/core/fluid/blackoil/SinglePvtLiveGas.cpp
+list(REMOVE_ITEM C_FILES_CORE ${REMOVE_FILES} ${REMOVE_FILESMX} )
-opm/core/fluid/blackoil/SinglePvtLiveOil.cpp
+list(REMOVE_ITEM CPP_FILES_CORE ${REMOVE_FILES} ${REMOVE_FILESAGMG})
-opm/core/fluid/blackoil/SinglePvtInterface.cpp
+add_library(opmcore ${C_FILES_CORE} ${CPP_FILES_CORE} )
-opm/core/fluid/BlackoilPropertiesBasic.cpp
+
-opm/core/fluid/BlackoilPropertiesFromDeck.cpp
-opm/core/fluid/IncompPropertiesBasic.cpp
-opm/core/fluid/IncompPropertiesFromDeck.cpp
-opm/core/fluid/PvtPropertiesBasic.cpp
-opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
-opm/core/fluid/RockBasic.cpp
-opm/core/fluid/RockFromDeck.cpp
-opm/core/fluid/SaturationPropsBasic.cpp
-opm/core/fluid/SaturationPropsFromDeck.cpp
-opm/core/utility/MonotCubicInterpolator.cpp
-opm/core/utility/parameters/Parameter.cpp
-opm/core/utility/parameters/ParameterGroup.cpp
-opm/core/utility/parameters/ParameterTools.cpp
-opm/core/utility/parameters/ParameterXML.cpp
-opm/core/utility/parameters/tinyxml/tinystr.cpp
-opm/core/utility/parameters/tinyxml/tinyxml.cpp
-opm/core/utility/parameters/tinyxml/tinyxmlerror.cpp
-opm/core/utility/parameters/tinyxml/tinyxmlparser.cpp
-opm/core/utility/cart_grid.c
-opm/core/utility/cpgpreprocess/geometry.c
-opm/core/utility/cpgpreprocess/preprocess.c
-opm/core/utility/cpgpreprocess/readvector.cpp
-opm/core/utility/cpgpreprocess/cgridinterface.c
-opm/core/utility/cpgpreprocess/sparsetable.c
-opm/core/utility/cpgpreprocess/facetopology.c
-opm/core/utility/cpgpreprocess/uniquepoints.c
-opm/core/utility/StopWatch.cpp
-opm/core/utility/newwells.c
-opm/core/utility/writeVtkData.cpp
-opm/core/GridManager.cpp
-opm/core/linalg/sparse_sys.c
-opm/core/pressure/cfsh.c
-opm/core/pressure/flow_bc.c
-opm/core/pressure/well.c
-opm/core/pressure/fsh_common_impl.c
-opm/core/pressure/fsh.c
-opm/core/pressure/tpfa/ifs_tpfa.c
-opm/core/pressure/tpfa/compr_source.c
-opm/core/pressure/tpfa/cfs_tpfa.c
-opm/core/pressure/tpfa/compr_bc.c
-opm/core/pressure/tpfa/compr_quant.c
-opm/core/pressure/tpfa/compr_quant_general.c
-opm/core/pressure/tpfa/cfs_tpfa_residual.c
-opm/core/pressure/tpfa/trans_tpfa.c
-opm/core/pressure/msmfem/coarse_conn.c
-opm/core/pressure/msmfem/partition.c
-opm/core/pressure/msmfem/hash_set.c
-opm/core/pressure/msmfem/ifsh_ms.c
-opm/core/pressure/msmfem/dfs.c
-opm/core/pressure/msmfem/coarse_sys.c
-opm/core/pressure/ifsh.c
-opm/core/pressure/IncompTpfa.cpp
-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/reorder/TransportModelInterface.cpp
-opm/core/transport/reorder/TransportModelTwophase.cpp
-opm/core/transport/reorder/reordersequence.cpp
-opm/core/transport/reorder/nlsolvers.c
-opm/core/transport/reorder/tarjan.c
-opm/core/linalg/call_umfpack.c
-opm/core/pressure/IncompTpfa.cpp
-opm/core/linalg/LinearSolverInterface.cpp
-opm/core/linalg/LinearSolverIstl.cpp
-opm/core/linalg/LinearSolverUmfpack.cpp
-)
 
 target_link_libraries(opmcore
   ${UMFPACK_LIBRARIES} ${LAPACK_LINKER_FLAGS} ${LAPACK_LIBRARIES} ${Boost_LIBRARIES}
  -lcholmod -lcamd -lccolamd -lmetis -ldunecommon
 )
 
+FILE(GLOB CPP_EXAMPLES "examples/*.cpp")
+FILE(GLOB CPP_tests "tests/*.cpp")
 add_executable(spu_2p examples/spu_2p.cpp)
+add_executable(sim_2p_incomp_reorder examples/sim_2p_incomp_reorder.cpp)
+add_executable(sim_wateroil examples/sim_wateroil.cpp)
+
+add_executable(pvt_test tests/pvt_test.cpp)
+add_executable(relperm_test tests/relperm_test.cpp)
 
 target_link_libraries(spu_2p
   opmcore
 )
+target_link_libraries(sim_2p_incomp_reorder
+  opmcore
+)
 
-set_target_properties(opmcore spu_2p PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
+target_link_libraries(sim_wateroil
+  opmcore
+)
+
+target_link_libraries(pvt_test
+  opmcore
+)
+target_link_libraries(relperm_test
+  opmcore
+)
+
+#set_target_properties(opmcore spu_2p PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
+#set_target_properties(opmcore sim_2p_incomp_reorder PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
+#set_target_properties(opmcore sim_wateroil PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
+#set_target_properties(opmcore pvt_test PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)
+set_target_properties(opmcore relperm_test PROPERTIES COMPILE_FLAGS -m64 LINKER_LANGUAGE CXX LINK_FLAGS -m64)

FindUmfPack.cmake

@@ -0,0 +1,56 @@
+if (UMFPACK_INCLUDES AND UMFPACK_LIBRARIES)
+  set(UMFPACK_FIND_QUIETLY TRUE)
+endif (UMFPACK_INCLUDES AND UMFPACK_LIBRARIES)
+
+find_package(BLAS)
+
+if(BLAS_FOUND)
+
+  find_path(UMFPACK_INCLUDES
+    NAMES
+    umfpack.h
+    PATHS
+    $ENV{UMFPACKDIR}
+    ${INCLUDE_INSTALL_DIR}
+    PATH_SUFFIXES
+    suitesparse
+  )
+
+  find_library(UMFPACK_LIBRARIES umfpack PATHS $ENV{UMFPACKDIR} ${LIB_INSTALL_DIR})
+
+  if(UMFPACK_LIBRARIES)
+
+    get_filename_component(UMFPACK_LIBDIR ${UMFPACK_LIBRARIES} PATH)
+
+    find_library(AMD_LIBRARY amd PATHS ${UMFPACK_LIBDIR} $ENV{UMFPACKDIR} ${LIB_INSTALL_DIR})
+    if (AMD_LIBRARY)
+      set(UMFPACK_LIBRARIES ${UMFPACK_LIBRARIES} ${AMD_LIBRARY})
+    #else (AMD_LIBRARY)
+    #  set(UMFPACK_LIBRARIES FALSE)
+    endif (AMD_LIBRARY)
+
+  endif(UMFPACK_LIBRARIES)
+
+  if(UMFPACK_LIBRARIES)
+
+    find_library(COLAMD_LIBRARY colamd PATHS ${UMFPACK_LIBDIR} $ENV{UMFPACKDIR} ${LIB_INSTALL_DIR})
+    if (COLAMD_LIBRARY)
+      set(UMFPACK_LIBRARIES ${UMFPACK_LIBRARIES} ${COLAMD_LIBRARY})
+    #else (COLAMD_LIBRARY)
+    #  set(UMFPACK_LIBRARIES FALSE)
+    endif (COLAMD_LIBRARY)
+
+  endif(UMFPACK_LIBRARIES)
+
+  if(UMFPACK_LIBRARIES)
+    set(UMFPACK_LIBRARIES ${UMFPACK_LIBRARIES} ${BLAS_LIBRARIES})
+  endif(UMFPACK_LIBRARIES)
+
+endif(BLAS_FOUND)
+
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(UMFPACK DEFAULT_MSG
+                                  UMFPACK_INCLUDES UMFPACK_LIBRARIES)
+
+mark_as_advanced(UMFPACK_INCLUDES UMFPACK_LIBRARIES AMD_LIBRARY COLAMD_LIBRARY)

Makefile.am

@@ -50,6 +50,8 @@ opm/core/fluid/RockCompressibility.cpp			\
 opm/core/fluid/RockFromDeck.cpp				\
 opm/core/fluid/SaturationPropsBasic.cpp			\
 opm/core/fluid/SaturationPropsFromDeck.cpp		\
+opm/core/fluid/SatFuncStone2.cpp			\
+opm/core/fluid/SatFuncSimple.cpp                        \
 opm/core/fluid/blackoil/BlackoilPvtProperties.cpp	\
 opm/core/fluid/blackoil/SinglePvtDead.cpp		\
 opm/core/fluid/blackoil/SinglePvtInterface.cpp		\
@@ -147,6 +149,8 @@ opm/core/fluid/RockCompressibility.hpp			\
 opm/core/fluid/RockFromDeck.hpp				\
 opm/core/fluid/SaturationPropsBasic.hpp			\
 opm/core/fluid/SaturationPropsFromDeck.hpp		\
+opm/core/fluid/SatFuncStone2.hpp			\
+opm/core/fluid/SatFuncSimple.hpp                        \
 opm/core/fluid/SimpleFluid2p.hpp			\
 opm/core/fluid/blackoil/BlackoilPhases.hpp		\
 opm/core/fluid/blackoil/BlackoilPvtProperties.hpp	\

opm/core/fluid/BlackoilPropertiesFromDeck.cpp

@@ -18,7 +18,7 @@
 */
 
 #include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
-
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 namespace Opm
 {
 
@@ -34,6 +34,21 @@ namespace Opm
 	}
     }
 
+    BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
+                                                           const UnstructuredGrid& grid,
+                                                           const parameter::ParameterGroup& param)
+    {
+        rock_.init(deck, grid);
+        int samples = param.getDefault("dead_tab_size", 1025);
+        pvt_.init(deck, samples);
+        satprops_.init(deck, grid, param);
+
+        if (pvt_.numPhases() != satprops_.numPhases()) {
+            THROW("BlackoilPropertiesBasic::BlackoilPropertiesBasic() - Inconsistent number of phases in pvt data ("
+                  << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_.numPhases() << ").");
+        }
+    }
+
     BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()
     {
     }

opm/core/fluid/BlackoilPropertiesFromDeck.hpp

@@ -26,6 +26,7 @@
 #include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp>
 #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 
 struct UnstructuredGrid;
 
@@ -43,8 +44,10 @@ namespace Opm
         ///                      mapping from cell indices (typically from a processed grid)
         ///                      to logical cartesian indices consistent with the deck.
         BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
-                                   const UnstructuredGrid& grid);
+                                   const UnstructuredGrid& grid);        
-
+        BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
+                                   const UnstructuredGrid& grid,
+                                   const parameter::ParameterGroup& param);
         /// Destructor.
         virtual ~BlackoilPropertiesFromDeck();
 

opm/core/fluid/SatFuncSimple.cpp

@@ -0,0 +1,213 @@
+#include <opm/core/fluid/SatFuncSimple.hpp>
+#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
+#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
+#include <opm/core/utility/buildUniformMonotoneTable.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <iostream>
+namespace Opm
+{
+
+
+
+
+    // ----------- Methods of SatFuncSet below -----------
+    void SatFuncSimple::init(const EclipseGridParser& deck,
+                             const int table_num,
+                             const PhaseUsage phase_usg){
+        init(deck, table_num, phase_usg, 200);
+    }
+
+    void SatFuncSimple::init(const EclipseGridParser& deck,
+                             const int table_num,
+                             const PhaseUsage phase_usg,
+                             const int samples)
+    {
+        phase_usage = phase_usg;
+        double swco = 0.0;
+        double swmax = 1.0;
+        if (phase_usage.phase_used[Aqua]) {
+            const SWOF::table_t& swof_table = deck.getSWOF().swof_;
+            const std::vector<double>& sw = swof_table[table_num][0];
+            const std::vector<double>& krw = swof_table[table_num][1];
+            const std::vector<double>& krow = swof_table[table_num][2];
+            const std::vector<double>& pcow = swof_table[table_num][3];
+            buildUniformMonotoneTable(sw, krw,  samples, krw_);
+            buildUniformMonotoneTable(sw, krow, samples, krow_);
+            buildUniformMonotoneTable(sw, pcow, samples, pcow_);
+            krocw_ = krow[0]; // At connate water -> ecl. SWOF
+            swco = sw[0];
+            smin_[phase_usage.phase_pos[Aqua]] = sw[0];
+            swmax = sw.back();
+            smax_[phase_usage.phase_pos[Aqua]] = sw.back();
+        }
+        if (phase_usage.phase_used[Vapour]) {
+            const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
+            const std::vector<double>& sg = sgof_table[table_num][0];
+            const std::vector<double>& krg = sgof_table[table_num][1];
+            const std::vector<double>& krog = sgof_table[table_num][2];
+            const std::vector<double>& pcog = sgof_table[table_num][3];
+            buildUniformMonotoneTable(sg, krg,  samples, krg_);
+            buildUniformMonotoneTable(sg, krog, samples, krog_);
+            buildUniformMonotoneTable(sg, pcog, samples, pcog_);
+            smin_[phase_usage.phase_pos[Vapour]] = sg[0];
+            if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
+                THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
+                      ", should equal (1.0 - connate water sat) = " << (1.0 - swco));
+            }
+            smax_[phase_usage.phase_pos[Vapour]] = sg.back();
+        }
+        // These only consider water min/max sats. Consider gas sats?
+        smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
+        smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
+    }
+
+
+    void SatFuncSimple::evalKr(const double* s, double* kr) const
+    {
+        if (phase_usage.num_phases == 3) {
+            // Stone-II relative permeability model.
+            double sw = s[Aqua];
+            double sg = s[Vapour];
+            double krw = krw_(sw);
+            double krg = krg_(sg);
+            double krow = krow_(sw + sg); // = 1 - so
+            double krog = krog_(sg);      // = 1 - so - sw
+            double krocw = krocw_;
+            kr[Aqua] = krw;
+            kr[Vapour] = krg;
+            kr[Liquid] = krow;
+            if (kr[Liquid] < 0.0) {
+                kr[Liquid] = 0.0;
+            }
+            return;
+        }
+        // We have a two-phase situation. We know that oil is active.
+        if (phase_usage.phase_used[Aqua]) {
+            int wpos = phase_usage.phase_pos[Aqua];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sw = s[wpos];
+            double krw = krw_(sw);
+            double krow = krow_(sw);
+            kr[wpos] = krw;
+            kr[opos] = krow;
+        } else {
+            ASSERT(phase_usage.phase_used[Vapour]);
+            int gpos = phase_usage.phase_pos[Vapour];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sg = s[gpos];
+            double krg = krg_(sg);
+            double krog = krog_(sg);
+            kr[gpos] = krg;
+            kr[opos] = krog;
+        }
+    }
+
+
+    void SatFuncSimple::evalKrDeriv(const double* s, double* kr, double* dkrds) const
+    {
+        const int np = phase_usage.num_phases;
+        std::fill(dkrds, dkrds + np*np, 0.0);
+
+        if (np == 3) {
+            // Stone-II relative permeability model.
+            double sw = s[Aqua];
+            double sg = s[Vapour];
+            double krw = krw_(sw);
+            double dkrww = krw_.derivative(sw);
+            double krg = krg_(sg);
+            double dkrgg = krg_.derivative(sg);
+            double krow = krow_(sw + sg);
+            double dkrow = krow_.derivative(sw + sg);
+            double krog = krog_(sg);
+            double dkrog = krog_.derivative(sg);
+            double krocw = krocw_;
+            kr[Aqua] = krw;
+            kr[Vapour] = krg;
+            kr[Liquid] =  krow;
+            //krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
+            if (kr[Liquid] < 0.0) {
+                kr[Liquid] = 0.0;
+            }
+            dkrds[Aqua + Aqua*np] = dkrww;
+            dkrds[Vapour + Vapour*np] = dkrgg;
+            //dkrds[Liquid + Aqua*np] = dkrow;
+            dkrds[Liquid + Liquid*np] = -dkrow;
+                    //krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
+            dkrds[Liquid + Vapour*np] = 0.0;
+                    //krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
+                    //+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
+            return;
+        }
+        // We have a two-phase situation. We know that oil is active.
+        if (phase_usage.phase_used[Aqua]) {
+            int wpos = phase_usage.phase_pos[Aqua];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sw = s[wpos];
+            double krw = krw_(sw);
+            double dkrww = krw_.derivative(sw);
+            double krow = krow_(sw);
+            double dkrow = krow_.derivative(sw);
+            kr[wpos] = krw;
+            kr[opos] = krow;
+            dkrds[wpos + wpos*np] = dkrww;
+            dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
+        } else {
+            ASSERT(phase_usage.phase_used[Vapour]);
+            int gpos = phase_usage.phase_pos[Vapour];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sg = s[gpos];
+            double krg = krg_(sg);
+            double dkrgg = krg_.derivative(sg);
+            double krog = krog_(sg);
+            double dkrog = krog_.derivative(sg);
+            kr[gpos] = krg;
+            kr[opos] = krog;
+            dkrds[gpos + gpos*np] = dkrgg;
+            dkrds[opos + gpos*np] = dkrog;
+        }
+
+    }
+
+
+    void SatFuncSimple::evalPc(const double* s, double* pc) const
+    {
+        pc[phase_usage.phase_pos[Liquid]] = 0.0;
+        if (phase_usage.phase_used[Aqua]) {
+            int pos = phase_usage.phase_pos[Aqua];
+            pc[pos] = pcow_(s[pos]);
+        }
+        if (phase_usage.phase_used[Vapour]) {
+            int pos = phase_usage.phase_pos[Vapour];
+            pc[pos] = pcog_(s[pos]);
+        }
+    }
+
+    void SatFuncSimple::evalPcDeriv(const double* s, double* pc, double* dpcds) const
+    {
+        // The problem of determining three-phase capillary pressures
+        // is very hard experimentally, usually one extends two-phase
+        // data (as for relative permeability).
+        // In our approach the derivative matrix is quite sparse, only
+        // the diagonal elements corresponding to non-oil phases are
+        // (potentially) nonzero.
+        const int np = phase_usage.num_phases;
+        std::fill(dpcds, dpcds + np*np, 0.0);
+        pc[phase_usage.phase_pos[Liquid]] = 0.0;
+        if (phase_usage.phase_used[Aqua]) {
+            int pos = phase_usage.phase_pos[Aqua];
+            pc[pos] = pcow_(s[pos]);
+            dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
+        }
+        if (phase_usage.phase_used[Vapour]) {
+            int pos = phase_usage.phase_pos[Vapour];
+            pc[pos] = pcog_(s[pos]);
+            dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
+        }
+    }
+
+
+
+
+} // namespace Opm

opm/core/fluid/SatFuncSimple.hpp

@@ -0,0 +1,32 @@
+#ifndef SATFUNCSIMPLE_HPP
+#define SATFUNCSIMPLE_HPP
+#include <opm/core/eclipse/EclipseGridParser.hpp>
+#include <opm/core/utility/UniformTableLinear.hpp>
+#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
+#include <vector>
+namespace Opm
+{
+    class SatFuncSimple: public BlackoilPhases
+    {
+    public:
+        void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg);
+        void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg,
+                  const int samples);
+        void evalKr(const double* s, double* kr) const;
+        void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
+        void evalPc(const double* s, double* pc) const;
+        void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
+        double smin_[PhaseUsage::MaxNumPhases];
+        double smax_[PhaseUsage::MaxNumPhases];
+    private:
+        PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
+        UniformTableLinear<double> krw_;
+        UniformTableLinear<double> krow_;
+        UniformTableLinear<double> pcow_;
+        UniformTableLinear<double> krg_;
+        UniformTableLinear<double> krog_;
+        UniformTableLinear<double> pcog_;
+        double krocw_; // = krow_(s_wc)
+    };
+} // namespace Opm
+#endif // SATFUNCSIMPLE_HPP

opm/core/fluid/SatFuncStone2.cpp

@@ -0,0 +1,209 @@
+#include <opm/core/fluid/SatFuncStone2.hpp>
+#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
+#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
+#include <opm/core/grid.h>
+#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
+#include <opm/core/utility/buildUniformMonotoneTable.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
+#include <iostream>
+namespace Opm
+{
+
+
+
+
+    // ----------- Methods of SatFuncSet below -----------
+    void SatFuncStone2::init(const EclipseGridParser& deck,
+                             const int table_num,
+                             const PhaseUsage phase_usg){
+        init(deck, table_num, phase_usg, 200);
+    }
+
+    void SatFuncStone2::init(const EclipseGridParser& deck,
+                             const int table_num,
+                             const PhaseUsage phase_usg,
+                             const int samples)
+    {
+        phase_usage = phase_usg;
+        double swco = 0.0;
+        double swmax = 1.0;
+        if (phase_usage.phase_used[Aqua]) {
+            const SWOF::table_t& swof_table = deck.getSWOF().swof_;
+            const std::vector<double>& sw = swof_table[table_num][0];
+            const std::vector<double>& krw = swof_table[table_num][1];
+            const std::vector<double>& krow = swof_table[table_num][2];
+            const std::vector<double>& pcow = swof_table[table_num][3];
+            buildUniformMonotoneTable(sw, krw,  samples, krw_);
+            buildUniformMonotoneTable(sw, krow, samples, krow_);
+            buildUniformMonotoneTable(sw, pcow, samples, pcow_);
+            krocw_ = krow[0]; // At connate water -> ecl. SWOF
+            swco = sw[0];
+            smin_[phase_usage.phase_pos[Aqua]] = sw[0];
+            swmax = sw.back();
+            smax_[phase_usage.phase_pos[Aqua]] = sw.back();
+        }
+        if (phase_usage.phase_used[Vapour]) {
+            const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
+            const std::vector<double>& sg = sgof_table[table_num][0];
+            const std::vector<double>& krg = sgof_table[table_num][1];
+            const std::vector<double>& krog = sgof_table[table_num][2];
+            const std::vector<double>& pcog = sgof_table[table_num][3];
+            buildUniformMonotoneTable(sg, krg,  samples, krg_);
+            buildUniformMonotoneTable(sg, krog, samples, krog_);
+            buildUniformMonotoneTable(sg, pcog, samples, pcog_);
+            smin_[phase_usage.phase_pos[Vapour]] = sg[0];
+            if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
+                THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
+                      ", should equal (1.0 - connate water sat) = " << (1.0 - swco));
+            }
+            smax_[phase_usage.phase_pos[Vapour]] = sg.back();
+        }
+        // These only consider water min/max sats. Consider gas sats?
+        smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
+        smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
+    }
+
+
+    void SatFuncStone2::evalKr(const double* s, double* kr) const
+    {
+        if (phase_usage.num_phases == 3) {
+            // Stone-II relative permeability model.
+            double sw = s[Aqua];
+            double sg = s[Vapour];
+            double krw = krw_(sw);
+            double krg = krg_(sg);
+            double krow = krow_(sw + sg); // = 1 - so
+            double krog = krog_(sg);      // = 1 - so - sw
+            double krocw = krocw_;
+            kr[Aqua] = krw;
+            kr[Vapour] = krg;
+            kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
+            if (kr[Liquid] < 0.0) {
+                kr[Liquid] = 0.0;
+            }
+            return;
+        }
+        // We have a two-phase situation. We know that oil is active.
+        if (phase_usage.phase_used[Aqua]) {
+            int wpos = phase_usage.phase_pos[Aqua];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sw = s[wpos];
+            double krw = krw_(sw);
+            double krow = krow_(sw);
+            kr[wpos] = krw;
+            kr[opos] = krow;
+        } else {
+            ASSERT(phase_usage.phase_used[Vapour]);
+            int gpos = phase_usage.phase_pos[Vapour];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sg = s[gpos];
+            double krg = krg_(sg);
+            double krog = krog_(sg);
+            kr[gpos] = krg;
+            kr[opos] = krog;
+        }
+    }
+
+
+    void SatFuncStone2::evalKrDeriv(const double* s, double* kr, double* dkrds) const
+    {
+        const int np = phase_usage.num_phases;
+        std::fill(dkrds, dkrds + np*np, 0.0);
+
+        if (np == 3) {
+            // Stone-II relative permeability model.
+            double sw = s[Aqua];
+            double sg = s[Vapour];
+            double krw = krw_(sw);
+            double dkrww = krw_.derivative(sw);
+            double krg = krg_(sg);
+            double dkrgg = krg_.derivative(sg);
+            double krow = krow_(sw + sg);
+            double dkrow = krow_.derivative(sw + sg);
+            double krog = krog_(sg);
+            double dkrog = krog_.derivative(sg);
+            double krocw = krocw_;
+            kr[Aqua] = krw;
+            kr[Vapour] = krg;
+            kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
+            if (kr[Liquid] < 0.0) {
+                kr[Liquid] = 0.0;
+            }
+            dkrds[Aqua + Aqua*np] = dkrww;
+            dkrds[Vapour + Vapour*np] = dkrgg;
+            dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
+            dkrds[Liquid + Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
+                    + krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
+            return;
+        }
+        // We have a two-phase situation. We know that oil is active.
+        if (phase_usage.phase_used[Aqua]) {
+            int wpos = phase_usage.phase_pos[Aqua];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sw = s[wpos];
+            double krw = krw_(sw);
+            double dkrww = krw_.derivative(sw);
+            double krow = krow_(sw);
+            double dkrow = krow_.derivative(sw);
+            kr[wpos] = krw;
+            kr[opos] = krow;
+            dkrds[wpos + wpos*np] = dkrww;
+            dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
+        } else {
+            ASSERT(phase_usage.phase_used[Vapour]);
+            int gpos = phase_usage.phase_pos[Vapour];
+            int opos = phase_usage.phase_pos[Liquid];
+            double sg = s[gpos];
+            double krg = krg_(sg);
+            double dkrgg = krg_.derivative(sg);
+            double krog = krog_(sg);
+            double dkrog = krog_.derivative(sg);
+            kr[gpos] = krg;
+            kr[opos] = krog;
+            dkrds[gpos + gpos*np] = dkrgg;
+            dkrds[opos + gpos*np] = dkrog;
+        }
+
+    }
+
+
+    void SatFuncStone2::evalPc(const double* s, double* pc) const
+    {
+        pc[phase_usage.phase_pos[Liquid]] = 0.0;
+        if (phase_usage.phase_used[Aqua]) {
+            int pos = phase_usage.phase_pos[Aqua];
+            pc[pos] = pcow_(s[pos]);
+        }
+        if (phase_usage.phase_used[Vapour]) {
+            int pos = phase_usage.phase_pos[Vapour];
+            pc[pos] = pcog_(s[pos]);
+        }
+    }
+
+    void SatFuncStone2::evalPcDeriv(const double* s, double* pc, double* dpcds) const
+    {
+        // The problem of determining three-phase capillary pressures
+        // is very hard experimentally, usually one extends two-phase
+        // data (as for relative permeability).
+        // In our approach the derivative matrix is quite sparse, only
+        // the diagonal elements corresponding to non-oil phases are
+        // (potentially) nonzero.
+        const int np = phase_usage.num_phases;
+        std::fill(dpcds, dpcds + np*np, 0.0);
+        pc[phase_usage.phase_pos[Liquid]] = 0.0;
+        if (phase_usage.phase_used[Aqua]) {
+            int pos = phase_usage.phase_pos[Aqua];
+            pc[pos] = pcow_(s[pos]);
+            dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
+        }
+        if (phase_usage.phase_used[Vapour]) {
+            int pos = phase_usage.phase_pos[Vapour];
+            pc[pos] = pcog_(s[pos]);
+            dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
+        }
+    }
+
+
+
+
+} // namespace Opm

opm/core/fluid/SatFuncStone2.hpp

@@ -0,0 +1,50 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#ifndef SATFUNCSTONE2_HPP
+#define SATFUNCSTONE2_HPP
+#include <opm/core/eclipse/EclipseGridParser.hpp>
+#include <opm/core/utility/UniformTableLinear.hpp>
+#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
+#include <vector>
+namespace Opm
+{
+    class SatFuncStone2: public BlackoilPhases
+    {
+    public:
+        void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg);
+        void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg,
+                  const int samples);
+        void evalKr(const double* s, double* kr) const;
+        void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
+        void evalPc(const double* s, double* pc) const;
+        void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
+        double smin_[PhaseUsage::MaxNumPhases];
+        double smax_[PhaseUsage::MaxNumPhases];
+    private:
+        PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
+        UniformTableLinear<double> krw_;
+        UniformTableLinear<double> krow_;
+        UniformTableLinear<double> pcow_;
+        UniformTableLinear<double> krg_;
+        UniformTableLinear<double> krog_;
+        UniformTableLinear<double> pcog_;
+        double krocw_; // = krow_(s_wc)
+    };
+} // namespace Opm
+#endif // SATFUNCSTONE2_HPP

opm/core/fluid/SaturationPropsFromDeck.cpp

@@ -88,7 +88,62 @@ namespace Opm
             satfuncset_[table].init(deck, table, phase_usage_);
         }
     }
+    void SaturationPropsFromDeck::init(const EclipseGridParser& deck,
+              const UnstructuredGrid& grid,
+              const parameter::ParameterGroup& param){
+        phase_usage_ = phaseUsageFromDeck(deck);
 
+        // Extract input data.
+        // Oil phase should be active.
+        if (!phase_usage_.phase_used[Liquid]) {
+            THROW("SaturationPropsFromDeck::init()   --  oil phase must be active.");
+        }
+
+        // Obtain SATNUM, if it exists, and create cell_to_func_.
+        // Otherwise, let the cell_to_func_ mapping be just empty.
+        int satfuncs_expected = 1;
+        if (deck.hasField("SATNUM")) {
+            const std::vector<int>& satnum = deck.getIntegerValue("SATNUM");
+            satfuncs_expected = *std::max_element(satnum.begin(), satnum.end());
+            const int num_cells = grid.number_of_cells;
+            cell_to_func_.resize(num_cells);
+            const int* gc = grid.global_cell;
+            for (int cell = 0; cell < num_cells; ++cell) {
+                const int deck_pos = (gc == NULL) ? cell : gc[cell];
+                cell_to_func_[cell] = satnum[deck_pos] - 1;
+            }
+        }
+
+        // Find number of tables, check for consistency.
+        enum { Uninitialized = -1 };
+        int num_tables = Uninitialized;
+        if (phase_usage_.phase_used[Aqua]) {
+            const SWOF::table_t& swof_table = deck.getSWOF().swof_;
+            num_tables = swof_table.size();
+            if (num_tables < satfuncs_expected) {
+                THROW("Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected);
+            }
+        }
+        if (phase_usage_.phase_used[Vapour]) {
+            const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
+            int num_sgof_tables = sgof_table.size();
+            if (num_sgof_tables < satfuncs_expected) {
+                THROW("Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected);
+            }
+            if (num_tables == Uninitialized) {
+                num_tables = num_sgof_tables;
+            } else if (num_tables != num_sgof_tables) {
+                THROW("Inconsistent number of tables in SWOF and SGOF.");
+            }
+        }
+
+        // Initialize tables.
+        int tab_size=param.getDefault("tab_size_kr",200);
+        satfuncset_.resize(num_tables);
+        for (int table = 0; table < num_tables; ++table) {
+            satfuncset_[table].init(deck, table, phase_usage_,tab_size);
+        }
+    }
 
 
 
@@ -192,206 +247,12 @@ namespace Opm
 	}
     }
 
-
     // Map the cell number to the correct function set.
-    const SaturationPropsFromDeck::SatFuncSet&
+    const SaturationPropsFromDeck::satfunc_t&
     SaturationPropsFromDeck::funcForCell(const int cell) const
     {
         return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
     }
-
-
-
-    // ----------- Methods of SatFuncSet below -----------
-
-
-    void SaturationPropsFromDeck::SatFuncSet::init(const EclipseGridParser& deck,
-                                                   const int table_num,
-                                                   const PhaseUsage phase_usg)
-    {
-        phase_usage = phase_usg;
-        const int samples = 200;
-        double swco = 0.0;
-        double swmax = 1.0;
-        if (phase_usage.phase_used[Aqua]) {
-            const SWOF::table_t& swof_table = deck.getSWOF().swof_;
-            const std::vector<double>& sw = swof_table[table_num][0];
-            const std::vector<double>& krw = swof_table[table_num][1];
-            const std::vector<double>& krow = swof_table[table_num][2];
-            const std::vector<double>& pcow = swof_table[table_num][3];
-            buildUniformMonotoneTable(sw, krw,  samples, krw_);
-            buildUniformMonotoneTable(sw, krow, samples, krow_);
-            buildUniformMonotoneTable(sw, pcow, samples, pcow_);
-            krocw_ = krow[0]; // At connate water -> ecl. SWOF
-            swco = sw[0];
-            smin_[phase_usage.phase_pos[Aqua]] = sw[0];
-            swmax = sw.back();
-            smax_[phase_usage.phase_pos[Aqua]] = sw.back();
-        }
-        if (phase_usage.phase_used[Vapour]) {
-            const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
-            const std::vector<double>& sg = sgof_table[table_num][0];
-            const std::vector<double>& krg = sgof_table[table_num][1];
-            const std::vector<double>& krog = sgof_table[table_num][2];
-            const std::vector<double>& pcog = sgof_table[table_num][3];
-            buildUniformMonotoneTable(sg, krg,  samples, krg_);
-            buildUniformMonotoneTable(sg, krog, samples, krog_);
-            buildUniformMonotoneTable(sg, pcog, samples, pcog_);
-            smin_[phase_usage.phase_pos[Vapour]] = sg[0];
-            if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
-                THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
-                      ", should equal (1.0 - connate water sat) = " << (1.0 - swco));
-            }
-            smax_[phase_usage.phase_pos[Vapour]] = sg.back();
-        }
-        // These only consider water min/max sats. Consider gas sats?
-        smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
-        smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
-    }
-
-
-    void SaturationPropsFromDeck::SatFuncSet::evalKr(const double* s, double* kr) const
-    {
-        if (phase_usage.num_phases == 3) {
-            // Stone-II relative permeability model.
-            double sw = s[Aqua];
-            double sg = s[Vapour];
-            double krw = krw_(sw);
-            double krg = krg_(sg);
-            double krow = krow_(sw + sg); // = 1 - so
-            double krog = krog_(sg);      // = 1 - so - sw
-            double krocw = krocw_;
-            kr[Aqua] = krw;
-            kr[Vapour] = krg;
-            kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
-            if (kr[Liquid] < 0.0) {
-                kr[Liquid] = 0.0;
-            }
-            return;
-        }
-        // We have a two-phase situation. We know that oil is active.
-        if (phase_usage.phase_used[Aqua]) {
-            int wpos = phase_usage.phase_pos[Aqua];
-            int opos = phase_usage.phase_pos[Liquid];
-            double sw = s[wpos];
-            double krw = krw_(sw);
-            double krow = krow_(sw);
-            kr[wpos] = krw;
-            kr[opos] = krow;
-        } else {
-            ASSERT(phase_usage.phase_used[Vapour]);
-            int gpos = phase_usage.phase_pos[Vapour];
-            int opos = phase_usage.phase_pos[Liquid];
-            double sg = s[gpos];
-            double krg = krg_(sg);
-            double krog = krog_(sg);
-            kr[gpos] = krg;
-            kr[opos] = krog;
-        }
-    }
-
-
-    void SaturationPropsFromDeck::SatFuncSet::evalKrDeriv(const double* s, double* kr, double* dkrds) const
-    {
-        const int np = phase_usage.num_phases;
-        std::fill(dkrds, dkrds + np*np, 0.0);
-
-        if (np == 3) {
-            // Stone-II relative permeability model.
-            double sw = s[Aqua];
-            double sg = s[Vapour];
-            double krw = krw_(sw);
-            double dkrww = krw_.derivative(sw);
-            double krg = krg_(sg);
-            double dkrgg = krg_.derivative(sg);
-            double krow = krow_(sw + sg);
-            double dkrow = krow_.derivative(sw + sg);
-            double krog = krog_(sg);
-            double dkrog = krog_.derivative(sg);
-            double krocw = krocw_;
-            kr[Aqua] = krw;
-            kr[Vapour] = krg;
-            kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
-            if (kr[Liquid] < 0.0) {
-                kr[Liquid] = 0.0;
-            }
-            dkrds[Aqua + Aqua*np] = dkrww;
-            dkrds[Vapour + Vapour*np] = dkrgg;
-            dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
-            dkrds[Liquid + Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
-                + krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
-            return;
-        }
-        // We have a two-phase situation. We know that oil is active.
-        if (phase_usage.phase_used[Aqua]) {
-            int wpos = phase_usage.phase_pos[Aqua];
-            int opos = phase_usage.phase_pos[Liquid];
-            double sw = s[wpos];
-            double krw = krw_(sw);
-            double dkrww = krw_.derivative(sw);
-            double krow = krow_(sw);
-            double dkrow = krow_.derivative(sw);
-            kr[wpos] = krw;
-            kr[opos] = krow;
-            dkrds[wpos + wpos*np] = dkrww;
-            dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
-        } else {
-            ASSERT(phase_usage.phase_used[Vapour]);
-            int gpos = phase_usage.phase_pos[Vapour];
-            int opos = phase_usage.phase_pos[Liquid];
-            double sg = s[gpos];
-            double krg = krg_(sg);
-            double dkrgg = krg_.derivative(sg);
-            double krog = krog_(sg);
-            double dkrog = krog_.derivative(sg);
-            kr[gpos] = krg;
-            kr[opos] = krog;
-            dkrds[gpos + gpos*np] = dkrgg;
-            dkrds[opos + gpos*np] = dkrog;
-        }
-
-    }
-
-
-    void SaturationPropsFromDeck::SatFuncSet::evalPc(const double* s, double* pc) const
-    {
-        pc[phase_usage.phase_pos[Liquid]] = 0.0;
-        if (phase_usage.phase_used[Aqua]) {
-            int pos = phase_usage.phase_pos[Aqua];
-            pc[pos] = pcow_(s[pos]);
-        }
-        if (phase_usage.phase_used[Vapour]) {
-            int pos = phase_usage.phase_pos[Vapour];
-            pc[pos] = pcog_(s[pos]);
-        }
-    }
-
-    void SaturationPropsFromDeck::SatFuncSet::evalPcDeriv(const double* s, double* pc, double* dpcds) const
-    {
-	// The problem of determining three-phase capillary pressures
-	// is very hard experimentally, usually one extends two-phase
-	// data (as for relative permeability).
-	// In our approach the derivative matrix is quite sparse, only
-	// the diagonal elements corresponding to non-oil phases are
-	// (potentially) nonzero.
-        const int np = phase_usage.num_phases;
-	std::fill(dpcds, dpcds + np*np, 0.0);
-        pc[phase_usage.phase_pos[Liquid]] = 0.0;
-        if (phase_usage.phase_used[Aqua]) {
-            int pos = phase_usage.phase_pos[Aqua];
-            pc[pos] = pcow_(s[pos]);
-            dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
-        }
-        if (phase_usage.phase_used[Vapour]) {
-            int pos = phase_usage.phase_pos[Vapour];
-            pc[pos] = pcog_(s[pos]);
-            dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
-        }
-    }
-
-
-
-
 } // namespace Opm
 
 

opm/core/fluid/SaturationPropsFromDeck.hpp

@@ -19,10 +19,12 @@
 
 #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
 #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
-
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
 #include <opm/core/utility/UniformTableLinear.hpp>
 #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
+#include <opm/core/fluid/SatFuncStone2.hpp>
+#include <opm/core/fluid/SatFuncSimple.hpp>
 #include <vector>
 
 struct UnstructuredGrid;
@@ -44,6 +46,10 @@ namespace Opm
         void init(const EclipseGridParser& deck,
                   const UnstructuredGrid& grid);
 
+        void init(const EclipseGridParser& deck,
+                  const UnstructuredGrid& grid,
+                  const parameter::ParameterGroup& param);
+
         /// \return   P, the number of phases.
         int numPhases() const;
 
@@ -88,30 +94,11 @@ namespace Opm
 
     private:
         PhaseUsage phase_usage_;
-        class SatFuncSet
+        typedef SatFuncSimple satfunc_t;
-        {
+        std::vector<satfunc_t> satfuncset_;
-        public:
-            void init(const EclipseGridParser& deck, const int table_num, PhaseUsage phase_usg);
-            void evalKr(const double* s, double* kr) const;
-            void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
-            void evalPc(const double* s, double* pc) const;
-            void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
-            double smin_[PhaseUsage::MaxNumPhases];
-            double smax_[PhaseUsage::MaxNumPhases];
-        private:
-            PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
-            UniformTableLinear<double> krw_;
-            UniformTableLinear<double> krow_;
-            UniformTableLinear<double> pcow_;
-            UniformTableLinear<double> krg_;
-            UniformTableLinear<double> krog_;
-            UniformTableLinear<double> pcog_;
-            double krocw_; // = krow_(s_wc)
-        };
-        std::vector<SatFuncSet> satfuncset_;
         std::vector<int> cell_to_func_; // = SATNUM - 1
 
-        const SatFuncSet& funcForCell(const int cell) const;
+        const satfunc_t& funcForCell(const int cell) const;
     };
 
 

opm/core/fluid/blackoil/BlackoilPvtProperties.cpp

@@ -39,7 +39,7 @@ namespace Opm
     }
 
 
-    void BlackoilPvtProperties::init(const EclipseGridParser& deck)
+    void BlackoilPvtProperties::init(const EclipseGridParser& deck,const int samples)
     {
         typedef std::vector<std::vector<std::vector<double> > > table_t;
         // If we need multiple regions, this class and the SinglePvt* classes must change.
@@ -78,7 +78,7 @@ namespace Opm
         // Oil PVT
         if (phase_usage_.phase_used[Liquid]) {
             if (deck.hasField("PVDO")) {
-                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_));
+                props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(deck.getPVDO().pvdo_, samples));
             } else if (deck.hasField("PVTO")) {
                 props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(deck.getPVTO().pvto_));
             } else if (deck.hasField("PVCDO")) {
@@ -90,7 +90,7 @@ namespace Opm
 	// Gas PVT
         if (phase_usage_.phase_used[Vapour]) {
             if (deck.hasField("PVDG")) {
-                props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_));
+                props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(deck.getPVDG().pvdg_, samples));
             } else if (deck.hasField("PVTG")) {
                 props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(deck.getPVTG().pvtg_));
             } else {

opm/core/fluid/blackoil/BlackoilPvtProperties.hpp

@@ -47,7 +47,7 @@ namespace Opm
         BlackoilPvtProperties();
 
         /// Initialize from deck.
-	void init(const EclipseGridParser& deck);
+    void init(const EclipseGridParser& deck,const int samples = 16);
 
         /// Number of active phases.
         int numPhases() const;

opm/core/fluid/blackoil/SinglePvtConstCompr.hpp

@@ -106,13 +106,18 @@ namespace Opm
                           double* output_B,
                           double* output_dBdp) const
         {
-            B(n, p, 0, output_B);
+            //B(n, p, 0, output_B);
             if (comp_) {
 // #pragma omp parallel for
                 for (int i = 0; i < n; ++i) {
-                    output_dBdp[i] = -comp_*output_B[i];
+                    //output_dBdp[i] = -comp_*output_B[i];
+                    double x = comp_*(p[i] - ref_press_);
+                    double d= (1.0 + x + 0.5*x*x);
+                    output_B[i] = ref_B_/d;//(1.0 + x + 0.5*x*x);
+                    output_dBdp[i] = (- ref_B_/(d*d))*(1 + x) *comp_;
                 }
             } else {
+                std::fill(output_B, output_B + n, ref_B_);
                 std::fill(output_dBdp, output_dBdp + n, 0.0);
             }
         }

opm/core/fluid/blackoil/SinglePvtDead.cpp

@@ -32,9 +32,8 @@ namespace Opm
     //------------------------------------------------------------------------
     // Member functions
     //-------------------------------------------------------------------------
-
     /// Constructor
-    SinglePvtDead::SinglePvtDead(const table_t& pvd_table)
+    SinglePvtDead::SinglePvtDead(const table_t& pvd_table,const int samples)
     {
 	const int region_number = 0;
 	if (pvd_table.size() != 1) {
@@ -51,7 +50,6 @@ namespace Opm
             B_inv[i] = 1.0 / pvd_table[region_number][1][i];
             visc[i]  = pvd_table[region_number][2][i];
 	}
-        int samples = 1025;
         buildUniformMonotoneTable(press, B_inv, samples, one_over_B_);
         buildUniformMonotoneTable(press, visc, samples, viscosity_);
 

opm/core/fluid/blackoil/SinglePvtDead.hpp

@@ -37,8 +37,7 @@ namespace Opm
     {
     public:
 	typedef std::vector<std::vector<std::vector<double> > > table_t;
-
+    SinglePvtDead(const table_t& pvd_table,const int samples = 16);
-	SinglePvtDead(const table_t& pvd_table);
 	virtual ~SinglePvtDead();
 
         /// Viscosity as a function of p and z.

opm/core/utility/UniformTableLinear.hpp

@@ -188,9 +188,14 @@ namespace Opm {
 	UniformTableLinear<T>
 	::derivative(const double xparam) const
 	{
-            // Implements ClosestValue policy.
+        // Implements derivative consistent
+        // with ClosestValue policy for function
+        double value;
+        if(xparam > xmax_ || xparam < xmin_){
+            value=0.0;
+        }else{
             double x = std::min(xparam, xmax_);
-            x = std::max(x, xmin_);
+            x =  std::max(x, xmin_);
 
             // Lookup is easy since we are uniform in x.
             double pos = (x - xmin_)/xdelta_;
@@ -200,7 +205,9 @@ namespace Opm {
                 // We are at xmax_
                 --left;
             }
-	    return (y_values_[left + 1] - y_values_[left])/xdelta_;
+            value = (y_values_[left + 1] - y_values_[left])/xdelta_;
+        }
+        return value;
 	}
 
 

opm/core/utility/initState_impl.hpp

@@ -512,11 +512,12 @@ namespace Opm
                            State& state)
     {
         const int num_phases = props.numPhases();
-        if (num_phases != 2) {
+        
-            THROW("initStateFromDeck(): currently handling only two-phase scenarios.");
-        }
         state.init(grid, num_phases);
         if (deck.hasField("EQUIL")) {
+            if (num_phases != 2) {
+                THROW("initStateFromDeck(): currently handling only two-phase scenarios.");
+            }
             // Set saturations depending on oil-water contact.
             const EQUIL& equil= deck.getEQUIL();
             if (equil.equil.size() != 1) {
@@ -535,11 +536,22 @@ namespace Opm
             const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
             const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
             const int num_cells = grid.number_of_cells;
-            for (int c = 0; c < num_cells; ++c) {
+            if(num_phases == 2){
-                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                for (int c = 0; c < num_cells; ++c) {
-                s[2*c] = sw_deck[c_deck];
+                    int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
-                s[2*c + 1] = 1.0 - s[2*c];
+                    s[2*c] = sw_deck[c_deck];
-                p[c] = p_deck[c_deck];
+                    s[2*c + 1] = 1.0 - s[2*c];
+                    p[c] = p_deck[c_deck];
+                }
+            }else{
+                const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
+                 for (int c = 0; c < num_cells; ++c) {
+                    int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                    s[2*c] = sw_deck[c_deck];
+                    s[2*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
+                    s[2*c + 2] = sg_deck[c_deck];
+                    p[c] = p_deck[c_deck];
+                 }
             }
         } else {
             THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE.");

opm/core/utility/miscUtilities.cpp

@@ -593,9 +593,10 @@ namespace Opm
     {
         int nw = well_bhp.size();
         ASSERT(nw == wells.number_of_wells);
-        if (props.numPhases() != 2) {
+        int np = props.numPhases();
-            THROW("WellReport for now assumes two phase flow.");
+        //if (props.numPhases() != 2) {
-        }
+        //    THROW("WellReport for now assumes two phase flow.");
+        //}
         const double* visc = props.viscosity();
         std::vector<double> data_now;
         data_now.reserve(1 + 3*nw);
@@ -605,7 +606,8 @@ namespace Opm
             double well_rate_total = 0.0;
             double well_rate_water = 0.0;
             for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) {
-                const double perf_rate = well_perfrates[perf]*(unit::day/unit::second);
+              const double perf_rate = unit::convert::to(well_perfrates[perf],
+                                                          unit::cubic(unit::meter)/unit::day);
                 well_rate_total += perf_rate;
                 if (perf_rate > 0.0) {
                     // Injection.
@@ -613,11 +615,14 @@ namespace Opm
                 } else {
                     // Production.
                     const int cell = wells.well_cells[perf];
-                    double mob[2];
+                    double mob[np];
                     props.relperm(1, &saturation[2*cell], &cell, mob, 0);
-                    mob[0] /= visc[0];
+                    double tmob=0;
-                    mob[1] /= visc[1];
+                    for(int i=0; i < np; ++i){
-                    const double fracflow = mob[0]/(mob[0] + mob[1]);
+                        mob[i] /= visc[i];
+                        tmob += mob[i];
+                    }                    
+                    const double fracflow = mob[0]/tmob;
                     well_rate_water += perf_rate*fracflow;
                 }
             }
@@ -646,9 +651,10 @@ namespace Opm
         // TODO: refactor, since this is almost identical to the other push().
         int nw = well_bhp.size();
         ASSERT(nw == wells.number_of_wells);
-        if (props.numPhases() != 2) {
+        int np = props.numPhases();
-            THROW("WellReport for now assumes two phase flow.");
+        //if (props.numPhases() != 2) {
-        }
+        //    THROW("WellReport for now assumes two phase flow.");
+        //}
         std::vector<double> data_now;
         data_now.reserve(1 + 3*nw);
         data_now.push_back(time/unit::day);
@@ -659,20 +665,25 @@ namespace Opm
             for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) {
                 const double perf_rate = well_perfrates[perf]*(unit::day/unit::second);
                 well_rate_total += perf_rate;
-                if (perf_rate > 0.0) {
+                if (perf_rate > 0.0) { 
                     // Injection.
                     well_rate_water += perf_rate*wells.comp_frac[0];
                 } else {
                     // Production.
                     const int cell = wells.well_cells[perf];
-                    double mob[2];
+                    double mob[np];
                     props.relperm(1, &s[2*cell], &cell, mob, 0);
-                    double visc[2];
+                    double visc[np];
                     props.viscosity(1, &p[cell], &z[2*cell], &cell, visc, 0);
-                    mob[0] /= visc[0];
+                    double tmob=0;
-                    mob[1] /= visc[1];
+                    for(int i=0; i < np; ++i){
-                    const double fracflow = mob[0]/(mob[0] + mob[1]);
+                        mob[i] /= visc[i];
+                        tmob += mob[i];
+                    }                    
+                    const double fracflow = mob[0]/(tmob);
                     well_rate_water += perf_rate*fracflow;
+                    //const double fracflow = mob[0]/(tmob);
+                    //well_rate_water += perf_rate*fracflow;
                 }
             }
             data_now.push_back(well_rate_total);

tests/Makefile.am

@@ -21,7 +21,9 @@ test_read_vag					\
 test_readpolymer				\
 test_sf2p					\
 test_writeVtkData				\
-unit_test
+unit_test                                       \
+pvt_test					\
+relperm_test                 
 
 
 bo_resprop_test_SOURCES = bo_resprop_test.cpp
@@ -54,6 +56,12 @@ test_writeVtkData_SOURCES = test_writeVtkData.cpp
 
 unit_test_SOURCES         = unit_test.cpp
 
+pvt_test_SOURCES          = pvt_test.cpp
+pvt_test_LDADD            = $(LDADD)
+
+relperm_test_SOURCES      = relperm_test.cpp
+relperm_test_LDADD       = $(LDADD)
+
 
 if UMFPACK
 noinst_PROGRAMS    += test_cfs_tpfa

tests/pvt_test.cpp

@@ -0,0 +1,138 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/eclipse/EclipseGridParser.hpp>
+#include <opm/core/eclipse/EclipseGridInspector.hpp>
+#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/grid.h>
+
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <iterator>
+#include <string>
+#include <vector>
+
+int main(int argc, char** argv)
+{
+    using namespace std;
+    // Parameters.
+    Opm::parameter::ParameterGroup param(argc, argv);
+
+    // Parser.
+    std::string ecl_file = param.get<std::string>("deck_filename");
+    std::string input_file = param.get<std::string>("input_filename");
+    std::string matrix_output = param.get<std::string>("matrix_output");
+    std::string rs_output = param.get<std::string>("rs_output");
+    std::string b_output = param.get<std::string>("b_output");
+    std::string mu_output = param.get<std::string>("mu_output");
+    Opm::EclipseGridParser deck(ecl_file);
+    UnstructuredGrid grid;
+    grid.number_of_cells = 1;
+    grid.global_cell = NULL;
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, param);
+    Opm::BlackoilPvtProperties pvt;
+    int samples = param.getDefault("dead_tab_size", 1025);
+    pvt.init(deck, samples);
+    
+  
+    const int n = 1;
+    std::fstream inos(input_file.c_str());
+    if(!inos.good()){
+        std::cout << "Could not open :" << input_file << std::endl;
+        exit(3);
+    }
+    std::fstream aos(matrix_output.c_str(), std::fstream::out | std::fstream::trunc);
+    aos << setiosflags(ios::scientific) << setprecision(12);
+    if(!(aos.good())){
+        std::cout << "Could not open"<< matrix_output << std::endl;
+        exit(3);
+    }
+    std::fstream bos(b_output.c_str(), std::fstream::out | std::fstream::trunc);
+    bos << setiosflags(ios::scientific) << setprecision(12);
+    if(!(bos.good())){
+        std::cout << "Could not open"<< b_output << std::endl;
+        exit(3);
+    }
+    std::fstream rsos(rs_output.c_str(), std::fstream::out | std::fstream::trunc);
+    rsos << setiosflags(ios::scientific) << setprecision(12);
+    if(!(rsos.good())){
+        std::cout << "Could not open"<< rs_output << std::endl;
+        exit(3);
+    }
+    std::fstream muos(mu_output.c_str(), std::fstream::out | std::fstream::trunc);
+    if(!(muos.good())){
+        std::cout << "Could not open"<< rs_output << std::endl;
+        exit(3);
+    }
+
+
+
+    int np=props.numPhases();
+    while((inos.good()) && (!inos.eof())){
+        double p[n];
+        double z[np*n];
+        int cells[n] = { 0 };
+        inos >> p[0];      
+        for(int i=0; i < np; ++i){
+            inos >> z[i];
+        }
+
+        if(inos.good()){
+            double A[np*np*n];
+            double dA[np*np*n];
+            props.matrix(n, p, z, cells, A, dA);
+            std::copy(A, A + np*np, std::ostream_iterator<double>(aos, " "));
+            std::copy(dA, dA + np*np, std::ostream_iterator<double>(aos, " "));
+            aos << std::endl;
+            double mu[np];
+            //double dmu[np];//not implemented
+            props.viscosity(n, p, z, cells, mu, 0);
+            std::copy(mu, mu + np, std::ostream_iterator<double>(muos, " "));
+            //std::copy(dmu, dmu + np, std::ostream_iterator<double>(muos, " "));
+            aos << std::endl;
+
+            double b[np];
+            double dbdp[np];
+            pvt.dBdp(n, p, z, b, dbdp);
+            std::copy(b, b + np, std::ostream_iterator<double>(bos, " "));
+            std::copy(dbdp, dbdp + np, std::ostream_iterator<double>(bos, " "));
+            bos << std::endl;
+
+            double rs[np];
+            double drs[np];
+            //pvt.R(n, p, z, rs);
+            pvt.dRdp(n, p, z, rs,drs);
+            std::copy(rs, rs + np, std::ostream_iterator<double>(rsos, " "));
+            std::copy(drs, drs + np, std::ostream_iterator<double>(rsos, " "));
+            rsos << std::endl;
+        }
+    }
+    if (param.anyUnused()) {
+        std::cout << "--------------------   Unused parameters:   --------------------\n";
+        param.displayUsage();
+        std::cout << "----------------------------------------------------------------" << std::endl;
+    }
+}

tests/relperm_test.cpp

@@ -0,0 +1,89 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#if HAVE_CONFIG_H
+#include "config.h"
+#endif // HAVE_CONFIG_H
+
+#include <opm/core/utility/parameters/ParameterGroup.hpp>
+#include <opm/core/eclipse/EclipseGridParser.hpp>
+#include <opm/core/eclipse/EclipseGridInspector.hpp>
+#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
+#include <opm/core/grid.h>
+
+#include <algorithm>
+#include <iostream>
+#include <iterator>
+#include <string>
+#include <vector>
+#include <fstream>
+
+int main(int argc, char** argv)
+{
+    // Parameters.
+    Opm::parameter::ParameterGroup param(argc, argv);
+
+    // Parser.
+    std::string ecl_file = param.get<std::string>("deck_filename");
+    std::string input_file = param.get<std::string>("input_filename");
+    std::string relperm_output = param.get<std::string>("relperm_output");
+    //std::string relperm_output = param.get<std::string>("relperm_outout");
+    Opm::EclipseGridParser deck(ecl_file);
+    UnstructuredGrid grid;
+    grid.number_of_cells = 1;
+    grid.global_cell = NULL;
+    Opm::BlackoilPropertiesFromDeck props(deck, grid, param);
+  
+    std::fstream inos(input_file.c_str());//, std::fstream::in);
+    if(!inos.good()){
+        std::cout << "Could not open :" << input_file << std::endl;
+        exit(3);
+    }
+    std::fstream kros(relperm_output.c_str(), std::fstream::out | std::fstream::trunc);
+    if(!kros.good()){
+        std::cout << "Could not open :" << input_file << std::endl;
+        exit(3);
+    }
+    int np=props.numPhases();
+    while((inos.good()) && (!inos.eof())){
+      double s[np];
+      for(int i=0; i < np; ++i){
+        inos >> s[i];
+      }
+      if(inos.good()){
+          double kr[np];
+          double dkr[np*np];
+          int cell[1];
+          cell[0]=1;
+          props.relperm(1,s, cell, kr, dkr);
+          std::copy(s, s + np, std::ostream_iterator<double>(kros, " "));
+          kros << " ";
+          std::copy(kr, kr + np, std::ostream_iterator<double>(kros, " "));
+          kros << " ";
+          std::copy(dkr, dkr + np*np, std::ostream_iterator<double>(kros, " "));
+          kros << "\n";
+      }
+    }
+    if (param.anyUnused()) {
+        std::cout << "--------------------   Unused parameters:   --------------------\n";
+        param.displayUsage();
+        std::cout << "----------------------------------------------------------------" << std::endl;
+    }
+}