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 9e87bce7..a477f733 100644 --- a/Makefile.am +++ b/Makefile.am @@ -53,6 +53,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 \ @@ -150,6 +152,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..5c834fe8 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); + const 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..13f5799d 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.hpp @@ -26,6 +26,7 @@ #include #include #include +#include struct UnstructuredGrid; @@ -38,12 +39,24 @@ namespace Opm { public: /// Initialize from deck and grid. - /// \param deck Deck input parser - /// \param grid Grid to which property object applies, needed for the + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the /// 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); + + /// Initialize from deck, grid and parameters. + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + /// \param[in] param Parameters. Accepted parameters include: + /// dead_tab_size (1025) number of uniform sample points for dead-oil pvt tables. + /// tab_size_kr (200) number of uniform sample points for saturation tables. + 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..740af43c --- /dev/null +++ b/opm/core/fluid/SatFuncSimple.cpp @@ -0,0 +1,214 @@ +#include +#include +#include +#include +#include +#include +#include +#include +namespace Opm +{ + + + + + 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) { + // A simplified 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) { + // A simplified 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..a3388848 --- /dev/null +++ b/opm/core/fluid/SatFuncSimple.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 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..a03bfb66 --- /dev/null +++ b/opm/core/fluid/SatFuncStone2.cpp @@ -0,0 +1,209 @@ +#include +#include +#include +#include +#include +#include +#include +#include +namespace Opm +{ + + + + + 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..3481b874 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.cpp +++ b/opm/core/fluid/SaturationPropsFromDeck.cpp @@ -90,6 +90,65 @@ namespace Opm } + /// Initialize from deck, grid and parameters + 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. + const 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); + } + } + /// \return P, the number of phases. @@ -192,206 +251,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..39b4ba05 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..191c6fa8 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..557bc1e8 100644 --- a/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp +++ b/opm/core/fluid/blackoil/SinglePvtConstCompr.hpp @@ -106,13 +106,16 @@ namespace Opm double* output_B, double* output_dBdp) const { - 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]; + double x = comp_*(p[i] - ref_press_); + double d = (1.0 + x + 0.5*x*x); + output_B[i] = ref_B_/d; + 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..1a0993be 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..90a7d695 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..2a0fc55f 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -188,19 +188,25 @@ namespace Opm { UniformTableLinear ::derivative(const double xparam) const { - // Implements ClosestValue policy. - double x = std::min(xparam, xmax_); - x = std::max(x, xmin_); - - // Lookup is easy since we are uniform in x. - double pos = (x - xmin_)/xdelta_; - double posi = std::floor(pos); - int left = int(posi); - if (left == int(y_values_.size()) - 1) { - // We are at xmax_ - --left; + // 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_); + // Lookup is easy since we are uniform in x. + double pos = (x - xmin_)/xdelta_; + double posi = std::floor(pos); + int left = int(posi); + if (left == int(y_values_.size()) - 1) { + // We are at xmax_ + --left; + } + value = (y_values_[left + 1] - y_values_[left])/xdelta_; } - return (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..7d01ade0 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -512,11 +512,11 @@ 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(): EQUIL-based init 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 +535,27 @@ 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 if (num_phases == 3) { + if (!deck.hasField("SGAS")) { + THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found)."); + } + 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[3*c] = sw_deck[c_deck]; + s[3*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); + s[3*c + 2] = sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } else { + THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases."); } } 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..aeafbfca 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -40,14 +40,14 @@ namespace Opm /// @param[out] porevol the pore volume by cell. void computePorevolume(const UnstructuredGrid& grid, const double* porosity, - std::vector& porevol) + std::vector& porevol) { - int num_cells = grid.number_of_cells; - porevol.resize(num_cells); - std::transform(porosity, porosity + num_cells, - grid.cell_volumes, - porevol.begin(), - std::multiplies()); + int num_cells = grid.number_of_cells; + porevol.resize(num_cells); + std::transform(porosity, porosity + num_cells, + grid.cell_volumes, + porevol.begin(), + std::multiplies()); } @@ -98,20 +98,20 @@ namespace Opm /// For each phase p, we compute /// sat_vol_p = sum_i s_p_i pv_i void computeSaturatedVol(const std::vector& pv, - const std::vector& s, - double* sat_vol) + const std::vector& s, + double* sat_vol) { - const int num_cells = pv.size(); - const int np = s.size()/pv.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and pv vectors do not match."); - } - std::fill(sat_vol, sat_vol + np, 0.0); - for (int c = 0; c < num_cells; ++c) { - for (int p = 0; p < np; ++p) { - sat_vol[p] += pv[c]*s[np*c + p]; - } - } + const int num_cells = pv.size(); + const int np = s.size()/pv.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and pv vectors do not match."); + } + std::fill(sat_vol, sat_vol + np, 0.0); + for (int c = 0; c < num_cells; ++c) { + for (int p = 0; p < np; ++p) { + sat_vol[p] += pv[c]*s[np*c + p]; + } + } } @@ -123,28 +123,28 @@ namespace Opm /// For each phase p, we compute /// aver_sat_p = (sum_i s_p_i pv_i) / (sum_i pv_i). void computeAverageSat(const std::vector& pv, - const std::vector& s, - double* aver_sat) + const std::vector& s, + double* aver_sat) { - const int num_cells = pv.size(); - const int np = s.size()/pv.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and pv vectors do not match."); - } - double tot_pv = 0.0; - // Note that we abuse the output array to accumulate the - // saturated pore volumes. - std::fill(aver_sat, aver_sat + np, 0.0); - for (int c = 0; c < num_cells; ++c) { - tot_pv += pv[c]; - for (int p = 0; p < np; ++p) { - aver_sat[p] += pv[c]*s[np*c + p]; - } - } - // Must divide by pore volumes to get saturations. - for (int p = 0; p < np; ++p) { - aver_sat[p] /= tot_pv; - } + const int num_cells = pv.size(); + const int np = s.size()/pv.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and pv vectors do not match."); + } + double tot_pv = 0.0; + // Note that we abuse the output array to accumulate the + // saturated pore volumes. + std::fill(aver_sat, aver_sat + np, 0.0); + for (int c = 0; c < num_cells; ++c) { + tot_pv += pv[c]; + for (int p = 0; p < np; ++p) { + aver_sat[p] += pv[c]*s[np*c + p]; + } + } + // Must divide by pore volumes to get saturations. + for (int p = 0; p < np; ++p) { + aver_sat[p] /= tot_pv; + } } @@ -161,38 +161,38 @@ namespace Opm /// where P = s.size()/src.size(). /// @param[out] produced must also point to a valid array with P elements. void computeInjectedProduced(const IncompPropertiesInterface& props, - const std::vector& s, - const std::vector& src, - const double dt, - double* injected, - double* produced) + const std::vector& s, + const std::vector& src, + const double dt, + double* injected, + double* produced) { - const int num_cells = src.size(); - const int np = s.size()/src.size(); - if (int(s.size()) != num_cells*np) { - THROW("Sizes of s and src vectors do not match."); - } - std::fill(injected, injected + np, 0.0); - std::fill(produced, produced + np, 0.0); - const double* visc = props.viscosity(); - std::vector mob(np); - for (int c = 0; c < num_cells; ++c) { - if (src[c] > 0.0) { - injected[0] += src[c]*dt; - } else if (src[c] < 0.0) { - const double flux = -src[c]*dt; - const double* sat = &s[np*c]; - props.relperm(1, sat, &c, &mob[0], 0); - double totmob = 0.0; - for (int p = 0; p < np; ++p) { - mob[p] /= visc[p]; - totmob += mob[p]; - } - for (int p = 0; p < np; ++p) { - produced[p] += (mob[p]/totmob)*flux; - } - } - } + const int num_cells = src.size(); + const int np = s.size()/src.size(); + if (int(s.size()) != num_cells*np) { + THROW("Sizes of s and src vectors do not match."); + } + std::fill(injected, injected + np, 0.0); + std::fill(produced, produced + np, 0.0); + const double* visc = props.viscosity(); + std::vector mob(np); + for (int c = 0; c < num_cells; ++c) { + if (src[c] > 0.0) { + injected[0] += src[c]*dt; + } else if (src[c] < 0.0) { + const double flux = -src[c]*dt; + const double* sat = &s[np*c]; + props.relperm(1, sat, &c, &mob[0], 0); + double totmob = 0.0; + for (int p = 0; p < np; ++p) { + mob[p] /= visc[p]; + totmob += mob[p]; + } + for (int p = 0; p < np; ++p) { + produced[p] += (mob[p]/totmob)*flux; + } + } + } } @@ -203,9 +203,9 @@ namespace Opm /// @param[in] s saturation values (for all phases) /// @param[out] totmob total mobilities. void computeTotalMobility(const Opm::IncompPropertiesInterface& props, - const std::vector& cells, - const std::vector& s, - std::vector& totmob) + const std::vector& cells, + const std::vector& s, + std::vector& totmob) { std::vector pmobc; @@ -231,10 +231,10 @@ namespace Opm /// @param[out] totmob total mobility /// @param[out] omega fractional-flow weighted fluid densities. void computeTotalMobilityOmega(const Opm::IncompPropertiesInterface& props, - const std::vector& cells, - const std::vector& s, - std::vector& totmob, - std::vector& omega) + const std::vector& cells, + const std::vector& s, + std::vector& totmob, + std::vector& omega) { std::vector pmobc; @@ -331,33 +331,33 @@ namespace Opm /// (+) positive inflow of first phase (water) /// (-) negative total outflow of both phases void computeTransportSource(const UnstructuredGrid& grid, - const std::vector& src, - const std::vector& faceflux, - const double inflow_frac, + const std::vector& src, + const std::vector& faceflux, + const double inflow_frac, const Wells* wells, const std::vector& well_perfrates, - std::vector& transport_src) + std::vector& transport_src) { - int nc = grid.number_of_cells; - transport_src.resize(nc); + int nc = grid.number_of_cells; + transport_src.resize(nc); // Source term and boundary contributions. - for (int c = 0; c < nc; ++c) { - transport_src[c] = 0.0; - transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; - for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) { - int f = grid.cell_faces[hf]; - const int* f2c = &grid.face_cells[2*f]; - double bdy_influx = 0.0; - if (f2c[0] == c && f2c[1] == -1) { - bdy_influx = -faceflux[f]; - } else if (f2c[0] == -1 && f2c[1] == c) { - bdy_influx = faceflux[f]; - } - if (bdy_influx != 0.0) { - transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; - } - } - } + for (int c = 0; c < nc; ++c) { + transport_src[c] = 0.0; + transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; + for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* f2c = &grid.face_cells[2*f]; + double bdy_influx = 0.0; + if (f2c[0] == c && f2c[1] == -1) { + bdy_influx = -faceflux[f]; + } else if (f2c[0] == -1 && f2c[1] == c) { + bdy_influx = faceflux[f]; + } + if (bdy_influx != 0.0) { + transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; + } + } + } // Well contributions. if (wells) { @@ -392,52 +392,52 @@ namespace Opm /// @param[in] face_flux signed per-face fluxes /// @param[out] cell_velocity the estimated velocities. void estimateCellVelocity(const UnstructuredGrid& grid, - const std::vector& face_flux, - std::vector& cell_velocity) + const std::vector& face_flux, + std::vector& cell_velocity) { - const int dim = grid.dimensions; - cell_velocity.clear(); - cell_velocity.resize(grid.number_of_cells*dim, 0.0); - for (int face = 0; face < grid.number_of_faces; ++face) { - int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; - const double* fc = &grid.face_centroids[face*dim]; - double flux = face_flux[face]; - for (int i = 0; i < 2; ++i) { - if (c[i] >= 0) { - const double* cc = &grid.cell_centroids[c[i]*dim]; - for (int d = 0; d < dim; ++d) { - double v_contrib = fc[d] - cc[d]; - v_contrib *= flux/grid.cell_volumes[c[i]]; - cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; - } - } - } - } + const int dim = grid.dimensions; + cell_velocity.clear(); + cell_velocity.resize(grid.number_of_cells*dim, 0.0); + for (int face = 0; face < grid.number_of_faces; ++face) { + int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; + const double* fc = &grid.face_centroids[face*dim]; + double flux = face_flux[face]; + for (int i = 0; i < 2; ++i) { + if (c[i] >= 0) { + const double* cc = &grid.cell_centroids[c[i]*dim]; + for (int d = 0; d < dim; ++d) { + double v_contrib = fc[d] - cc[d]; + v_contrib *= flux/grid.cell_volumes[c[i]]; + cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; + } + } + } + } } /// Extract a vector of water saturations from a vector of /// interleaved water and oil saturations. void toWaterSat(const std::vector& sboth, - std::vector& sw) + std::vector& sw) { - int num = sboth.size()/2; - sw.resize(num); - for (int i = 0; i < num; ++i) { - sw[i] = sboth[2*i]; - } + int num = sboth.size()/2; + sw.resize(num); + for (int i = 0; i < num; ++i) { + sw[i] = sboth[2*i]; + } } /// Make a a vector of interleaved water and oil saturations from /// a vector of water saturations. void toBothSat(const std::vector& sw, - std::vector& sboth) + std::vector& sboth) { - int num = sw.size(); - sboth.resize(2*num); - for (int i = 0; i < num; ++i) { - sboth[2*i] = sw[i]; - sboth[2*i + 1] = 1.0 - sw[i]; - } + int num = sw.size(); + sboth.resize(2*num); + for (int i = 0; i < num; ++i) { + sboth[2*i] = sw[i]; + sboth[2*i + 1] = 1.0 - sw[i]; + } } @@ -450,30 +450,30 @@ namespace Opm if (np != 2) { THROW("wellsToSrc() requires a 2 phase case."); } - src.resize(num_cells); - for (int w = 0; w < wells.number_of_wells; ++w) { + src.resize(num_cells); + for (int w = 0; w < wells.number_of_wells; ++w) { const int cur = wells.ctrls[w]->current; - if (wells.ctrls[w]->num != 1) { - MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored."); - } - if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) { - THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE."); - } - if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) { - THROW("In wellsToSrc(): well has multiple perforations."); - } + if (wells.ctrls[w]->num != 1) { + MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored."); + } + if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) { + THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE."); + } + if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) { + THROW("In wellsToSrc(): well has multiple perforations."); + } for (int p = 0; p < np; ++p) { if (wells.ctrls[w]->distr[np*cur + p] != 1.0) { THROW("In wellsToSrc(): well not controlled on total rate."); } } - double flow = wells.ctrls[w]->target[cur]; + double flow = wells.ctrls[w]->target[cur]; if (wells.type[w] == INJECTOR) { flow *= wells.comp_frac[np*w + 0]; // Obtaining water rate for inflow source. } - const double cell = wells.well_cells[wells.well_connpos[w]]; - src[cell] = flow; - } + const double cell = wells.well_cells[wells.well_connpos[w]]; + src[cell] = flow; + } } @@ -593,8 +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(); + const int max_np = 3; + if (np > max_np) { + THROW("WellReport for now assumes #phases <= " << max_np); } const double* visc = props.viscosity(); std::vector data_now; @@ -605,7 +607,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 +616,14 @@ namespace Opm } else { // Production. const int cell = wells.well_cells[perf]; - double mob[2]; + double mob[max_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,8 +652,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(); + const int max_np = 3; + if (np > max_np) { + THROW("WellReport for now assumes #phases <= " << max_np); } std::vector data_now; data_now.reserve(1 + 3*nw); @@ -657,7 +665,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. @@ -665,13 +674,16 @@ namespace Opm } else { // Production. const int cell = wells.well_cells[perf]; - double mob[2]; - props.relperm(1, &s[2*cell], &cell, mob, 0); - double visc[2]; - 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 mob[max_np]; + props.relperm(1, &s[np*cell], &cell, mob, 0); + double visc[max_np]; + props.viscosity(1, &p[cell], &z[np*cell], &cell, visc, 0); + 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; } } diff --git a/tests/Makefile.am b/tests/Makefile.am index f23270e7..b1e70503 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -12,6 +12,8 @@ noinst_PROGRAMS = \ bo_resprop_test \ monotcubicinterpolator_test \ param_test \ +pvt_test \ +relperm_test \ sparsetable_test \ sparsevector_test \ test_cartgrid \ @@ -32,6 +34,10 @@ monotcubicinterpolator_test_SOURCES = monotcubicinterpolator_test.cpp param_test_SOURCES = param_test.cpp param_test_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) +pvt_test_SOURCES = pvt_test.cpp + +relperm_test_SOURCES = relperm_test.cpp + sparsetable_test_SOURCES = sparsetable_test.cpp sparsetable_test_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) diff --git a/tests/pvt_test.cpp b/tests/pvt_test.cpp new file mode 100644 index 00000000..a175885a --- /dev/null +++ b/tests/pvt_test.cpp @@ -0,0 +1,146 @@ +/* + 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; + grid.dimensions = 3; + grid.cartdims[0] = 1; + grid.cartdims[1] = 1; + grid.cartdims[2] = 1; + 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); + } + + + + const int np = props.numPhases(); + const int max_np = 3; + if (np > max_np) { + THROW("Max #phases is 3."); + } + while((inos.good()) && (!inos.eof())){ + double p[n]; + double z[max_np*n]; + int cells[n] = { 0 }; + inos >> p[0]; + for(int i=0; i < np; ++i){ + inos >> z[i]; + } + + if(inos.good()){ + double A[max_np*max_np*n]; + double dA[max_np*max_np*n]; + props.matrix(n, p, z, cells, A, dA); + std::copy(A, A + np*np*n, std::ostream_iterator(aos, " ")); + std::copy(dA, dA + np*np*n, std::ostream_iterator(aos, " ")); + aos << std::endl; + double mu[max_np]; + //double dmu[max_np];//not implemented + props.viscosity(n, p, z, cells, mu, 0); + std::copy(mu, mu + np*n, std::ostream_iterator(muos, " ")); + //std::copy(dmu, dmu + np*n, std::ostream_iterator(muos, " ")); + aos << std::endl; + + double b[max_np]; + double dbdp[max_np]; + pvt.dBdp(n, p, z, b, dbdp); + std::copy(b, b + np*n, std::ostream_iterator(bos, " ")); + std::copy(dbdp, dbdp + np*n, std::ostream_iterator(bos, " ")); + bos << std::endl; + + double rs[max_np]; + double drs[max_np]; + //pvt.R(n, p, z, rs); + pvt.dRdp(n, p, z, rs,drs); + std::copy(rs, rs + np*n, std::ostream_iterator(rsos, " ")); + std::copy(drs, drs + np*n, 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..4c777071 --- /dev/null +++ b/tests/relperm_test.cpp @@ -0,0 +1,97 @@ +/* + 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; + grid.dimensions = 3; + grid.cartdims[0] = 1; + grid.cartdims[1] = 1; + grid.cartdims[2] = 1; + 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); + } + const int np = props.numPhases(); + const int max_np = 3; + if (np > max_np) { + THROW("Max #phases is 3."); + } + while((inos.good()) && (!inos.eof())){ + double s[max_np]; + for(int i=0; i < np; ++i){ + inos >> s[i]; + } + if(inos.good()){ + double kr[max_np]; + double dkr[max_np*max_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; + } +}