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

Merge remote-tracking branch 'upstream/master' into ert

Browse Source 
Conflicts:
	Makefile.am
	opm/core/grid/cpgpreprocess/preprocess.h
	tests/Makefile.am

This brings ert branch up-to-date with current Github master branch.
This commit is contained in:
Bård Skaflestad
2012-10-12 00:43:51 +02:00
parent 2e99dd48dc 09b2ecc5de
commit 31174cf727
41 changed files with 2888 additions and 386 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@@ -1,6 +1,7 @@
*~
.\#*
*.o
*.mod
*.lo
*.la
.libs
@@ -34,6 +35,7 @@ tutorials/tutorial[1-4]
# Ignoring executables
*_test
examples/compute_tof
examples/scaneclipsedeck
examples/spu_2p
examples/reorder-qfs
@@ -42,8 +44,10 @@ examples/sim_2p_incomp_reorder
examples/sim_2p_comp_reorder
examples/sim_wateroil
examples/wells_example
tests/test_agmg
tests/test_cfs_tpfa
tests/test_jacsys
tests/test_read_grid
tests/test_readvector
tests/test_sf2p
tests/bo_fluid_p_and_z_deps
@@ -53,4 +57,5 @@ tests/test_column_extract
tests/test_lapack
tests/test_read_vag
tests/test_readpolymer
tests/test_wells
tests/test_writeVtkData

429
Makefile.am
View File

