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

Merge remote-tracking branch 'upstream/master'

Browse Source
This commit is contained in:
Bård Skaflestad 2012-09-04 13:58:43 +02:00
commit 8b9ddc639c
22 changed files with 1205 additions and 501 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

120
CMakeLists.txt
View File

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

56
FindUmfPack.cmake Normal file
View File

@ -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)

4
Makefile.am
View File

@ -53,6 +53,8 @@ opm/core/fluid/RockCompressibility.cpp \
opm/core/fluid/RockFromDeck.cpp \ opm/core/fluid/RockFromDeck.cpp \
opm/core/fluid/SaturationPropsBasic.cpp \ opm/core/fluid/SaturationPropsBasic.cpp \
opm/core/fluid/SaturationPropsFromDeck.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/BlackoilPvtProperties.cpp \
opm/core/fluid/blackoil/SinglePvtDead.cpp \ opm/core/fluid/blackoil/SinglePvtDead.cpp \
opm/core/fluid/blackoil/SinglePvtInterface.cpp \ opm/core/fluid/blackoil/SinglePvtInterface.cpp \
@ -150,6 +152,8 @@ opm/core/fluid/RockCompressibility.hpp \
opm/core/fluid/RockFromDeck.hpp \ opm/core/fluid/RockFromDeck.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \ opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \ opm/core/fluid/SaturationPropsFromDeck.hpp \
opm/core/fluid/SatFuncStone2.hpp \
opm/core/fluid/SatFuncSimple.hpp \
opm/core/fluid/SimpleFluid2p.hpp \ opm/core/fluid/SimpleFluid2p.hpp \
opm/core/fluid/blackoil/BlackoilPhases.hpp \ opm/core/fluid/blackoil/BlackoilPhases.hpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \ opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \

17
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
View File

@ -18,7 +18,7 @@
*/ */
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp> #include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
namespace Opm 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() BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()
{ {
} }

19
opm/core/fluid/BlackoilPropertiesFromDeck.hpp
View File

