diff --git a/CMakeLists.txt b/CMakeLists.txt index a6e3b690..eec84d9c 100644 --- a/CMakeLists.txt +++ b/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 -opm/core/eclipse/EclipseGridInspector.cpp -opm/core/eclipse/EclipseGridParser.cpp -opm/core/fluid/blackoil/BlackoilPvtProperties.cpp -opm/core/fluid/blackoil/SinglePvtDead.cpp -opm/core/fluid/blackoil/SinglePvtLiveGas.cpp -opm/core/fluid/blackoil/SinglePvtLiveOil.cpp -opm/core/fluid/blackoil/SinglePvtInterface.cpp -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 -) +FILE(GLOB_RECURSE C_FILES_CORE "opm/core/*.c") +FILE(GLOB_RECURSE CPP_FILES_CORE "opm/core/*.cpp") +FILE(GLOB_RECURSE REMOVE_FILES "processgrid.c") +FILE(GLOB_RECURSE REMOVE_FILESMX "mx*.c") +FILE(GLOB_RECURSE REMOVE_FILESAGMG "*AGMG.cpp" "*test*") +list(REMOVE_ITEM C_FILES_CORE ${REMOVE_FILES} ${REMOVE_FILESMX} ) +list(REMOVE_ITEM CPP_FILES_CORE ${REMOVE_FILES} ${REMOVE_FILESAGMG}) +add_library(opmcore ${C_FILES_CORE} ${CPP_FILES_CORE} ) + 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) diff --git a/FindUmfPack.cmake b/FindUmfPack.cmake new file mode 100644 index 00000000..5fe5db16 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/Makefile.am b/Makefile.am index 712d7c65..1835e89e 100644 --- a/Makefile.am +++ b/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 \ diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 004e6f88..67943fc7 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -18,7 +18,7 @@ */ #include - +#include 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() { } diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp index a2a6c3de..1d0a0d38 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -26,6 +26,7 @@ #include #include #include +#include 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(); diff --git a/opm/core/fluid/SatFuncSimple.cpp b/opm/core/fluid/SatFuncSimple.cpp new file mode 100644 index 00000000..90ca0e56 --- /dev/null +++ b/opm/core/fluid/SatFuncSimple.cpp @@ -0,0 +1,213 @@ +#include +#include +#include +#include +#include +#include +#include +#include +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& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& 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& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& 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 diff --git a/opm/core/fluid/SatFuncSimple.hpp b/opm/core/fluid/SatFuncSimple.hpp new file mode 100644 index 00000000..93ba7a31 --- /dev/null +++ b/opm/core/fluid/SatFuncSimple.hpp @@ -0,0 +1,32 @@ +#ifndef SATFUNCSIMPLE_HPP +#define SATFUNCSIMPLE_HPP +#include +#include +#include +#include +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 krw_; + UniformTableLinear krow_; + UniformTableLinear pcow_; + UniformTableLinear krg_; + UniformTableLinear krog_; + UniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; +} // namespace Opm +#endif // SATFUNCSIMPLE_HPP diff --git a/opm/core/fluid/SatFuncStone2.cpp b/opm/core/fluid/SatFuncStone2.cpp new file mode 100644 index 00000000..c43238c1 --- /dev/null +++ b/opm/core/fluid/SatFuncStone2.cpp @@ -0,0 +1,209 @@ +#include +#include +#include +#include +#include +#include +#include +#include +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& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& 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& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& 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 diff --git a/opm/core/fluid/SatFuncStone2.hpp b/opm/core/fluid/SatFuncStone2.hpp new file mode 100644 index 00000000..ae4a6d38 --- /dev/null +++ b/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 . +*/ +#ifndef SATFUNCSTONE2_HPP +#define SATFUNCSTONE2_HPP +#include +#include +#include +#include +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 krw_; + UniformTableLinear krow_; + UniformTableLinear pcow_; + UniformTableLinear krg_; + UniformTableLinear krog_; + UniformTableLinear pcog_; + double krocw_; // = krow_(s_wc) + }; +} // namespace Opm +#endif // SATFUNCSTONE2_HPP diff --git a/opm/core/fluid/SaturationPropsFromDeck.cpp b/opm/core/fluid/SaturationPropsFromDeck.cpp index 36cf8364..582c29db 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.cpp +++ b/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& 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& sw = swof_table[table_num][0]; - const std::vector& krw = swof_table[table_num][1]; - const std::vector& krow = swof_table[table_num][2]; - const std::vector& 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& sg = sgof_table[table_num][0]; - const std::vector& krg = sgof_table[table_num][1]; - const std::vector& krog = sgof_table[table_num][2]; - const std::vector& 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 diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index b83a1271..bd70d802 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -19,10 +19,12 @@ #ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED #define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED - +#include #include #include #include +#include +#include #include 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 - { - 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 krw_; - UniformTableLinear krow_; - UniformTableLinear pcow_; - UniformTableLinear krg_; - UniformTableLinear krog_; - UniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) - }; - std::vector satfuncset_; + typedef SatFuncSimple satfunc_t; + std::vector satfuncset_; std::vector cell_to_func_; // = SATNUM - 1 - const SatFuncSet& funcForCell(const int cell) const; + const satfunc_t& funcForCell(const int cell) const; }; diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp index 9fdaf445..6278d1bd 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.cpp +++ b/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 > > 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 { diff --git a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp b/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp index 627f266b..99fd5b8f 100644 --- a/opm/core/fluid/blackoil/BlackoilPvtProperties.hpp +++ b/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; diff --git a/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp b/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp index 1bc5638f..3044bdb2 100644 --- a/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp +++ b/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); } } diff --git a/opm/core/fluid/blackoil/SinglePvtDead.cpp b/opm/core/fluid/blackoil/SinglePvtDead.cpp index 3d1a2740..b37fdc2d 100644 --- a/opm/core/fluid/blackoil/SinglePvtDead.cpp +++ b/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_); diff --git a/opm/core/fluid/blackoil/SinglePvtDead.hpp b/opm/core/fluid/blackoil/SinglePvtDead.hpp index 137e2e30..5a8f6971 100644 --- a/opm/core/fluid/blackoil/SinglePvtDead.hpp +++ b/opm/core/fluid/blackoil/SinglePvtDead.hpp @@ -37,8 +37,7 @@ namespace Opm { public: typedef std::vector > > table_t; - - SinglePvtDead(const table_t& pvd_table); + SinglePvtDead(const table_t& pvd_table,const int samples = 16); virtual ~SinglePvtDead(); /// Viscosity as a function of p and z. diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 11c2d4ca..60b3433d 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -188,9 +188,14 @@ namespace Opm { UniformTableLinear ::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; } diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index c50db8a9..f6670434 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/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& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& p_deck = deck.getFloatingPointValue("PRESSURE"); const int num_cells = grid.number_of_cells; - 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 - s[2*c]; - p[c] = p_deck[c_deck]; + if(num_phases == 2){ + 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 - s[2*c]; + p[c] = p_deck[c_deck]; + } + }else{ + const std::vector& 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."); diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index ac9552dc..0b8d2010 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/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) { - THROW("WellReport for now assumes two phase flow."); - } + int np = props.numPhases(); + //if (props.numPhases() != 2) { + // THROW("WellReport for now assumes two phase flow."); + //} const double* visc = props.viscosity(); std::vector 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]; - mob[1] /= visc[1]; - const double fracflow = mob[0]/(mob[0] + mob[1]); + double tmob=0; + for(int i=0; i < np; ++i){ + 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) { - THROW("WellReport for now assumes two phase flow."); - } + int np = props.numPhases(); + //if (props.numPhases() != 2) { + // THROW("WellReport for now assumes two phase flow."); + //} std::vector 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]; - mob[1] /= visc[1]; - const double fracflow = mob[0]/(mob[0] + mob[1]); + double tmob=0; + for(int i=0; i < np; ++i){ + 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); diff --git a/tests/Makefile.am b/tests/Makefile.am index a82a1dc3..4181eaea 100644 --- a/tests/Makefile.am +++ b/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 diff --git a/tests/pvt_test.cpp b/tests/pvt_test.cpp new file mode 100644 index 00000000..b8cc3528 --- /dev/null +++ b/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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +int main(int argc, char** argv) +{ + using namespace std; + // Parameters. + Opm::parameter::ParameterGroup param(argc, argv); + + // Parser. + std::string ecl_file = param.get("deck_filename"); + std::string input_file = param.get("input_filename"); + std::string matrix_output = param.get("matrix_output"); + std::string rs_output = param.get("rs_output"); + std::string b_output = param.get("b_output"); + std::string mu_output = param.get("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(aos, " ")); + std::copy(dA, dA + np*np, std::ostream_iterator(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(muos, " ")); + //std::copy(dmu, dmu + np, std::ostream_iterator(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(bos, " ")); + std::copy(dbdp, dbdp + np, std::ostream_iterator(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(rsos, " ")); + std::copy(drs, drs + np, std::ostream_iterator(rsos, " ")); + rsos << std::endl; + } + } + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } +} diff --git a/tests/relperm_test.cpp b/tests/relperm_test.cpp new file mode 100644 index 00000000..f521a483 --- /dev/null +++ b/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 . +*/ + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +int main(int argc, char** argv) +{ + // Parameters. + Opm::parameter::ParameterGroup param(argc, argv); + + // Parser. + std::string ecl_file = param.get("deck_filename"); + std::string input_file = param.get("input_filename"); + std::string relperm_output = param.get("relperm_output"); + //std::string relperm_output = param.get("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(kros, " ")); + kros << " "; + std::copy(kr, kr + np, std::ostream_iterator(kros, " ")); + kros << " "; + std::copy(dkr, dkr + np*np, std::ostream_iterator(kros, " ")); + kros << "\n"; + } + } + if (param.anyUnused()) { + std::cout << "-------------------- Unused parameters: --------------------\n"; + param.displayUsage(); + std::cout << "----------------------------------------------------------------" << std::endl; + } +}