@@ -11,23 +11,27 @@ lib_LTLIBRARIES = lib/libopmcore.la
# ----------------------------------------------------------------------
# Build-time flags needed to build libopmcore.la
AM_CPPFLAGS = \
AM_CPPFLAGS = \
$(ERT_CPPFLAGS) \
$(OPM_BOOST_CPPFLAGS)
$(OPM_BOOST_CPPFLAGS) \
${SUPERLU_CPPFLAGS}
# ----------------------------------------------------------------------
# Link-time flags needed both to successfully link the library and to
# (transitively) convey inter-library dependency information.
lib_libopmcore_la_LDFLAGS = \
$(ERT_LDFLAGS) \
lib_libopmcore_la_LDFLAGS = \
-R $(OPM_BOOST_LIBDIR) \
$(OPM_BOOST_LDFLAGS) \
$(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB) \
$(BOOST_DATE_TIME_LIB) \
$(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
${SUPERLU_LDFLAGS}
lib_libopmcore_la_LIBADD = \
$(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB) \
$(BOOST_DATE_TIME_LIB) \
$(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
$(ERT_LIBS) \
$(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS)
$(LAPACK_LIBS) ${SUPERLU_LIBS} $(BLAS_LIBS) $(LIBS)
# ----------------------------------------------------------------------
# Library constituents. SOURCES followed by HEADERS.
@@ -56,80 +60,82 @@ opm/core/fluid/SaturationPropsFromDeck.cpp \
opm/core/fluid/SatFuncGwseg.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/SinglePvtDeadSpline.cpp \
opm/core/fluid/blackoil/SinglePvtInterface.cpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.cpp \
opm/core/grid.c \
opm/core/grid/cart_grid.c \
opm/core/grid/cornerpoint_grid.c \
opm/core/grid/cpgpreprocess/facetopology.c \
opm/core/grid/cpgpreprocess/geometry.c \
opm/core/grid/cpgpreprocess/preprocess.c \
opm/core/grid/cpgpreprocess/uniquepoints.c \
opm/core/linalg/LinearSolverFactory.cpp \
opm/core/linalg/LinearSolverInterface.cpp \
opm/core/linalg/sparse_sys.c \
opm/core/newwells.c \
opm/core/pressure/CompressibleTpfa.cpp \
opm/core/pressure/FlowBCManager.cpp \
opm/core/pressure/IncompTpfa.cpp \
opm/core/pressure/cfsh.c \
opm/core/pressure/flow_bc.c \
opm/core/pressure/fsh.c \
opm/core/pressure/fsh_common_impl.c \
opm/core/pressure/ifsh.c \
opm/core/pressure/mimetic/hybsys.c \
opm/core/pressure/mimetic/hybsys_global.c \
opm/core/pressure/mimetic/mimetic.c \
opm/core/pressure/msmfem/coarse_conn.c \
opm/core/pressure/msmfem/coarse_sys.c \
opm/core/pressure/msmfem/dfs.c \
opm/core/pressure/msmfem/hash_set.c \
opm/core/pressure/msmfem/ifsh_ms.c \
opm/core/pressure/msmfem/partition.c \
opm/core/pressure/tpfa/cfs_tpfa.c \
opm/core/pressure/tpfa/cfs_tpfa_residual.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/compr_source.c \
opm/core/pressure/tpfa/ifs_tpfa.c \
opm/core/pressure/tpfa/trans_tpfa.c \
opm/core/pressure/well.c \
opm/core/simulator/SimulatorCompressibleTwophase.cpp \
opm/core/simulator/SimulatorIncompTwophase.cpp \
opm/core/simulator/SimulatorReport.cpp \
opm/core/simulator/SimulatorTimer.cpp \
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \
opm/core/transport/reorder/TransportModelInterface.cpp \
opm/core/transport/reorder/TransportModelTwophase.cpp \
opm/core/transport/reorder/nlsolvers.c \
opm/core/transport/reorder/reordersequence.cpp \
opm/core/transport/reorder/tarjan.c \
opm/core/transport/spu_explicit.c \
opm/core/transport/spu_implicit.c \
opm/core/transport/transport_source.c \
opm/core/utility/MonotCubicInterpolator.cpp \
opm/core/utility/StopWatch.cpp \
opm/core/utility/miscUtilities.cpp \
opm/core/utility/miscUtilitiesBlackoil.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/writeVtkData.cpp \
opm/core/vag_format/vag.cpp \
opm/core/wells/InjectionSpecification.cpp \
opm/core/wells/ProductionSpecification.cpp \
opm/core/wells/WellCollection.cpp \
opm/core/wells/WellsGroup.cpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp \
opm/core/fluid/blackoil/SinglePvtDead.cpp \
opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp \
opm/core/fluid/blackoil/SinglePvtInterface.cpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.cpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.cpp \
opm/core/grid.c \
opm/core/grid/cart_grid.c \
opm/core/grid/cornerpoint_grid.c \
opm/core/grid/cpgpreprocess/facetopology.c \
opm/core/grid/cpgpreprocess/geometry.c \
opm/core/grid/cpgpreprocess/preprocess.c \
opm/core/grid/cpgpreprocess/uniquepoints.c \
opm/core/linalg/LinearSolverFactory.cpp \
opm/core/linalg/LinearSolverInterface.cpp \
opm/core/linalg/sparse_sys.c \
opm/core/newwells.c \
opm/core/pressure/CompressibleTpfa.cpp \
opm/core/pressure/FlowBCManager.cpp \
opm/core/pressure/IncompTpfa.cpp \
opm/core/pressure/cfsh.c \
opm/core/pressure/flow_bc.c \
opm/core/pressure/fsh.c \
opm/core/pressure/fsh_common_impl.c \
opm/core/pressure/ifsh.c \
opm/core/pressure/mimetic/hybsys.c \
opm/core/pressure/mimetic/hybsys_global.c \
opm/core/pressure/mimetic/mimetic.c \
opm/core/pressure/msmfem/coarse_conn.c \
opm/core/pressure/msmfem/coarse_sys.c \
opm/core/pressure/msmfem/dfs.c \
opm/core/pressure/msmfem/hash_set.c \
opm/core/pressure/msmfem/ifsh_ms.c \
opm/core/pressure/msmfem/partition.c \
opm/core/pressure/tpfa/cfs_tpfa.c \
opm/core/pressure/tpfa/cfs_tpfa_residual.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/compr_source.c \
opm/core/pressure/tpfa/ifs_tpfa.c \
opm/core/pressure/tpfa/trans_tpfa.c \
opm/core/pressure/well.c \
opm/core/simulator/SimulatorCompressibleTwophase.cpp \
opm/core/simulator/SimulatorIncompTwophase.cpp \
opm/core/simulator/SimulatorReport.cpp \
opm/core/simulator/SimulatorTimer.cpp \
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp \
opm/core/transport/reorder/TransportModelInterface.cpp \
opm/core/transport/reorder/TransportModelTracerTof.cpp \
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp \
opm/core/transport/reorder/TransportModelTwophase.cpp \
opm/core/transport/reorder/nlsolvers.c \
opm/core/transport/reorder/reordersequence.cpp \
opm/core/transport/reorder/tarjan.c \
opm/core/transport/spu_explicit.c \
opm/core/transport/spu_implicit.c \
opm/core/transport/transport_source.c \
opm/core/utility/MonotCubicInterpolator.cpp \
opm/core/utility/StopWatch.cpp \
opm/core/utility/miscUtilities.cpp \
opm/core/utility/miscUtilitiesBlackoil.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/writeVtkData.cpp \
opm/core/vag_format/vag.cpp \
opm/core/wells/InjectionSpecification.cpp \
opm/core/wells/ProductionSpecification.cpp \
opm/core/wells/WellCollection.cpp \
opm/core/wells/WellsGroup.cpp \
opm/core/wells/WellsManager.cpp
nobase_include_HEADERS = \
@@ -155,137 +161,140 @@ opm/core/fluid/RockFromDeck.hpp \
opm/core/fluid/SatFuncGwseg.hpp \
opm/core/fluid/SatFuncStone2.hpp \
opm/core/fluid/SatFuncSimple.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \
opm/core/fluid/SaturationPropsFromDeck_impl.hpp \
opm/core/fluid/SaturationPropsInterface.hpp \
opm/core/fluid/SimpleFluid2p.hpp \
opm/core/fluid/blackoil/BlackoilPhases.hpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \
opm/core/fluid/blackoil/SinglePvtDead.hpp \
opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp \
opm/core/fluid/blackoil/SinglePvtInterface.hpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.hpp \
opm/core/fluid/blackoil/phaseUsageFromDeck.hpp \
opm/core/grid.h \
opm/core/grid/cart_grid.h \
opm/core/grid/cornerpoint_grid.h \
opm/core/grid/cpgpreprocess/facetopology.h \
opm/core/grid/cpgpreprocess/geometry.h \
opm/core/grid/cpgpreprocess/grdecl.h \
opm/core/grid/cpgpreprocess/preprocess.h \
opm/core/grid/cpgpreprocess/uniquepoints.h \
opm/core/linalg/LinearSolverFactory.hpp \
opm/core/linalg/LinearSolverInterface.hpp \
opm/core/linalg/blas_lapack.h \
opm/core/linalg/sparse_sys.h \
opm/core/newwells.h \
opm/core/pressure/CompressibleTpfa.hpp \
opm/core/pressure/FlowBCManager.hpp \
opm/core/pressure/HybridPressureSolver.hpp \
opm/core/pressure/IncompTpfa.hpp \
opm/core/pressure/TPFACompressiblePressureSolver.hpp \
opm/core/pressure/TPFAPressureSolver.hpp \
opm/core/pressure/flow_bc.h \
opm/core/pressure/fsh.h \
opm/core/pressure/fsh_common_impl.h \
opm/core/pressure/mimetic/hybsys.h \
opm/core/pressure/mimetic/hybsys_global.h \
opm/core/pressure/mimetic/mimetic.h \
opm/core/pressure/msmfem/coarse_conn.h \
opm/core/pressure/msmfem/coarse_sys.h \
opm/core/pressure/msmfem/dfs.h \
opm/core/pressure/msmfem/hash_set.h \
opm/core/pressure/msmfem/ifsh_ms.h \
opm/core/pressure/msmfem/partition.h \
opm/core/pressure/tpfa/cfs_tpfa.h \
opm/core/pressure/tpfa/cfs_tpfa_residual.h \
opm/core/pressure/tpfa/compr_bc.h \
opm/core/pressure/tpfa/compr_quant.h \
opm/core/pressure/tpfa/compr_quant_general.h \
opm/core/pressure/tpfa/compr_source.h \
opm/core/pressure/tpfa/ifs_tpfa.h \
opm/core/pressure/tpfa/trans_tpfa.h \
opm/core/simulator/BlackoilState.hpp \
opm/core/simulator/SimulatorCompressibleTwophase.hpp \
opm/core/simulator/SimulatorReport.hpp \
opm/core/simulator/SimulatorIncompTwophase.hpp \
opm/core/simulator/SimulatorTimer.hpp \
opm/core/simulator/TwophaseState.hpp \
opm/core/simulator/WellState.hpp \
opm/core/transport/CSRMatrixBlockAssembler.hpp \
opm/core/transport/CSRMatrixUmfpackSolver.hpp \
opm/core/transport/GravityColumnSolver.hpp \
opm/core/transport/GravityColumnSolver_impl.hpp \
opm/core/transport/ImplicitAssembly.hpp \
opm/core/transport/ImplicitTransport.hpp \
opm/core/transport/JacobianSystem.hpp \
opm/core/transport/NormSupport.hpp \
opm/core/transport/SimpleFluid2pWrapper.hpp \
opm/core/transport/SinglePointUpwindTwoPhase.hpp \
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp \
opm/core/transport/reorder/TransportModelInterface.hpp \
opm/core/transport/reorder/TransportModelTwophase.hpp \
opm/core/transport/reorder/nlsolvers.h \
opm/core/transport/reorder/reordersequence.h \
opm/core/transport/reorder/tarjan.h \
opm/core/transport/spu_explicit.h \
opm/core/transport/spu_implicit.h \
opm/core/transport/transport_source.h \
opm/core/utility/Average.hpp \
opm/core/utility/ColumnExtract.hpp \
opm/core/utility/DataMap.hpp \
opm/core/utility/ErrorMacros.hpp \
opm/core/utility/Factory.hpp \
opm/core/utility/MonotCubicInterpolator.hpp \
opm/core/utility/NonuniformTableLinear.hpp \
opm/core/utility/RootFinders.hpp \
opm/core/utility/SparseTable.hpp \
opm/core/utility/SparseVector.hpp \
opm/core/utility/StopWatch.hpp \
opm/core/utility/UniformTableLinear.hpp \
opm/core/utility/Units.hpp \
opm/core/utility/buildUniformMonotoneTable.hpp \
opm/core/utility/initState.hpp \
opm/core/utility/initState_impl.hpp \
opm/core/utility/linInt.hpp \
opm/core/utility/linearInterpolation.hpp \
opm/core/utility/miscUtilities.hpp \
opm/core/utility/miscUtilitiesBlackoil.hpp \
opm/core/utility/parameters/Parameter.hpp \
opm/core/utility/parameters/ParameterGroup.hpp \
opm/core/utility/parameters/ParameterGroup_impl.hpp \
opm/core/utility/parameters/ParameterMapItem.hpp \
opm/core/utility/parameters/ParameterRequirement.hpp \
opm/core/utility/parameters/ParameterStrings.hpp \
opm/core/utility/parameters/ParameterTools.hpp \
opm/core/utility/parameters/ParameterXML.hpp \
opm/core/utility/parameters/tinyxml/tinystr.h \
opm/core/utility/parameters/tinyxml/tinyxml.h \
opm/core/utility/writeVtkData.hpp \
opm/core/vag_format/vag.hpp \
opm/core/well.h \
opm/core/wells/InjectionSpecification.hpp \
opm/core/wells/ProductionSpecification.hpp \
opm/core/wells/WellCollection.hpp \
opm/core/wells/WellsGroup.hpp \
opm/core/fluid/SaturationPropsBasic.hpp \
opm/core/fluid/SaturationPropsFromDeck.hpp \
opm/core/fluid/SaturationPropsFromDeck_impl.hpp \
opm/core/fluid/SaturationPropsInterface.hpp \
opm/core/fluid/SimpleFluid2p.hpp \
opm/core/fluid/blackoil/BlackoilPhases.hpp \
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp \
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp \
opm/core/fluid/blackoil/SinglePvtDead.hpp \
opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp \
opm/core/fluid/blackoil/SinglePvtInterface.hpp \
opm/core/fluid/blackoil/SinglePvtLiveGas.hpp \
opm/core/fluid/blackoil/SinglePvtLiveOil.hpp \
opm/core/fluid/blackoil/phaseUsageFromDeck.hpp \
opm/core/grid.h \
opm/core/grid/cart_grid.h \
opm/core/grid/cornerpoint_grid.h \
opm/core/grid/cpgpreprocess/facetopology.h \
opm/core/grid/cpgpreprocess/geometry.h \
opm/core/grid/cpgpreprocess/grdecl.h \
opm/core/grid/cpgpreprocess/preprocess.h \
opm/core/grid/cpgpreprocess/uniquepoints.h \
opm/core/linalg/LinearSolverFactory.hpp \
opm/core/linalg/LinearSolverInterface.hpp \
opm/core/linalg/blas_lapack.h \
opm/core/linalg/sparse_sys.h \
opm/core/newwells.h \
opm/core/pressure/CompressibleTpfa.hpp \
opm/core/pressure/FlowBCManager.hpp \
opm/core/pressure/HybridPressureSolver.hpp \
opm/core/pressure/IncompTpfa.hpp \
opm/core/pressure/TPFACompressiblePressureSolver.hpp \
opm/core/pressure/TPFAPressureSolver.hpp \
opm/core/pressure/flow_bc.h \
opm/core/pressure/fsh.h \
opm/core/pressure/fsh_common_impl.h \
opm/core/pressure/mimetic/hybsys.h \
opm/core/pressure/mimetic/hybsys_global.h \
opm/core/pressure/mimetic/mimetic.h \
opm/core/pressure/msmfem/coarse_conn.h \
opm/core/pressure/msmfem/coarse_sys.h \
opm/core/pressure/msmfem/dfs.h \
opm/core/pressure/msmfem/hash_set.h \
opm/core/pressure/msmfem/ifsh_ms.h \
opm/core/pressure/msmfem/partition.h \
opm/core/pressure/tpfa/cfs_tpfa.h \
opm/core/pressure/tpfa/cfs_tpfa_residual.h \
opm/core/pressure/tpfa/compr_bc.h \
opm/core/pressure/tpfa/compr_quant.h \
opm/core/pressure/tpfa/compr_quant_general.h \
opm/core/pressure/tpfa/compr_source.h \
opm/core/pressure/tpfa/ifs_tpfa.h \
opm/core/pressure/tpfa/trans_tpfa.h \
opm/core/simulator/BlackoilState.hpp \
opm/core/simulator/SimulatorCompressibleTwophase.hpp \
opm/core/simulator/SimulatorReport.hpp \
opm/core/simulator/SimulatorIncompTwophase.hpp \
opm/core/simulator/SimulatorTimer.hpp \
opm/core/simulator/TwophaseState.hpp \
opm/core/simulator/WellState.hpp \
opm/core/transport/CSRMatrixBlockAssembler.hpp \
opm/core/transport/CSRMatrixUmfpackSolver.hpp \
opm/core/transport/GravityColumnSolver.hpp \
opm/core/transport/GravityColumnSolver_impl.hpp \
opm/core/transport/ImplicitAssembly.hpp \
opm/core/transport/ImplicitTransport.hpp \
opm/core/transport/JacobianSystem.hpp \
opm/core/transport/NormSupport.hpp \
opm/core/transport/SimpleFluid2pWrapper.hpp \
opm/core/transport/SinglePointUpwindTwoPhase.hpp \
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp \
opm/core/transport/reorder/TransportModelInterface.hpp \
opm/core/transport/reorder/TransportModelTracerTof.hpp \
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp \
opm/core/transport/reorder/TransportModelTwophase.hpp \
opm/core/transport/reorder/nlsolvers.h \
opm/core/transport/reorder/reordersequence.h \
opm/core/transport/reorder/tarjan.h \
opm/core/transport/spu_explicit.h \
opm/core/transport/spu_implicit.h \
opm/core/transport/transport_source.h \
opm/core/utility/Average.hpp \
opm/core/utility/ColumnExtract.hpp \
opm/core/utility/DataMap.hpp \
opm/core/utility/ErrorMacros.hpp \
opm/core/utility/Factory.hpp \
opm/core/utility/MonotCubicInterpolator.hpp \
opm/core/utility/NonuniformTableLinear.hpp \
opm/core/utility/RootFinders.hpp \
opm/core/utility/SparseTable.hpp \
opm/core/utility/SparseVector.hpp \
opm/core/utility/StopWatch.hpp \
opm/core/utility/UniformTableLinear.hpp \
opm/core/utility/Units.hpp \
opm/core/utility/buildUniformMonotoneTable.hpp \
opm/core/utility/initState.hpp \
opm/core/utility/initState_impl.hpp \
opm/core/utility/linInt.hpp \
opm/core/utility/linearInterpolation.hpp \
opm/core/utility/miscUtilities.hpp \
opm/core/utility/miscUtilitiesBlackoil.hpp \
opm/core/utility/parameters/Parameter.hpp \
opm/core/utility/parameters/ParameterGroup.hpp \
opm/core/utility/parameters/ParameterGroup_impl.hpp \
opm/core/utility/parameters/ParameterMapItem.hpp \
opm/core/utility/parameters/ParameterRequirement.hpp \
opm/core/utility/parameters/ParameterStrings.hpp \
opm/core/utility/parameters/ParameterTools.hpp \
opm/core/utility/parameters/ParameterXML.hpp \
opm/core/utility/parameters/tinyxml/tinystr.h \
opm/core/utility/parameters/tinyxml/tinyxml.h \
opm/core/utility/writeVtkData.hpp \
opm/core/vag_format/vag.hpp \
opm/core/well.h \
opm/core/wells/InjectionSpecification.hpp \
opm/core/wells/ProductionSpecification.hpp \
opm/core/wells/WellCollection.hpp \
opm/core/wells/WellsGroup.hpp \
opm/core/wells/WellsManager.hpp
# ----------------------------------------------------------------------
# Optional library constituents.
if UMFPACK
lib_libopmcore_la_SOURCES += \
opm/core/linalg/call_umfpack.c \
lib_libopmcore_la_SOURCES += \
opm/core/linalg/call_umfpack.c \
opm/core/linalg/LinearSolverUmfpack.cpp
nobase_include_HEADERS += \
opm/core/linalg/call_umfpack.h \
nobase_include_HEADERS += \
opm/core/linalg/call_umfpack.h \
opm/core/linalg/LinearSolverUmfpack.hpp
endif
if HAVE_ERT
lib_libopmcore_la_SOURCES += \
opm/core/utility/writeECLData.cpp
@@ -296,26 +305,26 @@ endif
if DUNE_ISTL
lib_libopmcore_la_SOURCES += \
lib_libopmcore_la_SOURCES += \
opm/core/linalg/LinearSolverIstl.cpp
nobase_include_HEADERS += \
nobase_include_HEADERS += \
opm/core/linalg/LinearSolverIstl.hpp
endif
if BUILD_AGMG
nodist_lib_libopmcore_la_SOURCES += \
$(AGMG_SRCDIR)/dagmg.f90 \
nodist_lib_libopmcore_la_SOURCES += \
$(AGMG_SRCDIR)/dagmg.f90 \
$(AGMG_SRCDIR)/dagmg_mumps.f90
lib_libopmcore_la_SOURCES += \
lib_libopmcore_la_SOURCES += \
opm/core/linalg/LinearSolverAGMG.cpp
nobase_include_HEADERS += \
nobase_include_HEADERS += \
opm/core/linalg/LinearSolverAGMG.hpp
lib_libopmcore_la_LDFLAGS += \
lib_libopmcore_la_LDFLAGS += \
$(FCLIBS)
endif

15
README
View File

@@ -88,10 +88,25 @@ If you want to contribute, fork OPM/opm-core on github.
BUILDING
--------
There are two ways to build the opm-core library:
1. As a stand-alone library.
cd opm-core
autoreconf -i
./configure
make
If you want to install the library:
make install
or (if installing to /usr/local or similar)
sudo make install
2. As a dune module.
- Put the opm-core directory in the same directory
as the other dune modules to be built (e.g. dune-commmon,
dune-grid).
- Run dunecontrol as normal. For more information on
the dune build system, see
http://www.dune-project.org/doc/installation-notes.html
DOCUMENTATION

29
examples/Makefile.am
View File

@@ -3,12 +3,24 @@ ERT_INCLUDE_PATH = $(ERT_ROOT)/include
AM_CPPFLAGS = \
-I$(top_srcdir) \
$(OPM_BOOST_CPPFLAGS) \
-I$(ERT_INCLUDE_PATH)
-I$(ERT_INCLUDE_PATH) \
$(OPM_BOOST_CPPFLAGS)
# All targets link to the library
LDADD = \
$(top_builddir)/lib/libopmcore.la \
$(top_builddir)/lib/libopmcore.la
# Convenience definition for targets that use Boost.Filesystem directly.
# While libopmcore depends on (and references) Boost.Filesystem (through
# the $(BOOST_FILESYSTEM_LIB) macro) this indirect dependency is not
# sufficient to satisfy the requirements of targets that use the indirect
# libraries directly.
#
# Additional details at
# https://fedoraproject.org/wiki/UnderstandingDSOLinkChange
#
LINK_BOOST_FILESYSTEM = \
$(OPM_BOOST_LDFLAGS) \
$(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB)
@@ -18,6 +30,7 @@ $(BOOST_SYSTEM_LIB)
# Please keep the list sorted.
noinst_PROGRAMS = \
compute_tof \
refine_wells \
scaneclipsedeck \
sim_2p_comp_reorder \
@@ -32,10 +45,19 @@ wells_example
#
# Please maintain sort order from "noinst_PROGRAMS".
compute_tof_SOURCES = compute_tof.cpp
compute_tof_LDADD = $(LDADD) $(LINK_BOOST_FILESYSTEM)
refine_wells_SOURCES = refine_wells.cpp
sim_2p_comp_reorder_SOURCES = sim_2p_comp_reorder.cpp
sim_2p_comp_reorder_LDADD = $(LDADD) $(LINK_BOOST_FILESYSTEM)
sim_2p_incomp_reorder_SOURCES = sim_2p_incomp_reorder.cpp
sim_2p_incomp_reorder_LDADD = $(LDADD) $(LINK_BOOST_FILESYSTEM)
sim_wateroil_SOURCES = sim_wateroil.cpp
sim_wateroil_LDADD = $(LDADD) $(LINK_BOOST_FILESYSTEM)
wells_example_SOURCES = wells_example.cpp
# ----------------------------------------------------------------------
@@ -47,5 +69,6 @@ noinst_PROGRAMS += spu_2p
spu_2p_SOURCES = spu_2p.cpp
spu_2p_LDADD = \
$(LDADD) \
$(LINK_BOOST_FILESYSTEM) \
$(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS)
endif

254
examples/compute_tof.cpp Normal file
View File

@@ -0,0 +1,254 @@
/*
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/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
#include <algorithm>
#include <iostream>
#include <vector>
#include <numeric>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
using namespace Opm;
std::cout << "\n================ Test program for incompressible tof computations ===============\n\n";
parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
boost::scoped_ptr<EclipseGridParser> deck;
boost::scoped_ptr<GridManager> grid;
boost::scoped_ptr<IncompPropertiesInterface> props;
boost::scoped_ptr<Opm::WellsManager> wells;
TwophaseState state;
// bool check_well_controls = false;
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
deck.reset(new EclipseGridParser(deck_filename));
// Grid init
grid.reset(new GridManager(*deck));
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
// Wells init.
wells.reset(new Opm::WellsManager(*deck, *grid->c_grid(), props->permeability()));
// Gravity.
gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, *deck, gravity[2], state);
}
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0);
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Wells init.
wells.reset(new Opm::WellsManager());
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
}
// Warn if gravity but no density difference.
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
if (use_gravity) {
if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
}
const double *grav = use_gravity ? &gravity[0] : 0;
// Initialising src
std::vector<double> porevol;
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> src(num_cells, 0.0);
if (use_deck) {
// Do nothing, wells will be the driving force, not source terms.
} else {
const double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
*tot_porevol_init/unit::day;
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
}
// Boundary conditions.
FlowBCManager bcs;
if (param.getDefault("use_pside", false)) {
int pside = param.get<int>("pside");
double pside_pressure = param.get<double>("pside_pressure");
bcs.pressureSide(*grid->c_grid(), FlowBCManager::Side(pside), pside_pressure);
}
// Linear solver.
LinearSolverFactory linsolver(param);
// Pressure solver.
Opm::IncompTpfa psolver(*grid->c_grid(), *props, 0, linsolver,
0.0, 0.0, 0,
grav, wells->c_wells(), src, bcs.c_bcs());
// Choice of tof solver.
bool use_dg = param.getDefault("use_dg", false);
int dg_degree = -1;
if (use_dg) {
dg_degree = param.getDefault("dg_degree", 0);
}
// Write parameters used for later reference.
bool output = param.getDefault("output", true);
std::ofstream epoch_os;
std::string output_dir;
if (output) {
output_dir =
param.getDefault("output_dir", std::string("output"));
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
std::string filename = output_dir + "/epoch_timing.param";
epoch_os.open(filename.c_str(), std::fstream::trunc | std::fstream::out);
// open file to clean it. The file is appended to in SimulatorTwophase
filename = output_dir + "/step_timing.param";
std::fstream step_os(filename.c_str(), std::fstream::trunc | std::fstream::out);
step_os.close();
param.writeParam(output_dir + "/simulation.param");
}
// Init wells.
Opm::WellState well_state;
well_state.init(wells->c_wells(), state);
// Main solvers.
Opm::time::StopWatch pressure_timer;
double ptime = 0.0;
Opm::time::StopWatch transport_timer;
double ttime = 0.0;
Opm::time::StopWatch total_timer;
total_timer.start();
std::cout << "\n\n================ Starting main solvers ===============" << std::endl;
// Solve pressure.
pressure_timer.start();
psolver.solve(1.0, state, well_state);
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
// Process transport sources (to include bdy terms and well flows).
std::vector<double> transport_src;
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_state.perfRates(), transport_src);
// Solve time-of-flight.
std::vector<double> tof;
if (use_dg) {
Opm::TransportModelTracerTofDiscGal tofsolver(*grid->c_grid());
transport_timer.start();
tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], dg_degree, tof);
transport_timer.stop();
} else {
Opm::TransportModelTracerTof tofsolver(*grid->c_grid());
transport_timer.start();
tofsolver.solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
transport_timer.stop();
}
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
total_timer.stop();
// Output.
if (output) {
std::string tof_filename = output_dir + "/tof.txt";
std::ofstream tof_stream(tof_filename.c_str());
std::copy(tof.begin(), tof.end(), std::ostream_iterator<double>(tof_stream, "\n"));
}
std::cout << "\n\n================ End of simulation ===============\n"
<< "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime
<< "\n Transport time: " << ttime << std::endl;
}

43
m4/ax_boost_unit_test_framework.m4
View File

@@ -102,27 +102,32 @@ AC_DEFUN([AX_BOOST_UNIT_TEST_FRAMEWORK],
fi
else
link_unit_test_framework="no"
saved_ldflags="${LDFLAGS}"
saved_ldflags="${LDFLAGS}"
for ax_lib in boost_unit_test_framework-$ax_boost_user_unit_test_framework_lib $ax_boost_user_unit_test_framework_lib ; do
if test "x$link_unit_test_framework" = "xyes"; then
break;
fi
for unittest_library in `ls $BOOSTLIBDIR/lib${ax_lib}.so* $BOOSTLIBDIR/lib${ax_lib}.a* 2>/dev/null` ; do
if test -r $unittest_library ; then
libextension=`echo $unittest_library | sed 's,.*/,,' | sed -e 's;^lib\(boost_unit_test_framework.*\)\.so.*$;\1;' -e 's;^lib\(boost_unit_test_framework.*\)\.a*$;\1;'`
ax_lib=${libextension}
link_unit_test_framework="yes"
else
link_unit_test_framework="no"
fi
#if test "x$link_unit_test_framework" = "xyes"; then
# break;
#fi
#for unittest_library in `ls $BOOSTLIBDIR/lib${ax_lib}.so* $BOOSTLIBDIR/lib${ax_lib}.a* 2>/dev/null` ; do
#if test -r $unittest_library ; then
# libextension=`echo $unittest_library | sed 's,.*/,,' | sed -e 's;^lib\(boost_unit_test_framework.*\)\.so.*$;\1;' -e 's;^lib\(boost_unit_test_framework.*\)\.a*$;\1;'`
# ax_lib=${libextension}
# link_unit_test_framework="yes"
# else
# link_unit_test_framework="no"
# fi
#
# if test "x$link_unit_test_framework" = "xyes"; then
# BOOST_UNIT_TEST_FRAMEWORK_LIB="-l$ax_lib"
# AC_SUBST(BOOST_UNIT_TEST_FRAMEWORK_LIB)
# break
# fi
if test "x$link_unit_test_framework" = "xyes"; then
BOOST_UNIT_TEST_FRAMEWORK_LIB="-l$ax_lib"
AC_SUBST(BOOST_UNIT_TEST_FRAMEWORK_LIB)
break
fi
done
done
# done
AC_CHECK_LIB($ax_lib, [exit],
[BOOST_UNIT_TEST_FRAMEWORK_LIB="-l$ax_lib"; AC_SUBST([BOOST_UNIT_TEST_FRAMEWORK_LIB]) link_unit_test_framework="yes"; break],
[link_unit_test_framework="no"])
done
fi
if test "x$ax_lib" = "x"; then
AC_MSG_ERROR(Could not find a version of the library!)

6
m4/opm_boost_base.m4
View File

@@ -103,6 +103,7 @@ if test "x$want_boost" = "xyes"; then
for ac_boost_path_tmp in $libsubdirs; do
if test -d "$ac_boost_path"/"$ac_boost_path_tmp" ; then
OPM_BOOST_LDFLAGS="-L$ac_boost_path/$ac_boost_path_tmp"
OPM_BOOST_LIBDIR="${ac_boost_path}/${ac_boost_path_tmp}"
break
fi
done
@@ -114,6 +115,7 @@ if test "x$want_boost" = "xyes"; then
done
OPM_BOOST_LDFLAGS="-L$ac_boost_path_tmp/$libsubdir"
OPM_BOOST_CPPFLAGS="-I$ac_boost_path_tmp/include"
OPM_BOOST_LIBDIR="${ac_boost_path_tmp}/${libsubdir}"
break;
fi
done
@@ -123,6 +125,7 @@ if test "x$want_boost" = "xyes"; then
dnl --with-boost-libdir parameter
if test "$ac_boost_lib_path" != ""; then
OPM_BOOST_LDFLAGS="-L$ac_boost_lib_path"
OPM_BOOST_LIBDIR="${ac_boost_lib_path}"
fi
CPPFLAGS_SAVED="$CPPFLAGS"
@@ -191,6 +194,7 @@ if test "x$want_boost" = "xyes"; then
if ls "$best_path/$libsubdir/libboost_"* >/dev/null 2>&1 ; then break; fi
done
OPM_BOOST_LDFLAGS="-L$best_path/$libsubdir"
OPM_BOOST_LIBDIR="$best_path/$libsubdir"
fi
fi
@@ -207,6 +211,7 @@ if test "x$want_boost" = "xyes"; then
AC_MSG_NOTICE(We will use a staged boost library from $BOOST_ROOT)
OPM_BOOST_CPPFLAGS="-I$BOOST_ROOT"
OPM_BOOST_LDFLAGS="-L$BOOST_ROOT/stage/$libsubdir"
OPM_BOOST_LIBDIR="$BOOST_ROOT/stage/$libsubdir"
fi
fi
fi
@@ -248,6 +253,7 @@ if test "x$want_boost" = "xyes"; then
else
AC_SUBST([OPM_BOOST_CPPFLAGS])
AC_SUBST([OPM_BOOST_LDFLAGS])
AC_SUBST([OPM_BOOST_LIBDIR])
AC_DEFINE([OPM_HAVE_BOOST], [1], [define if the Boost library is available])
# execute ACTION-IF-FOUND (if present):
ifelse([$2], , :, [$2])

2
m4/opm_core.m4
View File

@@ -17,8 +17,8 @@ AX_BOOST_SYSTEM
AX_BOOST_DATE_TIME
AX_BOOST_FILESYSTEM
AX_BOOST_UNIT_TEST_FRAMEWORK
AX_DUNE_ISTL
OPM_PATH_SUPERLU
OPM_AGMG
# Checks for header files.

334
m4/opm_superlu.m4 Normal file
View File

@@ -0,0 +1,334 @@
## -*- autoconf -*-
# $Id: superlu.m4 5908 2010-02-23 16:46:05Z joe $
# searches for SuperLU headers and libs
# _slu_lib_path(SUPERLU_ROOT, HEADER)
#
# Try to find the subpath unter SUPERLU_ROOT containing HEADER. Try
# SUPERLU_ROOT/"include/superlu", SUPERLU_ROOT/"include", and
# SUPERLU_ROOT/"SRC", in that order. Set the subpath for the library to
# "lib". If HEADER was found in SUPERLU_ROOT/"SRC", check whether
# SUPERLU_ROOT/"lib" is a directory, and set the subpath for the library to
# the empty string "" if it isn't.
#
# Shell variables:
# my_include_path
# The subpath HEADER was found in: "include/superlu", "include", or
# "SRC". Contents is only meaningful for my_slu_found=yes.
# my_lib_path
# The subpath for the library: "lib" or "". Contents is only meaningful
# for my_slu_found=yes.
# my_slu_found
# Whether HEADER was found at all. Either "yes" or "no".
AC_DEFUN([_slu_lib_path],
[
my_include_path=include/superlu
my_lib_path=lib
my_slu_found=yes
if test ! -f "$1/$my_include_path/$2" ; then
#Try to find headers under superlu
my_include_path=include
if test ! -f "$1/$my_include_path/$2" ; then
my_include_path=SRC
if test ! -f "$1/$my_include_path/$2"; then
my_slu_found=no
else
if ! test -d "$1/$my_lib_path"; then
my_lib_path=""
fi
fi
fi
fi
]
)
# _slu_search_versions(SUPERLU_ROOT)
#
# Search for either "slu_ddefs.h" or "dsp_defs.h" using _slu_lib_path().
#
# Shell variables:
# my_slu_header
# The name of the header that was found: first of "slu_ddefs.h" or
# "dsp_defs.h". Contents is only meaningful for my_slu_found=yes.
# my_include_path
# The subpath the header was found in: "include/superlu", "include", or
# "SRC". Contents is only meaningful for my_slu_found=yes.
# my_lib_path
# The subpath for the library: "lib" or "". Contents is only meaningful
# for my_slu_found=yes.
# my_slu_found
# Whether any of the headers. Either "yes" or "no".
AC_DEFUN([_slu_search_versions],
[
my_slu_header=slu_ddefs.h
_slu_lib_path($1, $my_slu_header)
if test "$my_slu_found" != "yes"; then
my_slu_header="dsp_defs.h"
_slu_lib_path($1, $my_slu_header)
fi
]
)
# _slu_search_default()
#
# Search for SuperLU in the default locations "/usr" and "/usr/local".
#
# Shell variables:
# with_superlu
# Root of the SuperLU installation: first of "/usr" and "/usr/local".
# Contents is only meaningful for my_slu_found=yes.
# For other output variables see documentation of _slu_search_versions().
AC_DEFUN([_slu_search_default],
[
with_superlu=/usr
_slu_search_versions($with_superlu)
if test "$my_slu_found" = "no"; then
with_superlu=/usr/local
_slu_search_versions($with_superlu)
fi
]
)
# OPM_PATH_SUPERLU()
#
# REQUIRES: AC_PROG_CC, AX_BLAS
#
# Shell variables:
# with_superlu
# "no", "yes (version 4.3 or newer)", "yes (version 4.2 or older, post 2005)" or "yes (pre 2005)"
# direct_SUPERLU_CPPFLAGS
# direct_SUPERLU_LIBS
# CPPFLAGS and LIBS necessary to link against SuperLU. This variable
# contains no indirect references and is suitable for use inside
# configure. Guaranteed empty if SuperLU was not found.
# SUPERLU_CPPFLAGS
# SUPERLU_LIBS
# CPPFLAGS and LIBS necessary to link against SuperLU. This variable may
# contain indirect references and is suitable for use inside makefiles.
# Guaranteed empty if SuperLU was not found.
# HAVE_SUPERLU
# "0" or "1" depending on whether SuperLU was found.
#
# Substitutions:
# SUPERLU_LIBS
# SUPERLU_CPPFLAGS
# Substitutes the values of the corresponding shell variables.
# ALL_PKG_LIBS
# ALL_PKG_CPPFLAGS
# Adds references to SuperLU's substitutions.
#
# Defines:
# HAVE_SUPERLU
# ENABLE_SUPERLU or undefined. Whether SuperLU was found. The correct
# way to check this is "#if HAVE_SUPERLU": This way SuperLU features will
# be disabled unless ${SUPERLU_CPPFLAGS} was given when compiling.
# SUPERLU_POST_2005_VERSION
# 1 or undefined. A post-2005 version of SuperLU uses the header
# "slu_ddefs.h" while earlier versions use "dsp_defs.h".
# SUPERLU_MIN_VERSION_4_3
# 1 or undefined. SuperLU version 4.3 or newer uses the symbol
# "SLU_DOUBLE" while earlier versions use "DOUBLE".
# HAVE_MEM_USAGE_T_EXPANSIONS
# 1 or undefined. Whether "mem_usage_t.expansions" was found in
# "slu_ddefs.h" or "dsp_defs.h" as apropriate.
#
# Conditionals:
# SUPERLU
AC_DEFUN([OPM_PATH_SUPERLU],[
AC_REQUIRE([AC_PROG_CC])
# we need this for FLIBS
AC_REQUIRE([AC_F77_LIBRARY_LDFLAGS])
AC_REQUIRE([AX_BLAS])
#
# User hints ...
#
my_lib_path=""
my_include_path=""
AC_ARG_WITH([superlu],
[AC_HELP_STRING([--with-superlu],[user defined path to SuperLU library])],
[dnl
if test x"$withval" != xno ; then
# get absolute path
with_superlu=`eval cd $withval > /dev/null 2>&1 && pwd`
if test x"$withval" = xyes; then
# Search in default locations
_slu_search_default
else
# Search for the headers in the specified location
_slu_search_versions(["$with_superlu"])
fi
fi
], [dnl
# Search in default locations
_slu_search_default
])
AC_ARG_WITH([superlu-lib],
[AC_HELP_STRING([--with-superlu-lib],
[The name of the static SuperLU library to link to. By default
the static library with the name superlu.a is tried, but
only if shared linking has failed first. Giving this
options forces static linking for SuperLU.])],
[
if test x"$withval" = xno ; then
with_superlu_lib=
fi
])
AC_ARG_WITH([superlu-blaslib],
[AC_HELP_STRING([--with-superlu-blaslib],
[The name of the static blas library to link to. By default
we try to link to the systems blas library (see
--with-blas). Giving this options forces static linking
for SuperLU.])],
[
if test "$withval" = no ; then
with_superlu_blaslib=
fi
])
# store old values
ac_save_LDFLAGS="$LDFLAGS"
ac_save_CPPFLAGS="$CPPFLAGS"
ac_save_LIBS="$LIBS"
# do nothing if --without-superlu is used
if test x"$with_superlu" != x"no" ; then
# defaultpath
SUPERLU_LIB_PATH="$with_superlu/$my_lib_path"
SUPERLU_INCLUDE_PATH="$with_superlu/$my_include_path"
# set variables so that tests can use them
direct_SUPERLU_CPPFLAGS="-I$SUPERLU_INCLUDE_PATH -DENABLE_SUPERLU"
SUPERLU_CPPFLAGS="-I$SUPERLU_INCLUDE_PATH -DENABLE_SUPERLU"
CPPFLAGS="$CPPFLAGS $direct_SUPERLU_CPPFLAGS"
# check for central header
AC_CHECK_HEADER([$my_slu_header],
[HAVE_SUPERLU="1"],
[
HAVE_SUPERLU="0"
AC_MSG_WARN([$my_slu_header not found in $SUPERLU_INCLUDE_PATH with $CPPFLAGS])
])
# if header is found check for the libs
if test x$HAVE_SUPERLU = x1 ; then
HAVE_SUPERLU=0
# if neither --with-superlu-lib nor --with-superlu-blaslib was
# given, try to link dynamically or with properly names static libs
if test x"$with_superlu_lib$with_superlu_blaslib" = x; then
LDFLAGS="$ac_save_LDFLAGS -L$SUPERLU_LIB_PATH"
LIBS="$ac_save_LIBS"
AC_CHECK_LIB([superlu], [dgssvx],
[
direct_SUPERLU_LIBS="-L$SUPERLU_LIB_PATH -lsuperlu $BLAS_LIBS $FLIBS"
SUPERLU_LIBS="-L$SUPERLU_LIB_PATH -lsuperlu \${BLAS_LIBS} \${FLIBS}"
HAVE_SUPERLU="1"
], [], [$BLAS_LIBS $FLIBS])
fi
if test $HAVE_SUPERLU = 0 && test x"$with_superlu_lib" = x; then
# set the default
with_superlu_lib=superlu.a
fi
if test $HAVE_SUPERLU = 0 &&
test x"$with_superlu_blaslib" = x; then
# try system blas
LDFLAGS="$ac_save_LDFLAGS"
LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $BLAS_LIBS $FLIBS $ac_save_LIBS"
AC_CHECK_FUNC([dgssvx],
[
direct_SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $BLAS_LIBS $FLIBS"
SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib \${BLAS_LIBS} \${FLIBS}"
HAVE_SUPERLU="1"
])
fi
# No default for with_superlu_blaslib
if test $HAVE_SUPERLU = 0 &&
test x"$with_superlu_blaslib" != x; then
# try internal blas
LDFLAGS="$ac_save_LDFLAGS"
LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $SUPERLU_LIB_PATH/$with_superlu_blaslib $FLIBS $ac_save_LIBS"
AC_CHECK_FUNC([dgssvx],
[
direct_SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $SUPERLU_LIB_PATH/$with_superlu_blaslib $FLIBS"
SUPERLU_LIBS="$SUPERLU_LIB_PATH/$with_superlu_lib $SUPERLU_LIB_PATH/$with_superlu_blaslib \${FLIBS}"
HAVE_SUPERLU="1"
])
fi
# test whether SuperLU version is at least 4.3
if test $HAVE_SUPERLU = "1" ; then
AC_CHECK_DECL([SLU_DOUBLE], [SUPERLU_MIN_VERSION_4_3="1"], [], [#include <$my_slu_header>])
fi
fi
else # $with_superlu = no
HAVE_SUPERLU=0
fi
# Inform the user whether SuperLU was sucessfully found
AC_MSG_CHECKING([SuperLU])
if test x$HAVE_SUPERLU = x1 ; then
if test "$my_slu_header" = "slu_ddefs.h"; then
if test x$SUPERLU_MIN_VERSION_4_3 = x1 ; then
with_superlu="yes (version 4.3 or newer)"
else
with_superlu="yes (version 4.2 or older, post 2005)"
fi
else
with_superlu="yes (pre 2005)"
fi
else
with_superlu="no"
fi
AC_MSG_RESULT([$with_superlu])
# check for optional member
if test $HAVE_SUPERLU = 1 ; then
if test "$my_slu_header" = "slu_ddefs.h"; then
AC_CHECK_MEMBERS([mem_usage_t.expansions],[],[],[#include "slu_ddefs.h"])
else
AC_CHECK_MEMBERS([mem_usage_t.expansions],[],[],[#include "dsp_defs.h"])
fi
fi
# substitute variables
if test x$HAVE_SUPERLU = x0 ; then
SUPERLU_LIBS=
SUPERLU_CPPFLAGS=
fi
AC_SUBST([SUPERLU_LIBS])
AC_SUBST([SUPERLU_CPPFLAGS])
#DUNE_ADD_ALL_PKG([SUPERLU], [\${SUPERLU_CPPFLAGS}], [], [\${SUPERLU_LIBS}])
# tell automake
AM_CONDITIONAL(SUPERLU, test x$HAVE_SUPERLU = x1)
# tell the preprocessor
if test x$HAVE_SUPERLU = x1 ; then
AC_DEFINE([HAVE_SUPERLU], [ENABLE_SUPERLU], [Define to ENABLE_SUPERLU if SUPERLU is found])
if test "$my_slu_header" = "slu_ddefs.h"; then
AC_DEFINE([SUPERLU_POST_2005_VERSION], 1, [define to 1 if there is a header slu_ddefs.h in SuperLU])
if test x$SUPERLU_MIN_VERSION_4_3 = x1 ; then
AC_DEFINE([SUPERLU_MIN_VERSION_4_3], 1, [define to 1 if there SLU_DOUBLE imported by header slu_ddefs.h from SuperLU])
fi
fi
fi
# summary
#DUNE_ADD_SUMMARY_ENTRY([SuperLU],[$with_superlu])
# restore variables
LDFLAGS="$ac_save_LDFLAGS"
CPPFLAGS="$ac_save_CPPFLAGS"
LIBS="$ac_save_LIBS"
])

16
opm/core/eclipse/EclipseGridParserHelpers.hpp
View File

@@ -398,10 +398,18 @@ namespace
yv.push_back(table[k][i]);
}
}
// Interpolate
for (int i=0; i<int(indx.size()); ++i) {
table[k][indx[i]] = linearInterpolationExtrap(xv, yv, x[i]);
}
if (xv.empty()) {
// Nothing specified, the entire column is defaulted.
// We insert zeros.
for (int i=0; i<int(indx.size()); ++i) {
table[k][indx[i]] = 0.0;
}
} else {
// Interpolate
for (int i=0; i<int(indx.size()); ++i) {
table[k][indx[i]] = linearInterpolationExtrap(xv, yv, x[i]);
}
}
}
}
}

2
opm/core/eclipse/SpecialEclipseFields.hpp
View File

@@ -1047,7 +1047,7 @@ struct GCONINJE : public SpecialBase
gconinje_line.injector_type_ = readString(is);
gconinje_line.control_mode_ = readString(is);
std::vector<double> double_data(10, -1.0E20);
const int num_to_read = 10;
const int num_to_read = 4;
int num_read = readDefaultedVectorData(is, double_data, num_to_read);
gconinje_line.surface_flow_max_rate_ = double_data[0];
gconinje_line.resv_flow_max_rate_ = double_data[1];

2
opm/core/fluid/IncompPropertiesInterface.hpp
View File

@@ -62,7 +62,7 @@ namespace Opm
/// \return Array of P viscosity values.
virtual const double* viscosity() const = 0;
/// Densities of fluid phases at surface conditions.
/// Densities of fluid phases at reservoir conditions.
/// \return Array of P density values.
virtual const double* density() const = 0;

398
opm/core/grid.c
View File

@@ -18,7 +18,12 @@
*/
#include <opm/core/grid.h>
#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
void
@@ -135,3 +140,396 @@ allocate_grid(size_t ndims ,
return G;
}
#define GRID_NMETA 6
#define GRID_NDIMS 0
#define GRID_NCELLS 1
#define GRID_NFACES 2
#define GRID_NNODES 3
#define GRID_NFACENODES 4
#define GRID_NCELLFACES 5
static void
input_error(FILE *fp, const char * const err)
{
int save_errno = errno;
if (ferror(fp)) {
fprintf(stderr, "%s: %s\n", err, strerror(save_errno));
clearerr(fp);
}
else if (feof(fp)) {
fprintf(stderr, "%s: End-of-file\n", err);
}
errno = save_errno;
}
static struct UnstructuredGrid *
allocate_grid_from_file(FILE *fp, int *has_tag, int *has_indexmap)
{
struct UnstructuredGrid *G;
int save_errno;
unsigned long tmp;
size_t dimens[GRID_NMETA], i;
save_errno = errno;
i = 0;
while ((i < GRID_NMETA) && (fscanf(fp, " %lu", &tmp) == 1)) {
dimens[i] = tmp;
i += 1;
}
if (i == GRID_NMETA) {
if (fscanf(fp, "%d %d", has_tag, has_indexmap) == 2) {
G = allocate_grid(dimens[GRID_NDIMS] ,
dimens[GRID_NCELLS] ,
dimens[GRID_NFACES] ,
dimens[GRID_NFACENODES],
dimens[GRID_NCELLFACES],
dimens[GRID_NNODES] );
if (G != NULL) {
if (! *has_tag) {
free(G->cell_facetag);
G->cell_facetag = NULL;
}
if (*has_indexmap) {
G->global_cell =
malloc(dimens[GRID_NCELLS] * sizeof *G->global_cell);
/* Allocation failure checked elsewhere. */
}
G->number_of_cells = (int) dimens[GRID_NCELLS];
G->number_of_faces = (int) dimens[GRID_NFACES];
G->number_of_nodes = (int) dimens[GRID_NNODES];
G->dimensions = (int) dimens[GRID_NDIMS];
i = 0;
while ((i < dimens[GRID_NDIMS]) &&
(fscanf(fp, "%d", & G->cartdims[ i ]) == 1)) {
i += 1;
}
if (i < dimens[GRID_NDIMS]) {
input_error(fp, "Unable to read Cartesian dimensions");
destroy_grid(G);
G = NULL;
}
else {
/* Account for dimens[GRID_DIMS] < 3 */
size_t n = (sizeof G->cartdims) / (sizeof G->cartdims[0]);
for (; i < n; i++) { G->cartdims[ i ] = 1; }
}
}
}
else {
input_error(fp, "Unable to read grid predicates");
G = NULL;
}
}
else {
input_error(fp, "Unable to read grid dimensions");
G = NULL;
}
errno = save_errno;
return G;
}
static int
read_grid_nodes(FILE *fp, struct UnstructuredGrid *G)
{
int save_errno;
size_t i, n;
save_errno = errno;
n = G->dimensions;
n *= G->number_of_nodes;
i = 0;
while ((i < n) &&
(fscanf(fp, " %lf", & G->node_coordinates[ i ]) == 1)) {
i += 1;
}
if (i < n) {
input_error(fp, "Unable to read node coordinates");
}
errno = save_errno;
return i == n;
}
static int
read_grid_faces(FILE *fp, struct UnstructuredGrid *G)
{
int save_errno, ok;
size_t nf, nfn, i;
save_errno = errno;
nf = G->number_of_faces;
/* G->face_nodepos */
i = 0;
while ((i < nf + 1) &&
(fscanf(fp, " %d", & G->face_nodepos[ i ]) == 1)) {
i += 1;
}
ok = i == nf + 1;
if (! ok) {
input_error(fp, "Unable to read node indirection array");
}
else {
/* G->face_nodes */
nfn = G->face_nodepos[ nf ];
i = 0;
while ((i < nfn) && (fscanf(fp, " %d", & G->face_nodes[ i ]) == 1)) {
i += 1;
}
ok = i == nfn;
if (! ok) {
input_error(fp, "Unable to read face-nodes");
}
}
if (ok) {
/* G->face_cells */
i = 0;
while ((i < 2 * nf) && (fscanf(fp, " %d", & G->face_cells[ i ]) == 1)) {
i += 1;
}
ok = i == 2 * nf;
if (! ok) {
input_error(fp, "Unable to read neighbourship");
}
}
if (ok) {
/* G->face_areas */
i = 0;
while ((i < nf) && (fscanf(fp, " %lf", & G->face_areas[ i ]) == 1)) {
i += 1;
}
ok = i == nf;
if (! ok) {
input_error(fp, "Unable to read face areas");
}
}
if (ok) {
/* G->face_centroids */
size_t n;
n = G->dimensions;
n *= nf;
i = 0;
while ((i < n) && (fscanf(fp, " %lf", & G->face_centroids[ i ]) == 1)) {
i += 1;
}
ok = i == n;
if (! ok) {
input_error(fp, "Unable to read face centroids");
}
}
if (ok) {
/* G->face_normals */
size_t n;
n = G->dimensions;
n *= nf;
i = 0;
while ((i < n) && (fscanf(fp, " %lf", & G->face_normals[ i ]) == 1)) {
i += 1;
}
ok = i == n;
if (! ok) {
input_error(fp, "Unable to read face normals");
}
}
errno = save_errno;
return ok;
}
static int
read_grid_cells(FILE *fp, int has_tag, int has_indexmap,
struct UnstructuredGrid *G)
{
int save_errno, ok;
size_t nc, ncf, i;
save_errno = errno;
nc = G->number_of_cells;
/* G->cell_facepos */
i = 0;
while ((i < nc + 1) && (fscanf(fp, " %d", & G->cell_facepos[ i ]) == 1)) {
i += 1;
}
ok = i == nc + 1;
if (! ok) {
input_error(fp, "Unable to read face indirection array");
}
else {
/* G->cell_faces (and G->cell_facetag if applicable) */
ncf = G->cell_facepos[ nc ];
i = 0;
if (has_tag) {
assert (G->cell_facetag != NULL);
while ((i < ncf) &&
(fscanf(fp, " %d %d",
& G->cell_faces [ i ],
& G->cell_facetag[ i ]) == 2)) {
i += 1;
}
}
else {
while ((i < ncf) &&
(fscanf(fp, " %d", & G->cell_faces[ i ]) == 1)) {
i += 1;
}
}
ok = i == ncf;
if (! ok) {
input_error(fp, "Unable to read cell-faces");
}
}
if (ok) {
/* G->global_cell if applicable */
if (has_indexmap) {
i = 0;
if (G->global_cell != NULL) {
while ((i < nc) &&
(fscanf(fp, " %d", & G->global_cell[ i ]) == 1)) {
i += 1;
}
}
else {
int discard;
while ((i < nc) && (fscanf(fp, " %d", & discard) == 1)) {
i += 1;
}
}
}
else {
assert (G->global_cell == NULL);
i = nc;
}
ok = i == nc;
if (! ok) {
input_error(fp, "Unable to read global cellmap");
}
}
if (ok) {
/* G->cell_volumes */
i = 0;
while ((i < nc) && (fscanf(fp, " %lf", & G->cell_volumes[ i ]) == 1)) {
i += 1;
}
ok = i == nc;
if (! ok) {
input_error(fp, "Unable to read cell volumes");
}
}
if (ok) {
/* G->cell_centroids */
size_t n;
n = G->dimensions;
n *= nc;
i = 0;
while ((i < n) && (fscanf(fp, " %lf", & G->cell_centroids[ i ]) == 1)) {
i += 1;
}
ok = i == n;
if (! ok) {
input_error(fp, "Unable to read cell centroids");
}
}
errno = save_errno;
return ok;
}
struct UnstructuredGrid *
read_grid(const char *fname)
{
struct UnstructuredGrid *G;
FILE *fp;
int save_errno;
int has_tag, has_indexmap, ok;
save_errno = errno;
fp = fopen(fname, "rt");
if (fp != NULL) {
G = allocate_grid_from_file(fp, & has_tag, & has_indexmap);
ok = G != NULL;
if (ok) { ok = read_grid_nodes(fp, G); }
if (ok) { ok = read_grid_faces(fp, G); }
if (ok) { ok = read_grid_cells(fp, has_tag, has_indexmap, G); }
if (! ok) {
destroy_grid(G);
G = NULL;
}
fclose(fp);
}
else {
G = NULL;
}
errno = save_errno;
return G;
}

11
opm/core/grid.h
View File

@@ -273,6 +273,17 @@ allocate_grid(size_t ndims ,
size_t ncellfaces,
size_t nnodes );
/**
* Import a grid from a character representation stored in file.
*
* @param[in] fname File name.
* @return Fully formed UnstructuredGrid with all fields allocated and filled.
* Returns @c NULL in case of allocation failure.
*/
struct UnstructuredGrid *
read_grid(const char *fname);
#ifdef __cplusplus
}
#endif

4
opm/core/grid/cpgpreprocess/preprocess.h
View File

@@ -56,8 +56,8 @@ extern "C" {
struct grdecl {
int dims[3]; /**< Cartesian box dimensions. */
const double *coord; /**< Pillar end-points. */
const double *zcorn; /**< Explicit "active" map. May be NULL.*/
const int *actnum; /**< Corner-point depths. */
const double *zcorn; /**< Corner-point depths. */
const int *actnum; /**< Explicit "active" map. May be NULL.*/
const double *mapaxes; /**< 6 Element rotation vector - can be NULL. */
};

86
opm/core/newwells.c
View File

@@ -468,6 +468,8 @@ add_well(enum WellType type ,
struct WellMgmt *m;
assert (W != NULL);
nw = W->number_of_wells;
nperf_tot = W->well_connpos[nw];
@@ -532,6 +534,7 @@ append_well_controls(enum WellControlType type,
struct WellControlMgmt *m;
assert (W != NULL);
assert ((0 <= well_index) && (well_index < W->number_of_wells));
ctrl = W->ctrls[well_index];
np = W->number_of_phases;
@@ -567,8 +570,13 @@ void
set_current_control(int well_index, int current_control, struct Wells *W)
/* ---------------------------------------------------------------------- */
{
assert (W != NULL);
assert ((0 <= well_index) && (well_index < W->number_of_wells));
assert (W->ctrls[well_index] != NULL);
assert (current_control < W->ctrls[well_index]->num);
W->ctrls[well_index]->current = current_control;
}
@@ -578,7 +586,85 @@ void
clear_well_controls(int well_index, struct Wells *W)
/* ---------------------------------------------------------------------- */
{
assert (W != NULL);
assert ((0 <= well_index) && (well_index < W->number_of_wells));
if (W->ctrls[well_index] != NULL) {
W->ctrls[well_index]->num = 0;
}
}
/* ---------------------------------------------------------------------- */
struct Wells *
clone_wells(const struct Wells *W)
/* ---------------------------------------------------------------------- */
{
int c, np, nperf, ok, pos, w;
double target;
const int *cells;
const double *WI, *comp_frac, *distr;
enum WellControlType type;
const struct WellControls *ctrls;
struct Wells *ret;
if (W == NULL) {
ret = NULL;
}
else {
np = W->number_of_phases;
ret = create_wells(W->number_of_phases, W->number_of_wells,
W->well_connpos[ W->number_of_wells ]);
if (ret != NULL) {
pos = W->well_connpos[ 0 ];
ok = 1;
for (w = 0; ok && (w < W->number_of_wells); w++) {
nperf = W->well_connpos[w + 1] - pos;
cells = W->well_cells + pos;
WI = W->WI != NULL ? W->WI + pos : NULL;
comp_frac = W->comp_frac != NULL ? W->comp_frac + w*np : NULL;
ok = add_well(W->type[ w ], W->depth_ref[ w ], nperf,
comp_frac, cells, WI, W->name[ w ], ret);
/* Capacity should be sufficient from create_wells() */
assert (ok);
ctrls = W->ctrls[ w ];
if (ok) {
ok = well_controls_reserve(ctrls->num, np, ret->ctrls[ w ]);
for (c = 0; ok && (c < ctrls->num); c++) {
type = ctrls->type [ c ];
target = ctrls->target[ c ];
distr = & ctrls->distr [ c * np ];
ok = append_well_controls(type, target, distr, w, ret);
/* Capacity *should* be sufficient from
* well_controls_reserve() */
assert (ok);
}
}
if (ok) {
set_current_control(w, ctrls->current, ret);
}
pos = W->well_connpos[w + 1];
}
if (! ok) {
destroy_wells(ret);
ret = NULL;
}
}
}
return ret;
}

13
opm/core/newwells.h
View File

@@ -235,6 +235,19 @@ void
destroy_wells(struct Wells *W);
/**
* Create a deep-copy (i.e., clone) of an existing Wells object, including its
* controls.
*
* @param[in] W Existing Wells object.
* @return Complete clone of the input object. Dispose of resources using
* function destroy_wells() when no longer needed. Returns @c NULL in case of
* allocation failure.
*/
struct Wells *
clone_wells(const struct Wells *W);
#ifdef __cplusplus
}
#endif

17
opm/core/pressure/CompressibleTpfa.cpp
View File

@@ -632,6 +632,23 @@ namespace Opm
&state.pressure()[0],
&state.faceflux()[0],
&state.facepressure()[0]);
// Compute well perforation pressures (not done by the C code).
if (wells_ != 0) {
const int nw = wells_->number_of_wells;
const int np = props_.numPhases();
for (int w = 0; w < nw; ++w) {
const double* comp_frac = &wells_->comp_frac[np*w];
for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
const double bhp = well_state.bhp()[w];
double perf_p = bhp;
for (int phase = 0; phase < np; ++phase) {
perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
}
well_state.perfPress()[j] = perf_p;
}
}
}
}
} // namespace Opm