@ -26,6 +26,7 @@
#include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp> #include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp> #include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp> #include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
struct UnstructuredGrid; struct UnstructuredGrid;
@ -38,12 +39,24 @@ namespace Opm
{ {
public: public:
/// Initialize from deck and grid. /// Initialize from deck and grid.
/// \param deck Deck input parser /// \param[in] deck Deck input parser
/// \param grid Grid to which property object applies, needed for the /// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid) /// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck. /// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(const EclipseGridParser& 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. /// Destructor.
virtual ~BlackoilPropertiesFromDeck(); virtual ~BlackoilPropertiesFromDeck();

214
opm/core/fluid/SatFuncSimple.cpp Normal file
View File

@ -0,0 +1,214 @@
#include <opm/core/fluid/SatFuncSimple.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>
namespace Opm
{
void SatFuncSimple::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg)
{
init(deck, table_num, phase_usg, 200);
}
void SatFuncSimple::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
buildUniformMonotoneTable(sw, krw, samples, krw_);
buildUniformMonotoneTable(sw, krow, samples, krow_);
buildUniformMonotoneTable(sw, pcow, samples, pcow_);
krocw_ = krow[0]; // At connate water -> ecl. SWOF
swco = sw[0];
smin_[phase_usage.phase_pos[Aqua]] = sw[0];
swmax = sw.back();
smax_[phase_usage.phase_pos[Aqua]] = sw.back();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
buildUniformMonotoneTable(sg, krg, samples, krg_);
buildUniformMonotoneTable(sg, krog, samples, krog_);
buildUniformMonotoneTable(sg, pcog, samples, pcog_);
smin_[phase_usage.phase_pos[Vapour]] = sg[0];
if (std::fabs(sg.back() + swco - 1.0) > 1e-3) {
THROW("Gas maximum saturation in SGOF table = " << sg.back() <<
", should equal (1.0 - connate water sat) = " << (1.0 - swco));
}
smax_[phase_usage.phase_pos[Vapour]] = sg.back();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncSimple::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// 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

50
opm/core/fluid/SatFuncSimple.hpp Normal file
View File

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

209
opm/core/fluid/SatFuncStone2.cpp Normal file
View File

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

50
opm/core/fluid/SatFuncStone2.hpp Normal file
View File

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

255
opm/core/fluid/SaturationPropsFromDeck.cpp
View File

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

33
opm/core/fluid/SaturationPropsFromDeck.hpp
View File

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

6
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp
View File

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

2
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp
View File

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

7
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp
View File

@ -106,13 +106,16 @@ namespace Opm
double* output_B, double* output_B,
double* output_dBdp) const double* output_dBdp) const
{ {
B(n, p, 0, output_B);
if (comp_) { if (comp_) {
// #pragma omp parallel for // #pragma omp parallel for
for (int i = 0; i < n; ++i) { 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 { } else {
std::fill(output_B, output_B + n, ref_B_);
std::fill(output_dBdp, output_dBdp + n, 0.0); std::fill(output_dBdp, output_dBdp + n, 0.0);
} }
} }

4
opm/core/fluid/blackoil/SinglePvtDead.cpp
View File

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

3
opm/core/fluid/blackoil/SinglePvtDead.hpp
View File

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

30
opm/core/utility/UniformTableLinear.hpp
View File

@ -188,19 +188,25 @@ namespace Opm {
UniformTableLinear<T> UniformTableLinear<T>
::derivative(const double xparam) const ::derivative(const double xparam) const
{ {
// Implements ClosestValue policy. // Implements derivative consistent
double x = std::min(xparam, xmax_); // with ClosestValue policy for function
x = std::max(x, xmin_); double value;
if (xparam > xmax_ || xparam < xmin_) {
// Lookup is easy since we are uniform in x. value = 0.0;
double pos = (x - xmin_)/xdelta_; } else {
double posi = std::floor(pos); double x = std::min(xparam, xmax_);
int left = int(posi); x = std::max(x, xmin_);
if (left == int(y_values_.size()) - 1) { // Lookup is easy since we are uniform in x.
// We are at xmax_ double pos = (x - xmin_)/xdelta_;
--left; 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;
} }

32
opm/core/utility/initState_impl.hpp
View File

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

346
opm/core/utility/miscUtilities.cpp
View File

@ -40,14 +40,14 @@ namespace Opm
/// @param[out] porevol the pore volume by cell. /// @param[out] porevol the pore volume by cell.
void computePorevolume(const UnstructuredGrid& grid, void computePorevolume(const UnstructuredGrid& grid,
const double* porosity, const double* porosity,
std::vector<double>& porevol) std::vector<double>& porevol)
{ {
int num_cells = grid.number_of_cells; int num_cells = grid.number_of_cells;
porevol.resize(num_cells); porevol.resize(num_cells);
std::transform(porosity, porosity + num_cells, std::transform(porosity, porosity + num_cells,
grid.cell_volumes, grid.cell_volumes,
porevol.begin(), porevol.begin(),
std::multiplies<double>()); std::multiplies<double>());
} }
@ -98,20 +98,20 @@ namespace Opm
/// For each phase p, we compute /// For each phase p, we compute
/// sat_vol_p = sum_i s_p_i pv_i /// sat_vol_p = sum_i s_p_i pv_i
void computeSaturatedVol(const std::vector<double>& pv, void computeSaturatedVol(const std::vector<double>& pv,
const std::vector<double>& s, const std::vector<double>& s,
double* sat_vol) double* sat_vol)
{ {
const int num_cells = pv.size(); const int num_cells = pv.size();
const int np = s.size()/pv.size(); const int np = s.size()/pv.size();
if (int(s.size()) != num_cells*np) { if (int(s.size()) != num_cells*np) {
THROW("Sizes of s and pv vectors do not match."); THROW("Sizes of s and pv vectors do not match.");
} }
std::fill(sat_vol, sat_vol + np, 0.0); std::fill(sat_vol, sat_vol + np, 0.0);
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
sat_vol[p] += pv[c]*s[np*c + p]; sat_vol[p] += pv[c]*s[np*c + p];
} }
} }
} }
@ -123,28 +123,28 @@ namespace Opm
/// For each phase p, we compute /// For each phase p, we compute
/// aver_sat_p = (sum_i s_p_i pv_i) / (sum_i pv_i). /// aver_sat_p = (sum_i s_p_i pv_i) / (sum_i pv_i).
void computeAverageSat(const std::vector<double>& pv, void computeAverageSat(const std::vector<double>& pv,
const std::vector<double>& s, const std::vector<double>& s,
double* aver_sat) double* aver_sat)
{ {
const int num_cells = pv.size(); const int num_cells = pv.size();
const int np = s.size()/pv.size(); const int np = s.size()/pv.size();
if (int(s.size()) != num_cells*np) { if (int(s.size()) != num_cells*np) {
THROW("Sizes of s and pv vectors do not match."); THROW("Sizes of s and pv vectors do not match.");
} }
double tot_pv = 0.0; double tot_pv = 0.0;
// Note that we abuse the output array to accumulate the // Note that we abuse the output array to accumulate the
// saturated pore volumes. // saturated pore volumes.
std::fill(aver_sat, aver_sat + np, 0.0); std::fill(aver_sat, aver_sat + np, 0.0);
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
tot_pv += pv[c]; tot_pv += pv[c];
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
aver_sat[p] += pv[c]*s[np*c + p]; aver_sat[p] += pv[c]*s[np*c + p];
} }
} }
// Must divide by pore volumes to get saturations. // Must divide by pore volumes to get saturations.
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
aver_sat[p] /= tot_pv; aver_sat[p] /= tot_pv;
} }
} }
@ -161,38 +161,38 @@ namespace Opm
/// where P = s.size()/src.size(). /// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements. /// @param[out] produced must also point to a valid array with P elements.
void computeInjectedProduced(const IncompPropertiesInterface& props, void computeInjectedProduced(const IncompPropertiesInterface& props,
const std::vector<double>& s, const std::vector<double>& s,
const std::vector<double>& src, const std::vector<double>& src,
const double dt, const double dt,
double* injected, double* injected,
double* produced) double* produced)
{ {
const int num_cells = src.size(); const int num_cells = src.size();
const int np = s.size()/src.size(); const int np = s.size()/src.size();
if (int(s.size()) != num_cells*np) { if (int(s.size()) != num_cells*np) {
THROW("Sizes of s and src vectors do not match."); THROW("Sizes of s and src vectors do not match.");
} }
std::fill(injected, injected + np, 0.0); std::fill(injected, injected + np, 0.0);
std::fill(produced, produced + np, 0.0); std::fill(produced, produced + np, 0.0);
const double* visc = props.viscosity(); const double* visc = props.viscosity();
std::vector<double> mob(np); std::vector<double> mob(np);
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
if (src[c] > 0.0) { if (src[c] > 0.0) {
injected[0] += src[c]*dt; injected[0] += src[c]*dt;
} else if (src[c] < 0.0) { } else if (src[c] < 0.0) {
const double flux = -src[c]*dt; const double flux = -src[c]*dt;
const double* sat = &s[np*c]; const double* sat = &s[np*c];
props.relperm(1, sat, &c, &mob[0], 0); props.relperm(1, sat, &c, &mob[0], 0);
double totmob = 0.0; double totmob = 0.0;
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
mob[p] /= visc[p]; mob[p] /= visc[p];
totmob += mob[p]; totmob += mob[p];
} }
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
produced[p] += (mob[p]/totmob)*flux; produced[p] += (mob[p]/totmob)*flux;
} }
} }
} }
} }
@ -203,9 +203,9 @@ namespace Opm
/// @param[in] s saturation values (for all phases) /// @param[in] s saturation values (for all phases)
/// @param[out] totmob total mobilities. /// @param[out] totmob total mobilities.
void computeTotalMobility(const Opm::IncompPropertiesInterface& props, void computeTotalMobility(const Opm::IncompPropertiesInterface& props,
const std::vector<int>& cells, const std::vector<int>& cells,
const std::vector<double>& s, const std::vector<double>& s,
std::vector<double>& totmob) std::vector<double>& totmob)
{ {
std::vector<double> pmobc; std::vector<double> pmobc;
@ -231,10 +231,10 @@ namespace Opm
/// @param[out] totmob total mobility /// @param[out] totmob total mobility
/// @param[out] omega fractional-flow weighted fluid densities. /// @param[out] omega fractional-flow weighted fluid densities.
void computeTotalMobilityOmega(const Opm::IncompPropertiesInterface& props, void computeTotalMobilityOmega(const Opm::IncompPropertiesInterface& props,
const std::vector<int>& cells, const std::vector<int>& cells,
const std::vector<double>& s, const std::vector<double>& s,
std::vector<double>& totmob, std::vector<double>& totmob,
std::vector<double>& omega) std::vector<double>& omega)
{ {
std::vector<double> pmobc; std::vector<double> pmobc;
@ -331,33 +331,33 @@ namespace Opm
/// (+) positive inflow of first phase (water) /// (+) positive inflow of first phase (water)
/// (-) negative total outflow of both phases /// (-) negative total outflow of both phases
void computeTransportSource(const UnstructuredGrid& grid, void computeTransportSource(const UnstructuredGrid& grid,
const std::vector<double>& src, const std::vector<double>& src,
const std::vector<double>& faceflux, const std::vector<double>& faceflux,
const double inflow_frac, const double inflow_frac,
const Wells* wells, const Wells* wells,
const std::vector<double>& well_perfrates, const std::vector<double>& well_perfrates,
std::vector<double>& transport_src) std::vector<double>& transport_src)
{ {
int nc = grid.number_of_cells; int nc = grid.number_of_cells;
transport_src.resize(nc); transport_src.resize(nc);
// Source term and boundary contributions. // Source term and boundary contributions.
for (int c = 0; c < nc; ++c) { for (int c = 0; c < nc; ++c) {
transport_src[c] = 0.0; transport_src[c] = 0.0;
transport_src[c] += src[c] > 0.0 ? inflow_frac*src[c] : src[c]; 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) { for (int hf = grid.cell_facepos[c]; hf < grid.cell_facepos[c + 1]; ++hf) {
int f = grid.cell_faces[hf]; int f = grid.cell_faces[hf];
const int* f2c = &grid.face_cells[2*f]; const int* f2c = &grid.face_cells[2*f];
double bdy_influx = 0.0; double bdy_influx = 0.0;
if (f2c[0] == c && f2c[1] == -1) { if (f2c[0] == c && f2c[1] == -1) {
bdy_influx = -faceflux[f]; bdy_influx = -faceflux[f];
} else if (f2c[0] == -1 && f2c[1] == c) { } else if (f2c[0] == -1 && f2c[1] == c) {
bdy_influx = faceflux[f]; bdy_influx = faceflux[f];
} }
if (bdy_influx != 0.0) { if (bdy_influx != 0.0) {
transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx; transport_src[c] += bdy_influx > 0.0 ? inflow_frac*bdy_influx : bdy_influx;
} }
} }
} }
// Well contributions. // Well contributions.
if (wells) { if (wells) {
@ -392,52 +392,52 @@ namespace Opm
/// @param[in] face_flux signed per-face fluxes /// @param[in] face_flux signed per-face fluxes
/// @param[out] cell_velocity the estimated velocities. /// @param[out] cell_velocity the estimated velocities.
void estimateCellVelocity(const UnstructuredGrid& grid, void estimateCellVelocity(const UnstructuredGrid& grid,
const std::vector<double>& face_flux, const std::vector<double>& face_flux,
std::vector<double>& cell_velocity) std::vector<double>& cell_velocity)
{ {
const int dim = grid.dimensions; const int dim = grid.dimensions;
cell_velocity.clear(); cell_velocity.clear();
cell_velocity.resize(grid.number_of_cells*dim, 0.0); cell_velocity.resize(grid.number_of_cells*dim, 0.0);
for (int face = 0; face < grid.number_of_faces; ++face) { for (int face = 0; face < grid.number_of_faces; ++face) {
int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] }; int c[2] = { grid.face_cells[2*face], grid.face_cells[2*face + 1] };
const double* fc = &grid.face_centroids[face*dim]; const double* fc = &grid.face_centroids[face*dim];
double flux = face_flux[face]; double flux = face_flux[face];
for (int i = 0; i < 2; ++i) { for (int i = 0; i < 2; ++i) {
if (c[i] >= 0) { if (c[i] >= 0) {
const double* cc = &grid.cell_centroids[c[i]*dim]; const double* cc = &grid.cell_centroids[c[i]*dim];
for (int d = 0; d < dim; ++d) { for (int d = 0; d < dim; ++d) {
double v_contrib = fc[d] - cc[d]; double v_contrib = fc[d] - cc[d];
v_contrib *= flux/grid.cell_volumes[c[i]]; v_contrib *= flux/grid.cell_volumes[c[i]];
cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib; cell_velocity[c[i]*dim + d] += (i == 0) ? v_contrib : -v_contrib;
} }
} }
} }
} }
} }
/// Extract a vector of water saturations from a vector of /// Extract a vector of water saturations from a vector of
/// interleaved water and oil saturations. /// interleaved water and oil saturations.
void toWaterSat(const std::vector<double>& sboth, void toWaterSat(const std::vector<double>& sboth,
std::vector<double>& sw) std::vector<double>& sw)
{ {
int num = sboth.size()/2; int num = sboth.size()/2;
sw.resize(num); sw.resize(num);
for (int i = 0; i < num; ++i) { for (int i = 0; i < num; ++i) {
sw[i] = sboth[2*i]; sw[i] = sboth[2*i];
} }
} }
/// Make a a vector of interleaved water and oil saturations from /// Make a a vector of interleaved water and oil saturations from
/// a vector of water saturations. /// a vector of water saturations.
void toBothSat(const std::vector<double>& sw, void toBothSat(const std::vector<double>& sw,
std::vector<double>& sboth) std::vector<double>& sboth)
{ {
int num = sw.size(); int num = sw.size();
sboth.resize(2*num); sboth.resize(2*num);
for (int i = 0; i < num; ++i) { for (int i = 0; i < num; ++i) {
sboth[2*i] = sw[i]; sboth[2*i] = sw[i];
sboth[2*i + 1] = 1.0 - sw[i]; sboth[2*i + 1] = 1.0 - sw[i];
} }
} }
@ -450,30 +450,30 @@ namespace Opm
if (np != 2) { if (np != 2) {
THROW("wellsToSrc() requires a 2 phase case."); THROW("wellsToSrc() requires a 2 phase case.");
} }
src.resize(num_cells); src.resize(num_cells);
for (int w = 0; w < wells.number_of_wells; ++w) { for (int w = 0; w < wells.number_of_wells; ++w) {
const int cur = wells.ctrls[w]->current; const int cur = wells.ctrls[w]->current;
if (wells.ctrls[w]->num != 1) { if (wells.ctrls[w]->num != 1) {
MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored."); MESSAGE("In wellsToSrc(): well has more than one control, all but current control will be ignored.");
} }
if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) { if (wells.ctrls[w]->type[cur] != RESERVOIR_RATE) {
THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE."); THROW("In wellsToSrc(): well is something other than RESERVOIR_RATE.");
} }
if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) { if (wells.well_connpos[w+1] - wells.well_connpos[w] != 1) {
THROW("In wellsToSrc(): well has multiple perforations."); THROW("In wellsToSrc(): well has multiple perforations.");
} }
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
if (wells.ctrls[w]->distr[np*cur + p] != 1.0) { if (wells.ctrls[w]->distr[np*cur + p] != 1.0) {
THROW("In wellsToSrc(): well not controlled on total rate."); 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) { if (wells.type[w] == INJECTOR) {
flow *= wells.comp_frac[np*w + 0]; // Obtaining water rate for inflow source. flow *= wells.comp_frac[np*w + 0]; // Obtaining water rate for inflow source.
} }
const double cell = wells.well_cells[wells.well_connpos[w]]; const double cell = wells.well_cells[wells.well_connpos[w]];
src[cell] = flow; src[cell] = flow;
} }
} }
@ -593,8 +593,10 @@ namespace Opm
{ {
int nw = well_bhp.size(); int nw = well_bhp.size();
ASSERT(nw == wells.number_of_wells); ASSERT(nw == wells.number_of_wells);
if (props.numPhases() != 2) { int np = props.numPhases();
THROW("WellReport for now assumes two phase flow."); const int max_np = 3;
if (np > max_np) {
THROW("WellReport for now assumes #phases <= " << max_np);
} }
const double* visc = props.viscosity(); const double* visc = props.viscosity();
std::vector<double> data_now; std::vector<double> data_now;
@ -605,7 +607,8 @@ namespace Opm
double well_rate_total = 0.0; double well_rate_total = 0.0;
double well_rate_water = 0.0; double well_rate_water = 0.0;
for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { 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; well_rate_total += perf_rate;
if (perf_rate > 0.0) { if (perf_rate > 0.0) {
// Injection. // Injection.
@ -613,11 +616,14 @@ namespace Opm
} else { } else {
// Production. // Production.
const int cell = wells.well_cells[perf]; const int cell = wells.well_cells[perf];
double mob[2]; double mob[max_np];
props.relperm(1, &saturation[2*cell], &cell, mob, 0); props.relperm(1, &saturation[2*cell], &cell, mob, 0);
mob[0] /= visc[0]; double tmob = 0;
mob[1] /= visc[1]; for(int i = 0; i < np; ++i) {
const double fracflow = mob[0]/(mob[0] + mob[1]); mob[i] /= visc[i];
tmob += mob[i];
}
const double fracflow = mob[0]/tmob;
well_rate_water += perf_rate*fracflow; well_rate_water += perf_rate*fracflow;
} }
} }
@ -646,8 +652,10 @@ namespace Opm
// TODO: refactor, since this is almost identical to the other push(). // TODO: refactor, since this is almost identical to the other push().
int nw = well_bhp.size(); int nw = well_bhp.size();
ASSERT(nw == wells.number_of_wells); ASSERT(nw == wells.number_of_wells);
if (props.numPhases() != 2) { int np = props.numPhases();
THROW("WellReport for now assumes two phase flow."); const int max_np = 3;
if (np > max_np) {
THROW("WellReport for now assumes #phases <= " << max_np);
} }
std::vector<double> data_now; std::vector<double> data_now;
data_now.reserve(1 + 3*nw); data_now.reserve(1 + 3*nw);
@ -657,7 +665,8 @@ namespace Opm
double well_rate_total = 0.0; double well_rate_total = 0.0;
double well_rate_water = 0.0; double well_rate_water = 0.0;
for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) { 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; well_rate_total += perf_rate;
if (perf_rate > 0.0) { if (perf_rate > 0.0) {
// Injection. // Injection.
@ -665,13 +674,16 @@ namespace Opm
} else { } else {
// Production. // Production.
const int cell = wells.well_cells[perf]; const int cell = wells.well_cells[perf];
double mob[2]; double mob[max_np];
props.relperm(1, &s[2*cell], &cell, mob, 0); props.relperm(1, &s[np*cell], &cell, mob, 0);
double visc[2]; double visc[max_np];
props.viscosity(1, &p[cell], &z[2*cell], &cell, visc, 0); props.viscosity(1, &p[cell], &z[np*cell], &cell, visc, 0);
mob[0] /= visc[0]; double tmob = 0;
mob[1] /= visc[1]; for(int i = 0; i < np; ++i) {
const double fracflow = mob[0]/(mob[0] + mob[1]); mob[i] /= visc[i];
tmob += mob[i];
}
const double fracflow = mob[0]/(tmob);
well_rate_water += perf_rate*fracflow; well_rate_water += perf_rate*fracflow;
} }
} }

6
tests/Makefile.am
View File

@ -12,6 +12,8 @@ noinst_PROGRAMS = \
bo_resprop_test \ bo_resprop_test \
monotcubicinterpolator_test \ monotcubicinterpolator_test \
param_test \ param_test \
pvt_test \
relperm_test \
sparsetable_test \ sparsetable_test \
sparsevector_test \ sparsevector_test \
test_cartgrid \ test_cartgrid \
@ -32,6 +34,10 @@ monotcubicinterpolator_test_SOURCES = monotcubicinterpolator_test.cpp
param_test_SOURCES = param_test.cpp param_test_SOURCES = param_test.cpp
param_test_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) 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_SOURCES = sparsetable_test.cpp
sparsetable_test_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB) sparsetable_test_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB)

