diff --git a/.gitignore b/.gitignore index 5fe7b14c..a3a8f8a4 100644 --- a/.gitignore +++ b/.gitignore @@ -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 diff --git a/Makefile.am b/Makefile.am index 61d08142..787af5f6 100644 --- a/Makefile.am +++ b/Makefile.am @@ -11,9 +11,10 @@ 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 @@ -23,6 +24,7 @@ lib_libopmcore_la_LDFLAGS = \ -R $(OPM_BOOST_LIBDIR) \ $(OPM_BOOST_LDFLAGS) \ $(ERT_LDFLAGS) +${SUPERLU_LDFLAGS} lib_libopmcore_la_LIBADD = \ $(BOOST_FILESYSTEM_LIB) \ @@ -30,7 +32,7 @@ $(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. @@ -59,80 +61,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 = \ @@ -158,137 +162,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 @@ -299,26 +306,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 diff --git a/README b/README index f0314100..f2dfddcc 100644 --- a/README +++ b/README @@ -56,10 +56,10 @@ sudo add-apt-repository -y ppa:opm/ppa sudo apt-get update # parts of DUNE needed -sudo apt-get install libdune-common-dev libdune-istl-dev +sudo apt-get install libdune-common-dev libdune-istl-dev libdune-grid-dev # libraries necessary for OPM -sudo apt-get install -y libxml0-dev +sudo apt-get install -y libxml2-dev DEPENDENCIES FOR SUSE BASED DISTRIBUTIONS @@ -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 diff --git a/examples/Makefile.am b/examples/Makefile.am index 35fe4e75..145fbbdd 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -3,19 +3,34 @@ 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 +# 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) + # ---------------------------------------------------------------------- # Declare products (i.e., the example programs). # # Please keep the list sorted. noinst_PROGRAMS = \ +compute_tof \ refine_wells \ scaneclipsedeck \ sim_2p_comp_reorder \ @@ -30,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 # ---------------------------------------------------------------------- @@ -45,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 diff --git a/examples/compute_tof.cpp b/examples/compute_tof.cpp new file mode 100644 index 00000000..046dae38 --- /dev/null +++ b/examples/compute_tof.cpp @@ -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 . +*/ + + + +#if HAVE_CONFIG_H +#include "config.h" +#endif // HAVE_CONFIG_H + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + + +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 deck; + boost::scoped_ptr grid; + boost::scoped_ptr props; + boost::scoped_ptr 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("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 porevol; + computePorevolume(*grid->c_grid(), props->porosity(), porevol); + int num_cells = grid->c_grid()->number_of_cells; + std::vector 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("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("pside"); + double pside_pressure = param.get("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 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 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(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; +} diff --git a/m4/opm_core.m4 b/m4/opm_core.m4 index efb75c13..40e1b7c5 100644 --- a/m4/opm_core.m4 +++ b/m4/opm_core.m4 @@ -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. diff --git a/m4/opm_superlu.m4 b/m4/opm_superlu.m4 new file mode 100644 index 00000000..ed42bd41 --- /dev/null +++ b/m4/opm_superlu.m4 @@ -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" +]) diff --git a/opm/core/eclipse/EclipseGridParserHelpers.hpp b/opm/core/eclipse/EclipseGridParserHelpers.hpp index 6ab5fb78..7bc0e003 100644 --- a/opm/core/eclipse/EclipseGridParserHelpers.hpp +++ b/opm/core/eclipse/EclipseGridParserHelpers.hpp @@ -398,10 +398,18 @@ namespace yv.push_back(table[k][i]); } } - // Interpolate - for (int i=0; i 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]; diff --git a/opm/core/fluid/IncompPropertiesInterface.hpp b/opm/core/fluid/IncompPropertiesInterface.hpp index 4f8bfb51..d3aaaa7b 100644 --- a/opm/core/fluid/IncompPropertiesInterface.hpp +++ b/opm/core/fluid/IncompPropertiesInterface.hpp @@ -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; diff --git a/opm/core/grid/cpgpreprocess/preprocess.c b/opm/core/grid/cpgpreprocess/preprocess.c index 84e9e8fe..dc274160 100644 --- a/opm/core/grid/cpgpreprocess/preprocess.c +++ b/opm/core/grid/cpgpreprocess/preprocess.c @@ -590,6 +590,9 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum, c1 = i/2 + nx*(j/2 + ny*(k/2)); c2 = i/2 + nx*(j/2 + ny*((k+1)/2)); + assert (c1 < (nx * ny * nz)); + assert (c2 < (nx * ny * nz)); + if (((actnum == NULL) || (actnum[c1] && actnum[c2])) && (z2 < z1)) { diff --git a/opm/core/newwells.c b/opm/core/newwells.c index 68037201..c28e8566 100644 --- a/opm/core/newwells.c +++ b/opm/core/newwells.c @@ -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; +} diff --git a/opm/core/newwells.h b/opm/core/newwells.h index f16675b6..20cfce61 100644 --- a/opm/core/newwells.h +++ b/opm/core/newwells.h @@ -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 diff --git a/opm/core/simulator/SimulatorCompressibleTwophase.cpp b/opm/core/simulator/SimulatorCompressibleTwophase.cpp index 6fb475e9..e921db87 100644 --- a/opm/core/simulator/SimulatorCompressibleTwophase.cpp +++ b/opm/core/simulator/SimulatorCompressibleTwophase.cpp @@ -407,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 diff --git a/opm/core/simulator/SimulatorIncompTwophase.cpp b/opm/core/simulator/SimulatorIncompTwophase.cpp index 60f63806..911c64ce 100644 --- a/opm/core/simulator/SimulatorIncompTwophase.cpp +++ b/opm/core/simulator/SimulatorIncompTwophase.cpp @@ -485,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 diff --git a/opm/core/transport/reorder/TransportModelInterface.cpp b/opm/core/transport/reorder/TransportModelInterface.cpp index f8e138a3..8a90ae34 100644 --- a/opm/core/transport/reorder/TransportModelInterface.cpp +++ b/opm/core/transport/reorder/TransportModelInterface.cpp @@ -20,6 +20,7 @@ #include #include #include +#include #include #include @@ -31,7 +32,11 @@ void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& g std::vector sequence(grid.number_of_cells); std::vector 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) { diff --git a/opm/core/transport/reorder/TransportModelTracerTof.cpp b/opm/core/transport/reorder/TransportModelTracerTof.cpp new file mode 100644 index 00000000..8d93f248 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.cpp @@ -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 . +*/ + +#include +#include +#include +#include +#include +#include + +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& 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 diff --git a/opm/core/transport/reorder/TransportModelTracerTof.hpp b/opm/core/transport/reorder/TransportModelTracerTof.hpp new file mode 100644 index 00000000..270ac328 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTof.hpp @@ -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 . +*/ + +#ifndef OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED + +#include +#include +#include +#include +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& 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 diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp new file mode 100644 index 00000000..2d58ec91 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp @@ -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 . +*/ + +#include +#include +#include +#include +#include +#include +#include + +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 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& 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 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 diff --git a/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp new file mode 100644 index 00000000..0285e1f0 --- /dev/null +++ b/opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp @@ -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 . +*/ + +#ifndef OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED +#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED + +#include +#include +#include +#include +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& 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 rhs_; // single-cell right-hand-side + std::vector jac_; // single-cell jacobian + // Below: storage for quantities needed by solveSingleCell(). + std::vector coord_; + std::vector basis_; + std::vector basis_nb_; + std::vector grad_basis_; + std::vector velocity_; + }; + +} // namespace Opm + +#endif // OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 53e99773..73e10ac9 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -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) { diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index afa17f12..d407c033 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -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; diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 95cb0369..8aaca9e2 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -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, diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 6df65ec4..19289c5d 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -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 diff --git a/tests/Makefile.am b/tests/Makefile.am index 402c4e1b..bcf9decb 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -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) @@ -24,6 +24,7 @@ test_read_grid \ test_read_vag \ test_readpolymer \ test_sf2p \ +test_wells \ test_writeVtkData \ unit_test @@ -60,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 diff --git a/tests/test_wells.cpp b/tests/test_wells.cpp new file mode 100644 index 00000000..4f8ad7bb --- /dev/null +++ b/tests/test_wells.cpp @@ -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 . +*/ + +#include + +#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 + +#include + +#include +#include +#include + +BOOST_AUTO_TEST_CASE(Construction) +{ + const int nphases = 2; + const int nwells = 2; + const int nperfs = 2; + + boost::shared_ptr 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 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 W1(create_wells(nphases, nwells, nperfs), + destroy_wells); + boost::shared_ptr 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]); + } + } + } + } +}