77
opm/core/simulator/SimulatorCompressibleTwophase.cpp
View File

@@ -196,6 +196,7 @@ namespace Opm
if (!file) {
THROW("Failed to open " << fname.str());
}
file.precision(15);
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
@@ -319,15 +320,11 @@ namespace Opm
Opm::time::StopWatch step_timer;
Opm::time::StopWatch total_timer;
total_timer.start();
double init_satvol[2] = { 0.0 };
double satvol[2] = { 0.0 };
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double init_surfvol[2] = { 0.0 };
double inplace_surfvol[2] = { 0.0 };
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init
<< " " << init_satvol[1]/tot_porevol_init << std::endl;
Opm::computeSaturatedVol(porevol, state.surfacevol(), init_surfvol);
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport;
@@ -410,8 +407,8 @@ namespace Opm
// Optionally, check if well controls are satisfied.
if (check_well_controls_) {
Opm::computePhaseFlowRatesPerWell(*wells_,
fractional_flows,
well_state.perfRates(),
fractional_flows,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
@@ -434,9 +431,8 @@ namespace Opm
computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol);
}
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
wells_, well_state.perfRates(), transport_src);
// Process transport sources from well flows.
Opm::computeTransportSource(props_, wells_, well_state, transport_src);
// Solve transport.
transport_timer.start();
@@ -445,13 +441,20 @@ namespace Opm
stepsize /= double(num_transport_substeps_);
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
}
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &state.pressure()[0],
&initial_porevol[0], &porevol[0], &transport_src[0], stepsize,
state.saturation(), state.surfacevol());
Opm::computeInjectedProduced(props_,
state.pressure(), state.surfacevol(), state.saturation(),
transport_src, stepsize, injected, produced);
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
Opm::computeInjectedProduced(props_, state, transport_src, stepsize,
substep_injected, substep_produced);
injected[0] += substep_injected[0];
injected[1] += substep_injected[1];
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
if (gravity_ != 0 && use_segregation_split_) {
tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol());
}
@@ -462,35 +465,35 @@ namespace Opm
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
std::cout.precision(5);
const int width = 18;
std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n";
std::cout << " Saturated volumes: "
<< std::setw(width) << satvol[0]/tot_porevol_init
<< std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
std::cout << " Injected volumes: "
<< std::setw(width) << injected[0]/tot_porevol_init
<< std::setw(width) << injected[1]/tot_porevol_init << std::endl;
std::cout << " Produced volumes: "
<< std::setw(width) << produced[0]/tot_porevol_init
<< std::setw(width) << produced[1]/tot_porevol_init << std::endl;
std::cout << " Total inj volumes: "
<< std::setw(width) << tot_injected[0]/tot_porevol_init
<< std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
std::cout << " Total prod volumes: "
<< std::setw(width) << tot_produced[0]/tot_porevol_init
<< std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
std::cout << " In-place + prod - inj: "
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
std::cout << " Init - now - pr + inj: "
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
std::cout << "\nMass balance report.\n";
std::cout << " Injected surface volumes: "
<< std::setw(width) << injected[0]
<< std::setw(width) << injected[1] << std::endl;
std::cout << " Produced surface volumes: "
<< std::setw(width) << produced[0]
<< std::setw(width) << produced[1] << std::endl;
std::cout << " Total inj surface volumes: "
<< std::setw(width) << tot_injected[0]
<< std::setw(width) << tot_injected[1] << std::endl;
std::cout << " Total prod surface volumes: "
<< std::setw(width) << tot_produced[0]
<< std::setw(width) << tot_produced[1] << std::endl;
const double balance[2] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1] };
std::cout << " Initial - inplace + inj - prod: "
<< std::setw(width) << balance[0]
<< std::setw(width) << balance[1]
<< std::endl;
std::cout << " Relative mass error: "
<< std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
<< std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
<< std::endl;
std::cout.precision(8);