146
tests/pvt_test.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/eclipse/EclipseGridInspector.hpp>
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/grid.h>
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <iterator>
#include <string>
#include <vector>
int main(int argc, char** argv)
{
using namespace std;
// Parameters.
Opm::parameter::ParameterGroup param(argc, argv);
// Parser.
std::string ecl_file = param.get<std::string>("deck_filename");
std::string input_file = param.get<std::string>("input_filename");
std::string matrix_output = param.get<std::string>("matrix_output");
std::string rs_output = param.get<std::string>("rs_output");
std::string b_output = param.get<std::string>("b_output");
std::string mu_output = param.get<std::string>("mu_output");
Opm::EclipseGridParser deck(ecl_file);
UnstructuredGrid grid;
grid.number_of_cells = 1;
grid.global_cell = NULL;
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<double>(aos, " "));
std::copy(dA, dA + np*np*n, std::ostream_iterator<double>(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<double>(muos, " "));
//std::copy(dmu, dmu + np*n, std::ostream_iterator<double>(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<double>(bos, " "));
std::copy(dbdp, dbdp + np*n, std::ostream_iterator<double>(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<double>(rsos, " "));
std::copy(drs, drs + np*n, std::ostream_iterator<double>(rsos, " "));
rsos << std::endl;
}
}
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}

97
tests/relperm_test.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/eclipse/EclipseGridParser.hpp>
#include <opm/core/eclipse/EclipseGridInspector.hpp>
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/grid.h>
#include <algorithm>
#include <iostream>
#include <iterator>
#include <string>
#include <vector>
#include <fstream>
int main(int argc, char** argv)
{
// Parameters.
Opm::parameter::ParameterGroup param(argc, argv);
// Parser.
std::string ecl_file = param.get<std::string>("deck_filename");
std::string input_file = param.get<std::string>("input_filename");
std::string relperm_output = param.get<std::string>("relperm_output");
//std::string relperm_output = param.get<std::string>("relperm_outout");
Opm::EclipseGridParser deck(ecl_file);
UnstructuredGrid grid;
grid.number_of_cells = 1;
grid.global_cell = NULL;
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<double>(kros, " "));
kros << " ";
std::copy(kr, kr + np, std::ostream_iterator<double>(kros, " "));
kros << " ";
std::copy(dkr, dkr + np*np, std::ostream_iterator<double>(kros, " "));
kros << "\n";
}
}
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}