16
opm/core/simulator/SimulatorIncompTwophase.cpp
View File

@@ -244,6 +244,7 @@ namespace Opm
if (!file) {
THROW("Failed to open " << fname.str());
}
file.precision(15);
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
@@ -398,8 +399,6 @@ namespace Opm
total_timer.start();
double init_satvol[2] = { 0.0 };
double satvol[2] = { 0.0 };
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
@@ -486,8 +485,8 @@ namespace Opm
// Optionally, check if well controls are satisfied.
if (check_well_controls_) {
Opm::computePhaseFlowRatesPerWell(*wells_,
fractional_flows,
well_state.perfRates(),
fractional_flows,
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
@@ -521,10 +520,19 @@ namespace Opm
stepsize /= double(num_transport_substeps_);
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
}
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0],
stepsize, state.saturation());
Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize, injected, produced);
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize,
substep_injected, substep_produced);
injected[0] += substep_injected[0];
injected[1] += substep_injected[1];
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
if (use_segregation_split_) {
tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation());
}

8
opm/core/simulator/WellState.hpp
View File

@@ -49,7 +49,8 @@ namespace Opm
bhp_[w] = state.pressure()[cell];
}
}
perfrates_.resize(wells->well_connpos[nw]);
perfrates_.resize(wells->well_connpos[nw], 0.0);
perfpress_.resize(wells->well_connpos[nw], -1e100);
}
}
@@ -61,9 +62,14 @@ namespace Opm
std::vector<double>& perfRates() { return perfrates_; }
const std::vector<double>& perfRates() const { return perfrates_; }
/// One pressure per well connection.
std::vector<double>& perfPress() { return perfpress_; }
const std::vector<double>& perfPress() const { return perfpress_; }
private:
std::vector<double> bhp_;
std::vector<double> perfrates_;
std::vector<double> perfpress_;
};
} // namespace Opm

6
opm/core/transport/GravityColumnSolver_impl.hpp
View File

@@ -175,10 +175,10 @@ namespace Opm
}
// model_.sourceTerms(); // Not needed
// Solve.
const int num_rhs = 1;
int info = 0;
const MAT_SIZE_T num_rhs = 1, colSize = col_size;
MAT_SIZE_T info = 0;
// Solution will be written to rhs.
dgtsv_(&col_size, &num_rhs, DL, D, DU, &rhs[0], &col_size, &info);
dgtsv_(&colSize, &num_rhs, DL, D, DU, &rhs[0], &colSize, &info);
if (info != 0) {
THROW("Lapack reported error in dgtsv: " << info);
}

2
opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp
View File

@@ -152,7 +152,7 @@ namespace Opm
B_cell = 1.0/tm.A_[np*np*cell + 0];
double src_flux = -tm.source_[cell];
bool src_is_inflow = src_flux < 0.0;
influx = src_is_inflow ? src_flux : 0.0;
influx = src_is_inflow ? B_cell* src_flux : 0.0;
outflux = !src_is_inflow ? src_flux : 0.0;
comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell];
dtpv = tm.dt_/tm.porevolume0_[cell];

7
opm/core/transport/reorder/TransportModelInterface.cpp
View File

@@ -20,6 +20,7 @@
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/grid.h>
#include <opm/core/utility/StopWatch.hpp>
#include <vector>
#include <cassert>
@@ -31,10 +32,15 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g
std::vector<int> sequence(grid.number_of_cells);
std::vector<int> components(grid.number_of_cells + 1);
int ncomponents;
time::StopWatch clock;
clock.start();
compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
clock.stop();
std::cout << "Topological sort took: " << clock.secsSinceStart() << " seconds." << std::endl;
// Invoke appropriate solve method for each interdependent component.
for (int comp = 0; comp < ncomponents; ++comp) {
#if 0
#ifdef MATLAB_MEX_FILE
// \TODO replace this with general signal handling code, check if it costs performance.
if (interrupt_signal) {
@@ -42,6 +48,7 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g
"cells finished.\n", i, grid.number_of_cells);
break;
}
#endif
#endif
const int comp_size = components[comp + 1] - components[comp];
if (comp_size == 1) {

122
opm/core/transport/reorder/TransportModelTracerTof.cpp Normal file
View File

@@ -0,0 +1,122 @@
/*
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/>.
*/
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <numeric>
#include <cmath>
namespace Opm
{
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid)
: grid_(grid)
{
}
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void TransportModelTracerTof::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
}
#endif
tof.resize(grid_.number_of_cells);
std::fill(tof.begin(), tof.end(), 0.0);
tof_ = &tof[0];
reorderAndTransport(grid_, darcyflux);
}
void TransportModelTracerTof::solveSingleCell(const int cell)
{
// Compute flux terms.
// Sources have zero tof, and therefore do not contribute
// to upwind_term. Sinks on the other hand, must be added
// to the downwind_flux (note sign change resulting from
// different sign conventions: pos. source is injection,
// pos. flux is outflow).
double upwind_term = 0.0;
double downwind_flux = std::max(-source_[cell], 0.0);
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
int f = grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == grid_.face_cells[2*f]) {
flux = darcyflux_[f];
other = grid_.face_cells[2*f+1];
} else {
flux =-darcyflux_[f];
other = grid_.face_cells[2*f];
}
// Add flux to upwind_term or downwind_flux, if interior.
if (other != -1) {
if (flux < 0.0) {
upwind_term += flux*tof_[other];
} else {
downwind_flux += flux;
}
}
}
// Compute tof.
tof_[cell] = (porevolume_[cell] - upwind_term)/downwind_flux;
}
void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
}
}
} // namespace Opm

74
opm/core/transport/reorder/TransportModelTracerTof.hpp Normal file
View File

@@ -0,0 +1,74 @@
/*
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 OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <vector>
#include <map>
#include <ostream>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
/// Implements a first-order finite volume solver for
/// (single-phase) time-of-flight using reordering.
/// The equation solved is:
/// v \cdot \grad\tau = \phi
/// where v is the fluid velocity, \tau is time-of-flight and
/// \phi is the porosity. This is a boundary value problem, where
/// \tau is specified to be zero on all inflow boundaries.
class TransportModelTracerTof : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTof(const UnstructuredGrid& grid);
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
private:
const UnstructuredGrid& grid_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
double* tof_;
};
} // namespace Opm
#endif // OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED

671
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp Normal file
View File

@@ -0,0 +1,671 @@
/*
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/>.
*/
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/linalg/blas_lapack.h>
#include <algorithm>
#include <numeric>
#include <cmath>
namespace Opm
{
// --------------- Helpers for TransportModelTracerTofDiscGal ---------------
/// A class providing discontinuous Galerkin basis functions.
struct DGBasis
{
static int numBasisFunc(const int dimensions,
const int degree)
{
switch (dimensions) {
case 1:
return degree + 1;
case 2:
return (degree + 2)*(degree + 1)/2;
case 3:
return (degree + 3)*(degree + 2)*(degree + 1)/6;
default:
THROW("Dimensions must be 1, 2 or 3.");
}
}
/// Evaluate all nonzero basis functions at x,
/// writing to f_x. The array f_x must have
/// size numBasisFunc(grid.dimensions, degree).
///
/// The basis functions are the following
/// Degree 0: 1.
/// Degree 1: x - xc, y - yc, z - zc etc.
/// Further degrees await development.
static void eval(const UnstructuredGrid& grid,
const int cell,
const int degree,
const double* x,
double* f_x)
{
const int dim = grid.dimensions;
const double* cc = grid.cell_centroids + dim*cell;
// Note intentional fallthrough in this switch statement!
switch (degree) {
case 1:
for (int ix = 0; ix < dim; ++ix) {
f_x[1 + ix] = x[ix] - cc[ix];
}
case 0:
f_x[0] = 1;
break;
default:
THROW("Maximum degree is 1 for now.");
}
}
/// Evaluate gradients of all nonzero basis functions at x,
/// writing to grad_f_x. The array grad_f_x must have size
/// numBasisFunc(grid.dimensions, degree) * grid.dimensions.
/// The <grid.dimensions> components of the first basis function
/// gradient come before the components of the second etc.
static void evalGrad(const UnstructuredGrid& grid,
const int /*cell*/,
const int degree,
const double* /*x*/,
double* grad_f_x)
{
const int dim = grid.dimensions;
const int num_basis = numBasisFunc(dim, degree);
std::fill(grad_f_x, grad_f_x + num_basis*dim, 0.0);
if (degree > 1) {
THROW("Maximum degree is 1 for now.");
} else if (degree == 1) {
for (int ix = 0; ix < dim; ++ix) {
grad_f_x[dim*(ix + 1) + ix] = 1.0;
}
}
}
};
static void cross(const double* a, const double* b, double* res)
{
res[0] = a[1]*b[2] - a[2]*b[1];
res[1] = a[2]*b[0] - a[0]*b[2];
res[2] = a[0]*b[1] - a[1]*b[0];
}
static double triangleArea3d(const double* p0,
const double* p1,
const double* p2)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double cr[3];
cross(a, b, cr);
return 0.5*std::sqrt(cr[0]*cr[0] + cr[1]*cr[1] + cr[2]*cr[2]);
}
/// Calculates the determinant of a 3 x 3 matrix, represented as
/// three three-dimensional arrays.
static double determinantOf(const double* a0,
const double* a1,
const double* a2)
{
return
a0[0] * (a1[1] * a2[2] - a2[1] * a1[2]) -
a0[1] * (a1[0] * a2[2] - a2[0] * a1[2]) +
a0[2] * (a1[0] * a2[1] - a2[0] * a1[1]);
}
/// Computes the volume of a tetrahedron consisting of 4 vertices
/// with 3-dimensional coordinates
static double tetVolume(const double* p0,
const double* p1,
const double* p2,
const double* p3)
{
double a[3] = { p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2] };
double b[3] = { p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2] };
double c[3] = { p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2] };
return std::fabs(determinantOf(a, b, c) / 6.0);
}
/// A class providing numerical quadrature for cells.
/// In general: \int_{cell} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by cell volume,
/// so weights always sum to cell volume.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = cell volume, x_0 = cell centroid
/// Degree 2 method:
/// Based on subdivision of each cell face into triangles
/// with the face centroid as a common vertex, and then
/// subdividing the cell into tetrahedra with the cell
/// centroid as a common vertex. Then apply the tetrahedron
/// rule with the following 4 nodes (uniform weights):
/// a = 0.138196601125010515179541316563436
/// x_i has all barycentric coordinates = a, except for
/// the i'th coordinate which is = 1 - 3a.
/// This rule is from http://nines.cs.kuleuven.be/ecf,
/// it is the second degree 2 4-point rule for tets,
/// referenced to Stroud(1971).
/// The tetrahedra are numbered T_{i,j}, and are given by the
/// cell centroid C, the face centroid FC_i, and two nodes
/// of face i: FN_{i,j}, FN_{i,j+1}.
class CellQuadrature
{
public:
CellQuadrature(const UnstructuredGrid& grid,
const int cell,
const int degree)
: grid_(grid), cell_(cell), degree_(degree)
{
if (degree > 2) {
THROW("CellQuadrature exact for polynomial degrees > 1 not implemented.");
}
if (degree == 2) {
// Prepare subdivision.
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
int sumnodes = 0;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
sumnodes += grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
}
return 4*sumnodes;
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
if (degree_ < 2) {
std::copy(cc, cc + dim, coord);
return;
}
// Degree 2 case.
int tetindex = index / 4;
const int subindex = index % 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
const double a = 0.138196601125010515179541316563436;
// Barycentric coordinates of our point in the tet.
double baryc[4] = { a, a, a, a };
baryc[subindex] = 1.0 - 3.0*a;
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = baryc[0]*cc[dd] + baryc[1]*fc[dd] + baryc[2]*n0c[dd] + baryc[3]*n1c[dd];
}
return;
}
THROW("Should never reach this point.");
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.cell_volumes[cell_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* cc = grid_.cell_centroids + dim*cell_;
int tetindex = index / 4;
const double* nc = grid_.node_coordinates;
for (int hf = grid_.cell_facepos[cell_]; hf < grid_.cell_facepos[cell_ + 1]; ++hf) {
const int face = grid_.cell_faces[hf];
const int nfn = grid_.face_nodepos[face + 1] - grid_.face_nodepos[face];
if (nfn <= tetindex) {
// Our tet is not associated with this face.
tetindex -= nfn;
continue;
}
const double* fc = grid_.face_centroids + dim*face;
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face];
const int node0 = fnodes[tetindex];
const int node1 = fnodes[(tetindex + 1) % nfn];
const double* n0c = nc + dim*node0;
const double* n1c = nc + dim*node1;
return 0.25*tetVolume(cc, fc, n0c, n1c);
}
THROW("Should never reach this point.");
}
private:
const UnstructuredGrid& grid_;
const int cell_;
const int degree_;
};
/// A class providing numerical quadrature for faces.
/// In general: \int_{face} g(x) dx = \sum_{i=0}^{n-1} w_i g(x_i).
/// Note that this class does multiply weights by face area,
/// so weights always sum to face area.
/// Degree 1 method:
/// Midpoint (centroid) method.
/// n = 1, w_0 = face area, x_0 = face centroid
/// Degree 2 method:
/// Based on subdivision of the face into triangles,
/// with the centroid as a common vertex, and the triangle
/// edge midpoint rule.
/// Triangle i consists of the centroid C, nodes N_i and N_{i+1}.
/// Its area is A_i.
/// n = 2 * nn (nn = num nodes in face)
/// For i = 0..(nn-1):
/// w_i = 1/3 A_i.
/// w_{nn+i} = 1/3 A_{i-1} + 1/3 A_i
/// x_i = (N_i + N_{i+1})/2
/// x_{nn+i} = (C + N_i)/2
/// All N and A indices are interpreted cyclic, modulus nn.
class FaceQuadrature
{
public:
FaceQuadrature(const UnstructuredGrid& grid,
const int face,
const int degree)
: grid_(grid), face_(face), degree_(degree)
{
if (grid_.dimensions != 3) {
THROW("FaceQuadrature only implemented for 3D case.");
}
if (degree_ > 2) {
THROW("FaceQuadrature exact for polynomial degrees > 2 not implemented.");
}
}
int numQuadPts() const
{
if (degree_ < 2) {
return 1;
}
// Degree 2 case.
return 2 * (grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_]);
}
void quadPtCoord(const int index, double* coord) const
{
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
if (degree_ < 2) {
std::copy(fc, fc + dim, coord);
return;
}
// Degree 2 case.
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node0 + dd] + nc[dim*node1 + dd]);
}
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node = fnodes[index - nn];
for (int dd = 0; dd < dim; ++dd) {
coord[dd] = 0.5*(nc[dim*node + dd] + fc[dd]);
}
}
}
double quadPtWeight(const int index) const
{
if (degree_ < 2) {
return grid_.face_areas[face_];
}
// Degree 2 case.
const int dim = grid_.dimensions;
const double* fc = grid_.face_centroids + dim*face_;
const int nn = grid_.face_nodepos[face_ + 1] - grid_.face_nodepos[face_];
const int* fnodes = grid_.face_nodes + grid_.face_nodepos[face_];
const double* nc = grid_.node_coordinates;
if (index < nn) {
// Boundary edge midpoint.
const int node0 = fnodes[index];
const int node1 = fnodes[(index + 1)%nn];
const double area = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
return area / 3.0;
} else {
// Interiour edge midpoint.
// Recall that index is now in [nn, 2*nn).
const int node0 = fnodes[(index - 1) % nn];
const int node1 = fnodes[index - nn];
const int node2 = fnodes[(index + 1) % nn];
const double area0 = triangleArea3d(nc + dim*node1, nc + dim*node0, fc);
const double area1 = triangleArea3d(nc + dim*node2, nc + dim*node1, fc);
return (area0 + area1) / 3.0;
}
}
private:
const UnstructuredGrid& grid_;
const int face_;
const int degree_;
};
// Initial version: only a constant interpolation.
static void interpolateVelocity(const UnstructuredGrid& grid,
const int cell,
const double* darcyflux,
const double* /*x*/,
double* v)
{
const int dim = grid.dimensions;
std::fill(v, v + dim, 0.0);
const double* cc = grid.cell_centroids + cell*dim;
for (int hface = grid.cell_facepos[cell]; hface < grid.cell_facepos[cell+1]; ++hface) {
const int face = grid.cell_faces[hface];
const double* fc = grid.face_centroids + face*dim;
double flux = 0.0;
if (cell == grid.face_cells[2*face]) {
flux = darcyflux[face];
} else {
flux = -darcyflux[face];
}
for (int dd = 0; dd < dim; ++dd) {
v[dd] += flux * (fc[dd] - cc[dd]) / grid.cell_volumes[cell];
}
}
}
// --------------- Methods of TransportModelTracerTofDiscGal ---------------
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid)
: grid_(grid),
coord_(grid.dimensions),
velocity_(grid.dimensions)
{
}
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] degree Polynomial degree of DG basis functions used.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
source_ = source;
#ifndef NDEBUG
// Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) {
THROW("Sources do not sum to zero: " << cum_src);
}
#endif
degree_ = degree;
const int num_basis = DGBasis::numBasisFunc(grid_.dimensions, degree_);
tof_coeff.resize(num_basis*grid_.number_of_cells);
std::fill(tof_coeff.begin(), tof_coeff.end(), 0.0);
tof_coeff_ = &tof_coeff[0];
rhs_.resize(num_basis);
jac_.resize(num_basis*num_basis);
basis_.resize(num_basis);
basis_nb_.resize(num_basis);
grad_basis_.resize(num_basis*grid_.dimensions);
reorderAndTransport(grid_, darcyflux);
}
void TransportModelTracerTofDiscGal::solveSingleCell(const int cell)
{
// Residual:
// For each cell K, basis function b_j (spanning V_h),
// writing the solution u_h|K = \sum_i c_i b_i
// Res = - \int_K \sum_i c_i b_i v(x) \cdot \grad b_j dx
// + \int_{\partial K} F(u_h, u_h^{ext}, v(x) \cdot n) b_j ds
// - \int_K \phi b_j
// This is linear in c_i, so we do not need any nonlinear iterations.
// We assemble the jacobian and the right-hand side. The residual is
// equal to Res = Jac*c - rhs, and we compute rhs directly.
const int dim = grid_.dimensions;
const int num_basis = DGBasis::numBasisFunc(dim, degree_);
std::fill(rhs_.begin(), rhs_.end(), 0.0);
std::fill(jac_.begin(), jac_.end(), 0.0);
// Compute cell residual contribution.
// Note: Assumes that \int_K b_j = 0 for all j > 0
rhs_[0] += porevolume_[cell];
// Compute upstream residual contribution.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
int upstream_cell = -1;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
upstream_cell = grid_.face_cells[2*face+1];
} else {
flux = -darcyflux_[face];
upstream_cell = grid_.face_cells[2*face];
}
if (upstream_cell < 0) {
// This is an outer boundary. Assumed tof = 0 on inflow, so no contribution.
continue;
}
if (flux >= 0.0) {
// This is an outflow boundary.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} u_h^{ext} (v(x) \cdot n) b_j ds
// (where u_h^{ext} is the upstream unknown (tof)).
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::eval(grid_, upstream_cell, degree_, &coord_[0], &basis_nb_[0]);
const double tof_upstream = std::inner_product(basis_nb_.begin(), basis_nb_.end(),
tof_coeff_ + num_basis*upstream_cell, 0.0);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
rhs_[j] -= w * tof_upstream * normal_velocity * basis_[j];
}
}
}
// Compute cell jacobian contribution. We use Fortran ordering
// for jac_, i.e. rows cycling fastest.
{
CellQuadrature quad(grid_, cell, 2*degree_ - 1);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// b_i (v \cdot \grad b_j)
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
DGBasis::evalGrad(grid_, cell, degree_, &coord_[0], &grad_basis_[0]);
interpolateVelocity(grid_, cell, darcyflux_, &coord_[0], &velocity_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
for (int dd = 0; dd < dim; ++dd) {
jac_[j*num_basis + i] -= w * basis_[j] * grad_basis_[dim*i + dd] * velocity_[dd];
}
}
}
}
}
// Compute downstream jacobian contribution from faces.
for (int hface = grid_.cell_facepos[cell]; hface < grid_.cell_facepos[cell+1]; ++hface) {
const int face = grid_.cell_faces[hface];
double flux = 0.0;
if (cell == grid_.face_cells[2*face]) {
flux = darcyflux_[face];
} else {
flux = -darcyflux_[face];
}
if (flux <= 0.0) {
// This is an inflow boundary.
continue;
}
// Do quadrature over the face to compute
// \int_{\partial K} b_i (v(x) \cdot n) b_j ds
const double normal_velocity = flux / grid_.face_areas[face];
FaceQuadrature quad(grid_, face, 2*degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
// u^ext flux B (B = {b_j})
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * normal_velocity * basis_[j];
}
}
}
}
// Compute downstream jacobian contribution from sink terms.
// Contribution from inflow sources would be
// similar to the contribution from upstream faces, but
// it is zero since we let all external inflow be associated
// with a zero tof.
if (source_[cell] < 0.0) {
// A sink.
const double flux = -source_[cell]; // Sign convention for flux: outflux > 0.
const double flux_density = flux / grid_.cell_volumes[cell];
// Do quadrature over the cell to compute
// \int_{K} b_i flux b_j dx
CellQuadrature quad(grid_, cell, 2*degree_);
for (int quad_pt = 0; quad_pt < quad.numQuadPts(); ++quad_pt) {
quad.quadPtCoord(quad_pt, &coord_[0]);
DGBasis::eval(grid_, cell, degree_, &coord_[0], &basis_[0]);
const double w = quad.quadPtWeight(quad_pt);
for (int j = 0; j < num_basis; ++j) {
for (int i = 0; i < num_basis; ++i) {
jac_[j*num_basis + i] += w * basis_[i] * flux_density * basis_[j];
}
}
}
}
// Solve linear equation.
MAT_SIZE_T n = num_basis;
MAT_SIZE_T nrhs = 1;
MAT_SIZE_T lda = num_basis;
std::vector<MAT_SIZE_T> piv(num_basis);
MAT_SIZE_T ldb = num_basis;
MAT_SIZE_T info = 0;
dgesv_(&n, &nrhs, &jac_[0], &lda, &piv[0], &rhs_[0], &ldb, &info);
if (info != 0) {
// Print the local matrix and rhs.
std::cerr << "Failed solving single-cell system Ax = b in cell " << cell
<< " with A = \n";
for (int row = 0; row < n; ++row) {
for (int col = 0; col < n; ++col) {
std::cerr << " " << jac_[row + n*col];
}
std::cerr << '\n';
}
std::cerr << "and b = \n";
for (int row = 0; row < n; ++row) {
std::cerr << " " << rhs_[row] << '\n';
}
THROW("Lapack error: " << info << " encountered in cell " << cell);
}
// The solution ends up in rhs_, so we must copy it.
std::copy(rhs_.begin(), rhs_.end(), tof_coeff_ + num_basis*cell);
}
void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
}
}
} // namespace Opm

93
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp Normal file
View File

@@ -0,0 +1,93 @@
/*
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 OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <vector>
#include <map>
#include <ostream>
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
/// Implements a discontinuous Galerkin solver for
/// (single-phase) time-of-flight using reordering.
/// The equation solved is:
/// v \cdot \grad\tau = \phi
/// where v is the fluid velocity, \tau is time-of-flight and
/// \phi is the porosity. This is a boundary value problem, where
/// \tau is specified to be zero on all inflow boundaries.
/// The user may specify the polynomial degree of the basis function space
/// used, but only degrees 0 and 1 are supported so far.
class TransportModelTracerTofDiscGal : public TransportModelInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTofDiscGal(const UnstructuredGrid& grid);
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Source term. Sign convention is:
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[in] degree Polynomial degree of DG basis functions used.
/// \param[out] tof_coeff Array of time-of-flight solution coefficients.
/// The values are ordered by cell, meaning that
/// the K coefficients corresponding to the first
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
const int degree,
std::vector<double>& tof_coeff);
private:
virtual void solveSingleCell(const int cell);
virtual void solveMultiCell(const int num_cells, const int* cells);
private:
const UnstructuredGrid& grid_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
int degree_;
double* tof_coeff_;
std::vector<double> rhs_; // single-cell right-hand-side
std::vector<double> jac_; // single-cell jacobian
// Below: storage for quantities needed by solveSingleCell().
std::vector<double> coord_;
std::vector<double> basis_;
std::vector<double> basis_nb_;
std::vector<double> grad_basis_;
std::vector<double> velocity_;
};
} // namespace Opm
#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED

26
opm/core/transport/reorder/nlsolvers.c
View File

@@ -26,24 +26,8 @@ SOFTWARE.
#include <assert.h>
#include <math.h>
#ifdef MATLAB_MEX_FILE
#include "nlsolvers.h"
#include <mex.h>
extern int interrupt_signal;
#define print mexPrintf
#define malloc mxMalloc
#define calloc mxCalloc
#define realloc mxRealloc
#define free mxFree
#else /* MATLAB_MEX_FILE */
#define print printf
#include <opm/core/transport/reorder/nlsolvers.h>
#endif /* MATLAB_MEX_FILE */
static const char no_root_str[]=
" In %s:\n"
" With G(%5f) =% 5f, G(%5f) =% 5f, G(x) is not bracketed!\n";
@@ -100,7 +84,7 @@ ridders (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctr
if (G0*G1 > 0)
{
print(no_root_str, "ridder", s0, G0, s1, G1);
printf(no_root_str, "ridder", s0, G0, s1, G1);
return -1.0;
}
@@ -157,7 +141,7 @@ ridders (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *ctr
}
else
{
print("In ridder:\nG0=%10.10f, G1=%10.10f, "
printf("In ridder:\nG0=%10.10f, G1=%10.10f, "
"G3=%10.10f\n", G0, G1, G3);
}
s2 = 0.5*(s0+s1);
@@ -199,7 +183,7 @@ regulafalsi (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl
if (G0*G1 > 0)
{
print(no_root_str, "regulafalsi", s0, G0, s1, G1);
printf(no_root_str, "regulafalsi", s0, G0, s1, G1);
return -1.0;
}
@@ -282,7 +266,7 @@ bisection (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *c
if (G0*G1 > 0.0)
{
print(no_root_str, "bisection", s0, G0, s1, G1);
printf(no_root_str, "bisection", s0, G0, s1, G1);
return -1.0;
}
@@ -312,7 +296,7 @@ bisection (double (*G)(double, void*), void *data, struct NonlinearSolverCtrl *c
}
if (ctrl->iterations >= ctrl->maxiterations)
{
print("Warning: convergence criterion not met\n");
printf("Warning: convergence criterion not met\n");
}
ctrl->residual = Gn;
return sn;

4
opm/core/transport/reorder/reordersequence.cpp
View File

@@ -1,11 +1,11 @@
/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
#include <opm/core/grid.h>
#ifdef MATLAB_MEX_FILE
#include "grid.h"
#include "reordersequence.h"
#include "tarjan.h"
#else
#include <opm/core/grid.h>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/transport/reorder/tarjan.h>
#endif

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

@@ -553,6 +553,7 @@ namespace Opm
{
const int np = wells.number_of_phases;
const int nw = wells.number_of_wells;
ASSERT(int(flow_rates_per_well_cell.size()) == wells.well_connpos[nw]);
phase_flow_per_well.resize(nw * np);
for (int wix = 0; wix < nw; ++wix) {
for (int phase = 0; phase < np; ++phase) {

124
opm/core/utility/miscUtilitiesBlackoil.cpp
View File

@@ -21,7 +21,10 @@
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/grid.h>
#include <opm/core/newwells.h>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <functional>
@@ -31,53 +34,74 @@
namespace Opm
{
/// @brief Computes injected and produced volumes of all phases.
/// @brief Computes injected and produced surface volumes of all phases.
/// Note 1: assumes that only the first phase is injected.
/// Note 2: assumes that transport has been done with an
/// implicit method, i.e. that the current state
/// gives the mobilities used for the preceding timestep.
/// @param[in] props fluid and rock properties.
/// @param[in] p pressure (one value per cell)
/// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all P phases)
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
/// Note 3: Gives surface volume values, not reservoir volumes
/// (as the incompressible version of the function does).
/// Also, assumes that transport_src is given in surface volumes
/// for injector terms!
/// @param[in] props fluid and rock properties.
/// @param[in] state state variables (pressure, sat, surfvol)
/// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
const std::vector<double>& press,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& src,
const BlackoilState& state,
const std::vector<double>& transport_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.");
const int num_cells = transport_src.size();
if (props.numCells() != num_cells) {
THROW("Size of transport_src vector does not match number of cells in props.");
}
const int np = props.numPhases();
if (int(state.saturation().size()) != num_cells*np) {
THROW("Sizes of state vectors do not match number of cells.");
}
const std::vector<double>& press = state.pressure();
const std::vector<double>& s = state.saturation();
const std::vector<double>& z = state.surfacevol();
std::fill(injected, injected + np, 0.0);
std::fill(produced, produced + np, 0.0);
std::vector<double> visc(np);
std::vector<double> mob(np);
std::vector<double> A(np*np);
std::vector<double> prod_resv_phase(np);
std::vector<double> prod_surfvol(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;
if (transport_src[c] > 0.0) {
// Inflowing transport source is a surface volume flux
// for the first phase.
injected[0] += transport_src[c]*dt;
} else if (transport_src[c] < 0.0) {
// Outflowing transport source is a total reservoir
// volume flux.
const double flux = -transport_src[c]*dt;
const double* sat = &s[np*c];
props.relperm(1, sat, &c, &mob[0], 0);
props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0);
double totmob = 0.0;
for (int p = 0; p < np; ++p) {
mob[p] /= visc[p];
totmob += mob[p];
}
std::fill(prod_surfvol.begin(), prod_surfvol.end(), 0.0);
for (int p = 0; p < np; ++p) {
produced[p] += (mob[p]/totmob)*flux;
prod_resv_phase[p] = (mob[p]/totmob)*flux;
for (int q = 0; q < np; ++q) {
prod_surfvol[q] += prod_resv_phase[p]*A[q + np*p];
}
}
for (int p = 0; p < np; ++p) {
produced[p] += prod_surfvol[p];
}
}
}
@@ -251,4 +275,58 @@ namespace Opm
}
}
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (surface volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropertiesInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& transport_src)
{
int nc = props.numCells();
transport_src.clear();
transport_src.resize(nc, 0.0);
// Well contributions.
if (wells) {
const int nw = wells->number_of_wells;
const int np = wells->number_of_phases;
if (np != 2) {
THROW("computeTransportSource() requires a 2 phase case.");
}
std::vector<double> A(np*np);
for (int w = 0; w < nw; ++w) {
const double* comp_frac = wells->comp_frac + np*w;
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
const int perf_cell = wells->well_cells[perf];
double perf_rate = well_state.perfRates()[perf];
if (perf_rate > 0.0) {
// perf_rate is a total inflow reservoir rate, we want a surface water rate.
if (wells->type[w] != INJECTOR) {
std::cout << "**** Warning: crossflow in well "
<< w << " perf " << perf - wells->well_connpos[w]
<< " ignored. Reservoir rate was "
<< perf_rate/Opm::unit::day << " m^3/day." << std::endl;
perf_rate = 0.0;
} else {
ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
perf_rate *= comp_frac[0]; // Water reservoir volume rate.
props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0);
perf_rate *= A[0]; // Water surface volume rate.
}
}
transport_src[perf_cell] += perf_rate;
}
}
}
}
} // namespace Opm

51
opm/core/utility/miscUtilitiesBlackoil.hpp
View File

@@ -22,36 +22,40 @@
#include <vector>
struct UnstructuredGrid;
struct Wells;
namespace Opm
{
class BlackoilPropertiesInterface;
class BlackoilState;
class WellState;
/// @brief Computes injected and produced volumes of all phases.
/// @brief Computes injected and produced surface volumes of all phases.
/// Note 1: assumes that only the first phase is injected.
/// Note 2: assumes that transport has been done with an
/// implicit method, i.e. that the current state
/// gives the mobilities used for the preceding timestep.
/// @param[in] props fluid and rock properties.
/// @param[in] p pressure (one value per cell)
/// @param[in] z surface-volume values (for all P phases)
/// @param[in] s saturation values (for all P phases)
/// @param[in] src if < 0: total outflow, if > 0: first phase inflow.
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
/// Note 3: Gives surface volume values, not reservoir volumes
/// (as the incompressible version of the function does).
/// Also, assumes that transport_src is given in surface volumes
/// for injector terms!
/// @param[in] props fluid and rock properties.
/// @param[in] state state variables (pressure, sat, surfvol)
/// @param[in] transport_src if < 0: total resv outflow, if > 0: first phase surfv inflow
/// @param[in] dt timestep used
/// @param[out] injected must point to a valid array with P elements,
/// where P = s.size()/src.size().
/// @param[out] produced must also point to a valid array with P elements.
void computeInjectedProduced(const BlackoilPropertiesInterface& props,
const std::vector<double>& p,
const std::vector<double>& z,
const std::vector<double>& s,
const std::vector<double>& src,
const BlackoilState& state,
const std::vector<double>& transport_src,
const double dt,
double* injected,
double* produced);
/// @brief Computes total mobility for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
@@ -66,6 +70,7 @@ namespace Opm
const std::vector<double>& s,
std::vector<double>& totmob);
/// @brief Computes total mobility and omega for a set of saturation values.
/// @param[in] props rock and fluid properties
/// @param[in] cells cells with which the saturation values are associated
@@ -131,6 +136,22 @@ namespace Opm
const double* saturation,
double* surfacevol);
/// Compute two-phase transport source terms from well terms.
/// Note: Unlike the incompressible version of this function,
/// this version computes surface volume injection rates,
/// production rates are still total reservoir volumes.
/// \param[in] props Fluid and rock properties.
/// \param[in] wells Wells data structure.
/// \param[in] well_state Well pressures and fluxes.
/// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
/// (+) positive inflow of first (water) phase (surface volume),
/// (-) negative total outflow of both phases (reservoir volume).
void computeTransportSource(const BlackoilPropertiesInterface& props,
const Wells* wells,
const WellState& well_state,
std::vector<double>& transport_src);
} // namespace Opm
#endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED

5
opm/core/utility/parameters/ParameterGroup.cpp
View File

@@ -4,7 +4,7 @@
//
// Created: Tue Jun 2 19:13:17 2009
//
// Author(s): B<EFBFBD>rd Skaflestad <bard.skaflestad@sintef.no>
// Author(s): Bård Skaflestad <bard.skaflestad@sintef.no>
// Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $Date$
@@ -139,9 +139,6 @@ namespace Opm {
}
}
}
#ifdef MATLAB_MEX_FILE
fclose(is);
#endif
}
void ParameterGroup::writeParam(const std::string& param_filename) const {

13
opm/core/wells/WellsGroup.cpp
View File

@@ -572,6 +572,7 @@ namespace Opm
break;
}
const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase);
const double total_reinjected = - total_produced; // Production negative, injection positive
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
@@ -580,11 +581,11 @@ namespace Opm
const double children_guide_rate = children_[i]->injectionGuideRate(true);
#ifdef DIRTY_WELLCTRL_HACK
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
false);
#else
children_[i]->applyInjGroupControl(InjectionSpecification::RATE,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().reinjection_fraction_target_,
false);
#endif
}
@@ -600,15 +601,15 @@ namespace Opm
if (phaseUsage().phase_used[BlackoilPhases::Vapour]) {
total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour);
}
const double my_guide_rate = injectionGuideRate(true);
const double total_reinjected = - total_produced; // Production negative, injection positive
const double my_guide_rate = injectionGuideRate(true);
for (size_t i = 0; i < children_.size(); ++i) {
// Apply for all children.
// Note, we do _not_ want to call the applyProdGroupControl in this object,
// as that would check if we're under group control, something we're not.
const double children_guide_rate = children_[i]->injectionGuideRate(true);
children_[i]->applyInjGroupControl(InjectionSpecification::RESV,
(children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_,
(children_guide_rate / my_guide_rate) * total_reinjected * injSpec().voidage_replacment_fraction_,
false);
}
@@ -863,7 +864,7 @@ namespace Opm
return;
}
// We're a producer, so we need to negate the input
double ntarget = target;
double ntarget = -target;
double distr[3] = { 0.0, 0.0, 0.0 };
const int* phase_pos = phaseUsage().phase_pos;

8
opm/core/wells/WellsManager.cpp
View File

@@ -226,6 +226,14 @@ namespace Opm
/// Construct from existing wells object.
WellsManager::WellsManager(struct Wells* W)
: w_(clone_wells(W))
{
}
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,

9
opm/core/wells/WellsManager.hpp
View File

@@ -41,7 +41,14 @@ namespace Opm
{
public:
/// Default constructor -- no wells.
WellsManager();
WellsManager();
/// Construct from existing wells object.
/// WellsManager is not properly initialised in the sense that the logic to
/// manage control switching does not exist.
///
/// @param[in] W Existing wells object.
WellsManager(struct Wells* W);
/// Construct from input deck and grid.
/// The permeability argument may be zero if the input contain

7
tests/Makefile.am
View File

@@ -1,6 +1,6 @@
AM_CPPFLAGS = \
-I$(top_srcdir) \
$(OPM_BOOST_CPPFLAGS) \
$(OPM_BOOST_CPPFLAGS) ${SUPERLU_CPPFLAGS} \
$(ERT_CPPFLAGS)
AM_LDFLAGS = $(OPM_BOOST_LDFLAGS) $(ERT_LDFLAGS)
@@ -20,9 +20,11 @@ sparsevector_test \
test_cartgrid \
test_column_extract \
test_lapack \
test_read_grid \
test_read_vag \
test_readpolymer \
test_sf2p \
test_wells \
test_writeVtkData \
unit_test
@@ -59,6 +61,9 @@ test_sf2p_SOURCES = test_sf2p.cpp
test_writeVtkData_SOURCES = test_writeVtkData.cpp
test_wells_SOURCES = test_wells.cpp
test_wells_LDADD = $(LDADD) $(BOOST_UNIT_TEST_FRAMEWORK_LIB)
unit_test_SOURCES = unit_test.cpp

54
tests/test_lapack.cpp
View File

@@ -24,34 +24,36 @@
namespace {
struct BandMatrixCoeff
{
BandMatrixCoeff(int N, int ku, int kl) : ku_(ku), kl_(kl), nrow_(2*kl + ku + 1), N_(N) {
BandMatrixCoeff(MAT_SIZE_T N, MAT_SIZE_T ku, MAT_SIZE_T kl)
: ku_(ku), kl_(kl), nrow_(2 * kl + ku + 1), N_(N)
{
}
// compute the position where to store the coefficient of a matrix A_{i,j} (i,j=0,...,N-1)
// in a array which is sent to the band matrix solver of LAPACK.
int operator ()(int i, int j) {
int operator ()(MAT_SIZE_T i, MAT_SIZE_T j) {
return kl_ + ku_ + i - j + j*nrow_;
}
const int ku_;
const int kl_;
const int nrow_;
const int N_;
const MAT_SIZE_T ku_;
const MAT_SIZE_T kl_;
const MAT_SIZE_T nrow_;
const MAT_SIZE_T N_;
};
} //end anonymous namespace
int main()
{
const int N = 5;
const int nrhs = 1;
const MAT_SIZE_T N = 5;
const MAT_SIZE_T nrhs = 1;
double DU[N-1] = { 2.1, -1.0, 1.9, 8.0 };
double D[N] = { 3.0, 2.3, -5.0, -0.9, 7.1 };
double DL[N-1] = { 3.4, 3.6, 7.0, -6.0 };
double B[N] = { 2.7, -0.5, 2.6, 0.6, 2.7 };
// double B[N] = { 2.7, -0.5, 2.6, 0.6, 2.7 };
int info = 0;
MAT_SIZE_T info = 0;
dgtsv_(&N, &nrhs, DL, D, DU, B, &N, &info);
if (info == 0) {
for (int i = 0; i < N; ++i) {
@@ -64,12 +66,12 @@ int main()
std::cout << std::endl;
//test of dgbsv_
int ldb = N;
int lda = N;
std::vector<int> ipiv(N, 0);
MAT_SIZE_T ldb = N;
MAT_SIZE_T lda = N;
std::vector<MAT_SIZE_T> ipiv(N, 0);
std::vector<double> AA(N*N, 0.);
std::vector<double> BB(N, 0.);
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
AA[i + i*N] = 10.;
if (i > 0) {
AA[i + (i - 1)*N] = i;
@@ -78,31 +80,31 @@ int main()
BB[i] = i;
}
for (int i = 0; i < N; ++i) {
for (int j = 0; j < N; ++j) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
for (MAT_SIZE_T j = 0; j < N; ++j) {
std::cout << " " << AA[i + j*N];
}
std::cout << " " << std::endl;
}
std::cout << std::endl;
int kl = 1;
int ku = 1;
int nrowAB = 2*kl + ku + 1;
int ldab = nrowAB;
MAT_SIZE_T kl = 1;
MAT_SIZE_T ku = 1;
MAT_SIZE_T nrowAB = 2*kl + ku + 1;
MAT_SIZE_T ldab = nrowAB;
std::vector<double> AB(nrowAB*N, 0.);
BandMatrixCoeff bmc(N, ku, kl);
// Rewrite AA matrix in band format AB
for (int i = 0; i < N; ++i) {
for (int j = -1; j < 2; ++j) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
for (MAT_SIZE_T j = -1; j < 2; ++j) {
if (i + j > -1 && i + j < N)
AB[bmc(i, i + j)] = AA[i + N*(i + j)];
}
}
for (int i = 0; i < nrowAB; ++i) {
for (int j = 0; j < N; ++j) {
for (MAT_SIZE_T i = 0; i < nrowAB; ++i) {
for (MAT_SIZE_T j = 0; j < N; ++j) {
std::cout << " " << AB[i + j*nrowAB];
}
std::cout << " " << std::endl;
@@ -113,7 +115,7 @@ int main()
dgesv_(&N ,&nrhs, &AA[0], &lda, &ipiv[0], &BB[0], &ldb, &info);
if (info == 0) {
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
std::cout << BB[i] << ' ';
}
std::cout << std::endl;
@@ -121,14 +123,14 @@ int main()
std::cerr << "Something went wrong in debsv_()\n";
}
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
BB[i] = i;
}
dgbsv_(&N, &kl, &ku, &nrhs, &AB[0], &ldab, &ipiv[0], &BB[0], &ldb, &info);
if (info == 0) {
for (int i = 0; i < N; ++i) {
for (MAT_SIZE_T i = 0; i < N; ++i) {
std::cout << BB[i] << ' ';
}
std::cout << std::endl;

25
tests/test_read_grid.c Normal file
View File

@@ -0,0 +1,25 @@
/*
* test_read_grid.c
*
* Created on: Aug 28, 2012
* Author: bska
*/
#include <stdio.h>
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
int
main(void)
{
struct UnstructuredGrid *G1, *G2;
G1 = read_grid("cart_grid_2d.txt");
G2 = create_grid_cart2d(2, 2);
destroy_grid(G2);
destroy_grid(G1);
return 0;
}

205
tests/test_wells.cpp Normal file
View File

@@ -0,0 +1,205 @@
/*
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/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE WellsModuleTest
#include <boost/test/unit_test.hpp>
#include <opm/core/newwells.h>
#include <iostream>
#include <vector>
#include <boost/shared_ptr.hpp>
BOOST_AUTO_TEST_CASE(Construction)
{
const int nphases = 2;
const int nwells = 2;
const int nperfs = 2;
boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
destroy_wells);
if (W) {
int cells[] = { 0, 9 };
double WI = 1.0;
const double ifrac[] = { 1.0, 0.0 };
const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
&WI, "INJECTOR", W.get());
const double pfrac[] = { 0.0, 0.0 };
const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
&WI, "PRODUCER", W.get());
if (ok0 && ok1) {
BOOST_CHECK_EQUAL(W->number_of_phases, nphases);
BOOST_CHECK_EQUAL(W->number_of_wells , nwells );
BOOST_CHECK_EQUAL(W->well_connpos[0], 0);
BOOST_CHECK_EQUAL(W->well_connpos[1], 1);
BOOST_CHECK_EQUAL(W->well_connpos[W->number_of_wells], nperfs);
BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[0]], cells[0]);
BOOST_CHECK_EQUAL(W->well_cells[W->well_connpos[1]], cells[1]);
BOOST_CHECK_EQUAL(W->WI[W->well_connpos[0]], WI);
BOOST_CHECK_EQUAL(W->WI[W->well_connpos[1]], WI);
using std::string;
BOOST_CHECK_EQUAL(string(W->name[0]), string("INJECTOR"));
BOOST_CHECK_EQUAL(string(W->name[1]), string("PRODUCER"));
}
}
}
BOOST_AUTO_TEST_CASE(Controls)
{
const int nphases = 2;
const int nwells = 1;
const int nperfs = 2;
boost::shared_ptr<Wells> W(create_wells(nphases, nwells, nperfs),
destroy_wells);
if (W) {
int cells[] = { 0 , 9 };
double WI [] = { 1.0, 1.0 };
const double ifrac[] = { 1.0, 0.0 };
const bool ok = add_well(INJECTOR, 0.0, nperfs, &ifrac[0], &cells[0],
&WI[0], "INJECTOR", W.get());
if (ok) {
const double distr[] = { 1.0, 0.0 };
const bool ok1 = append_well_controls(BHP, 1, &distr[0],
0, W.get());
const bool ok2 = append_well_controls(SURFACE_RATE, 1,
&distr[0], 0, W.get());
if (ok1 && ok2) {
WellControls* ctrls = W->ctrls[0];
BOOST_CHECK_EQUAL(ctrls->num , 2);
BOOST_CHECK_EQUAL(ctrls->current, -1);
set_current_control(0, 0, W.get());
BOOST_CHECK_EQUAL(ctrls->current, 0);
set_current_control(0, 1, W.get());
BOOST_CHECK_EQUAL(ctrls->current, 1);
BOOST_CHECK_EQUAL(ctrls->type[0], BHP);
BOOST_CHECK_EQUAL(ctrls->type[1], SURFACE_RATE);
BOOST_CHECK_EQUAL(ctrls->target[0], 1.0);
BOOST_CHECK_EQUAL(ctrls->target[1], 1.0);
}
}
}
}
BOOST_AUTO_TEST_CASE(Copy)
{
const int nphases = 2;
const int nwells = 2;
const int nperfs = 2;
boost::shared_ptr<Wells> W1(create_wells(nphases, nwells, nperfs),
destroy_wells);
boost::shared_ptr<Wells> W2;
if (W1) {
int cells[] = { 0, 9 };
const double WI = 1.0;
const double ifrac[] = { 1.0, 0.0 };
const bool ok0 = add_well(INJECTOR, 0.0, 1, &ifrac[0], &cells[0],
&WI, "INJECTOR", W1.get());
const double pfrac[] = { 0.0, 0.0 };
const bool ok1 = add_well(PRODUCER, 0.0, 1, &pfrac[0], &cells[1],
&WI, "PRODUCER", W1.get());
bool ok = ok0 && ok1;
for (int w = 0; ok && (w < W1->number_of_wells); ++w) {
const double distr[] = { 1.0, 0.0 };
const bool okc1 = append_well_controls(BHP, 1, &distr[0],
w, W1.get());
const bool okc2 = append_well_controls(SURFACE_RATE, 1,
&distr[0], w,
W1.get());
ok = okc1 && okc2;
}
if (ok) {
W2.reset(clone_wells(W1.get()), destroy_wells);
}
}
if (W2) {
BOOST_CHECK_EQUAL(W2->number_of_phases, W1->number_of_phases);
BOOST_CHECK_EQUAL(W2->number_of_wells , W1->number_of_wells );
BOOST_CHECK_EQUAL(W2->well_connpos[0] , W1->well_connpos[0] );
for (int w = 0; w < W1->number_of_wells; ++w) {
using std::string;
BOOST_CHECK_EQUAL(string(W2->name[w]), string(W1->name[w]));
BOOST_CHECK_EQUAL( W2->type[w] , W1->type[w] );
BOOST_CHECK_EQUAL(W2->well_connpos[w + 1],
W1->well_connpos[w + 1]);
for (int j = W1->well_connpos[w];
j < W1->well_connpos[w + 1]; ++j) {
BOOST_CHECK_EQUAL(W2->well_cells[j], W1->well_cells[j]);
BOOST_CHECK_EQUAL(W2->WI [j], W1->WI [j]);
}
BOOST_CHECK(W1->ctrls[w] != 0);
BOOST_CHECK(W2->ctrls[w] != 0);
WellControls* c1 = W1->ctrls[w];
WellControls* c2 = W2->ctrls[w];
BOOST_CHECK_EQUAL(c2->num , c1->num );
BOOST_CHECK_EQUAL(c2->current, c1->current);
for (int c = 0; c < c1->num; ++c) {
BOOST_CHECK_EQUAL(c2->type [c], c1->type [c]);
BOOST_CHECK_EQUAL(c2->target[c], c1->target[c]);
for (int p = 0; p < W1->number_of_phases; ++p) {
BOOST_CHECK_EQUAL(c2->distr[c*W1->number_of_phases + p],
c1->distr[c*W1->number_of_phases + p]);
}
}
}
}
}