diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index e199d022..f4f275da 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -29,128 +29,123 @@ # originally generated with the command: # find opm -name '*.c*' -printf '\t%p\n' | sort list (APPEND MAIN_SOURCE_FILES - opm/core/grid/GridHelpers.cpp - opm/core/grid/GridManager.cpp - opm/core/grid/GridUtilities.cpp - opm/core/grid/grid.c + opm/core/flowdiagnostics/AnisotropicEikonal.cpp + opm/core/flowdiagnostics/DGBasis.cpp + opm/core/flowdiagnostics/FlowDiagnostics.cpp + opm/core/flowdiagnostics/TofDiscGalReorder.cpp + opm/core/flowdiagnostics/TofReorder.cpp + opm/core/grid/GridHelpers.cpp + opm/core/grid/GridManager.cpp + opm/core/grid/GridUtilities.cpp + 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/grid/grid.c opm/core/grid/grid_equal.cpp - 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/io/eclipse/EclipseGridInspector.cpp - opm/core/io/eclipse/EclipseWriter.cpp - opm/core/io/eclipse/EclipseReader.cpp - opm/core/io/eclipse/EclipseWriteRFTHandler.cpp - opm/core/io/eclipse/writeECLData.cpp - opm/core/io/OutputWriter.cpp - opm/core/io/vag/vag.cpp - opm/core/io/vtk/writeVtkData.cpp - opm/core/linalg/LinearSolverFactory.cpp - opm/core/linalg/LinearSolverInterface.cpp - opm/core/linalg/LinearSolverIstl.cpp - opm/core/linalg/LinearSolverUmfpack.cpp - opm/core/linalg/LinearSolverPetsc.cpp - opm/core/linalg/call_umfpack.c - opm/core/linalg/sparse_sys.c - opm/core/pressure/CompressibleTpfa.cpp - opm/core/pressure/FlowBCManager.cpp - opm/core/pressure/IncompTpfa.cpp - opm/core/pressure/IncompTpfaSinglePhase.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/TransTpfa.cpp - opm/core/pressure/tpfa/trans_tpfa.c - opm/core/pressure/legacy_well.c - opm/core/props/BlackoilPropertiesBasic.cpp - opm/core/props/BlackoilPropertiesFromDeck.cpp - opm/core/props/IncompPropertiesBasic.cpp - opm/core/props/IncompPropertiesFromDeck.cpp - opm/core/props/IncompPropertiesSinglePhase.cpp - opm/core/props/pvt/BlackoilPvtProperties.cpp - opm/core/props/pvt/PvtPropertiesBasic.cpp - opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp - opm/core/props/pvt/PvtDead.cpp - opm/core/props/pvt/PvtDeadSpline.cpp - opm/core/props/pvt/PvtInterface.cpp - opm/core/props/pvt/PvtLiveGas.cpp - opm/core/props/pvt/PvtLiveOil.cpp - opm/core/props/rock/RockBasic.cpp - opm/core/props/rock/RockCompressibility.cpp - opm/core/props/rock/RockFromDeck.cpp - opm/core/props/satfunc/SaturationPropsBasic.cpp - opm/core/props/satfunc/SaturationPropsFromDeck.cpp - opm/core/props/satfunc/RelpermDiagnostics.cpp - opm/core/simulator/AdaptiveSimulatorTimer.cpp - opm/core/simulator/BlackoilState.cpp - opm/core/simulator/TimeStepControl.cpp - opm/core/simulator/SimulatorCompressibleTwophase.cpp - opm/core/simulator/SimulatorIncompTwophase.cpp - opm/core/simulator/SimulatorOutput.cpp - opm/core/simulator/SimulatorReport.cpp - opm/core/simulator/SimulatorState.cpp - opm/core/simulator/SimulatorTimer.cpp - opm/core/flowdiagnostics/AnisotropicEikonal.cpp - opm/core/flowdiagnostics/DGBasis.cpp - opm/core/flowdiagnostics/FlowDiagnostics.cpp - opm/core/flowdiagnostics/TofReorder.cpp - opm/core/flowdiagnostics/TofDiscGalReorder.cpp - opm/core/transport/TransportSolverTwophaseInterface.cpp - opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp - opm/core/transport/implicit/transport_source.c - opm/core/transport/minimal/spu_explicit.c - opm/core/transport/minimal/spu_implicit.c - opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp - opm/core/transport/reorder/ReorderSolverInterface.cpp - opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp - opm/core/transport/reorder/reordersequence.cpp - opm/core/transport/reorder/tarjan.c - opm/core/utility/compressedToCartesian.cpp - opm/core/utility/Event.cpp - opm/core/utility/MonotCubicInterpolator.cpp - opm/core/utility/StopWatch.cpp - opm/core/utility/VelocityInterpolation.cpp - opm/core/utility/WachspressCoord.cpp - opm/core/utility/miscUtilities.cpp - opm/core/utility/miscUtilitiesBlackoil.cpp - opm/core/utility/NullStream.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/parameters/tinyxml/xmltest.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 - opm/core/wells/wells.c - opm/core/wells/well_controls.c + opm/core/io/OutputWriter.cpp + opm/core/io/eclipse/EclipseGridInspector.cpp + opm/core/io/eclipse/EclipseReader.cpp + opm/core/io/eclipse/EclipseWriteRFTHandler.cpp + opm/core/io/eclipse/EclipseWriter.cpp + opm/core/io/eclipse/writeECLData.cpp + opm/core/io/vag/vag.cpp + opm/core/io/vtk/writeVtkData.cpp + opm/core/linalg/LinearSolverFactory.cpp + opm/core/linalg/LinearSolverInterface.cpp + opm/core/linalg/LinearSolverIstl.cpp + opm/core/linalg/LinearSolverPetsc.cpp + opm/core/linalg/LinearSolverUmfpack.cpp + opm/core/linalg/call_umfpack.c + opm/core/linalg/sparse_sys.c + opm/core/pressure/CompressibleTpfa.cpp + opm/core/pressure/FlowBCManager.cpp + opm/core/pressure/IncompTpfa.cpp + opm/core/pressure/IncompTpfaSinglePhase.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/legacy_well.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/TransTpfa.cpp + 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/props/BlackoilPropertiesBasic.cpp + opm/core/props/BlackoilPropertiesFromDeck.cpp + opm/core/props/IncompPropertiesBasic.cpp + opm/core/props/IncompPropertiesFromDeck.cpp + opm/core/props/IncompPropertiesSinglePhase.cpp + opm/core/props/pvt/PvtPropertiesBasic.cpp + opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp + opm/core/props/rock/RockBasic.cpp + opm/core/props/rock/RockCompressibility.cpp + opm/core/props/rock/RockFromDeck.cpp + opm/core/props/satfunc/RelpermDiagnostics.cpp + opm/core/props/satfunc/SaturationPropsBasic.cpp + opm/core/props/satfunc/SaturationPropsFromDeck.cpp + opm/core/simulator/AdaptiveSimulatorTimer.cpp + opm/core/simulator/BlackoilState.cpp + opm/core/simulator/SimulatorCompressibleTwophase.cpp + opm/core/simulator/SimulatorIncompTwophase.cpp + opm/core/simulator/SimulatorOutput.cpp + opm/core/simulator/SimulatorReport.cpp + opm/core/simulator/SimulatorState.cpp + opm/core/simulator/SimulatorTimer.cpp + opm/core/simulator/TimeStepControl.cpp + opm/core/transport/TransportSolverTwophaseInterface.cpp + opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp + opm/core/transport/implicit/transport_source.c + opm/core/transport/minimal/spu_explicit.c + opm/core/transport/minimal/spu_implicit.c + opm/core/transport/reorder/ReorderSolverInterface.cpp + opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp + opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp + opm/core/transport/reorder/reordersequence.cpp + opm/core/transport/reorder/tarjan.c + opm/core/utility/Event.cpp + opm/core/utility/MonotCubicInterpolator.cpp + opm/core/utility/NullStream.cpp + opm/core/utility/StopWatch.cpp + opm/core/utility/VelocityInterpolation.cpp + opm/core/utility/WachspressCoord.cpp + opm/core/utility/compressedToCartesian.cpp + opm/core/utility/extractPvtTableIndex.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/parameters/tinyxml/xmltest.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 + opm/core/wells/well_controls.c + opm/core/wells/wells.c ) # originally generated with the command: @@ -182,7 +177,6 @@ list (APPEND TEST_SOURCE_FILES tests/test_linearsolver.cpp tests/test_parallel_linearsolver.cpp tests/test_param.cpp - tests/test_blackoilfluid.cpp tests/test_satfunc.cpp tests/test_shadow.cpp tests/test_equil.cpp @@ -279,191 +273,185 @@ list (APPEND PROGRAM_SOURCE_FILES # originally generated with the command: # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort list (APPEND PUBLIC_HEADER_FILES - opm/core/doxygen_main.hpp - opm/core/grid.h - opm/core/grid/CellQuadrature.hpp - opm/core/grid/ColumnExtract.hpp - opm/core/grid/FaceQuadrature.hpp - opm/core/grid/GridHelpers.hpp - opm/core/grid/GridManager.hpp - opm/core/grid/GridUtilities.hpp - opm/core/grid/MinpvProcessor.hpp - opm/core/grid/PinchProcessor.hpp - 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/preprocess.h - opm/core/grid/cpgpreprocess/uniquepoints.h - opm/core/io/eclipse/CornerpointChopper.hpp - opm/core/io/eclipse/EclipseIOUtil.hpp - opm/core/io/eclipse/EclipseGridInspector.hpp - opm/core/io/eclipse/EclipseUnits.hpp - opm/core/io/eclipse/EclipseWriter.hpp - opm/core/io/eclipse/EclipseReader.hpp - opm/core/io/eclipse/EclipseWriteRFTHandler.hpp - opm/core/io/eclipse/writeECLData.hpp - opm/core/io/OutputWriter.hpp - opm/core/io/vag/vag.hpp - opm/core/io/vtk/writeVtkData.hpp - opm/core/linalg/LinearSolverFactory.hpp - opm/core/linalg/LinearSolverInterface.hpp - opm/core/linalg/LinearSolverIstl.hpp - opm/core/linalg/LinearSolverUmfpack.hpp - opm/core/linalg/LinearSolverPetsc.hpp - opm/core/linalg/ParallelIstlInformation.hpp - opm/core/linalg/blas_lapack.h - opm/core/linalg/call_umfpack.h - opm/core/linalg/sparse_sys.h - opm/core/wells.h - opm/core/well_controls.h - opm/core/pressure/CompressibleTpfa.hpp - opm/core/pressure/FlowBCManager.hpp - opm/core/pressure/IncompTpfa.hpp - opm/core/pressure/IncompTpfaSinglePhase.hpp - opm/core/pressure/flow_bc.h - opm/core/pressure/fsh.h - opm/core/pressure/fsh_common_impl.h - opm/core/pressure/legacy_well.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/TransTpfa.hpp - opm/core/pressure/tpfa/TransTpfa_impl.hpp - opm/core/pressure/tpfa/trans_tpfa.h - opm/core/props/BlackoilPhases.hpp - opm/core/props/BlackoilPropertiesBasic.hpp - opm/core/props/BlackoilPropertiesFromDeck.hpp - opm/core/props/BlackoilPropertiesInterface.hpp - opm/core/props/IncompPropertiesBasic.hpp - opm/core/props/IncompPropertiesFromDeck.hpp - opm/core/props/IncompPropertiesInterface.hpp - opm/core/props/IncompPropertiesShadow.hpp - opm/core/props/IncompPropertiesShadow_impl.hpp - opm/core/props/IncompPropertiesSinglePhase.hpp - opm/core/props/phaseUsageFromDeck.hpp - opm/core/props/pvt/BlackoilPvtProperties.hpp - opm/core/props/pvt/PvtPropertiesBasic.hpp - opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp - opm/core/props/pvt/PvtConstCompr.hpp - opm/core/props/pvt/PvtDead.hpp - opm/core/props/pvt/PvtDeadSpline.hpp - opm/core/props/pvt/PvtInterface.hpp - opm/core/props/pvt/PvtLiveGas.hpp - opm/core/props/pvt/PvtLiveOil.hpp - opm/core/props/pvt/ThermalWaterPvtWrapper.hpp - opm/core/props/pvt/ThermalOilPvtWrapper.hpp - opm/core/props/pvt/ThermalGasPvtWrapper.hpp - opm/core/props/rock/RockBasic.hpp - opm/core/props/rock/RockCompressibility.hpp - opm/core/props/rock/RockFromDeck.hpp - opm/core/props/satfunc/SaturationPropsBasic.hpp - opm/core/props/satfunc/SaturationPropsFromDeck.hpp - opm/core/props/satfunc/SaturationPropsInterface.hpp - opm/core/props/satfunc/RelpermDiagnostics.hpp + opm/core/doxygen_main.hpp + opm/core/flowdiagnostics/AnisotropicEikonal.hpp + opm/core/flowdiagnostics/DGBasis.hpp + opm/core/flowdiagnostics/FlowDiagnostics.hpp + opm/core/flowdiagnostics/TofDiscGalReorder.hpp + opm/core/flowdiagnostics/TofReorder.hpp + opm/core/grid.h + opm/core/grid/CellQuadrature.hpp + opm/core/grid/ColumnExtract.hpp + opm/core/grid/FaceQuadrature.hpp + opm/core/grid/GridHelpers.hpp + opm/core/grid/GridManager.hpp + opm/core/grid/GridUtilities.hpp + opm/core/grid/MinpvProcessor.hpp + opm/core/grid/PinchProcessor.hpp + 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/preprocess.h + opm/core/grid/cpgpreprocess/uniquepoints.h + opm/core/io/OutputWriter.hpp + opm/core/io/eclipse/CornerpointChopper.hpp + opm/core/io/eclipse/EclipseGridInspector.hpp + opm/core/io/eclipse/EclipseIOUtil.hpp + opm/core/io/eclipse/EclipseReader.hpp + opm/core/io/eclipse/EclipseUnits.hpp + opm/core/io/eclipse/EclipseWriteRFTHandler.hpp + opm/core/io/eclipse/EclipseWriter.hpp + opm/core/io/eclipse/writeECLData.hpp + opm/core/io/vag/vag.hpp + opm/core/io/vtk/writeVtkData.hpp + opm/core/linalg/LinearSolverFactory.hpp + opm/core/linalg/LinearSolverInterface.hpp + opm/core/linalg/LinearSolverIstl.hpp + opm/core/linalg/LinearSolverPetsc.hpp + opm/core/linalg/LinearSolverUmfpack.hpp + opm/core/linalg/ParallelIstlInformation.hpp + opm/core/linalg/blas_lapack.h + opm/core/linalg/call_umfpack.h + opm/core/linalg/sparse_sys.h + opm/core/pressure/CompressibleTpfa.hpp + opm/core/pressure/FlowBCManager.hpp + opm/core/pressure/IncompTpfa.hpp + opm/core/pressure/IncompTpfaSinglePhase.hpp + opm/core/pressure/flow_bc.h + opm/core/pressure/fsh.h + opm/core/pressure/fsh_common_impl.h + opm/core/pressure/legacy_well.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/TransTpfa.hpp + opm/core/pressure/tpfa/TransTpfa_impl.hpp + 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/props/BlackoilPhases.hpp + opm/core/props/BlackoilPropertiesBasic.hpp + opm/core/props/BlackoilPropertiesFromDeck.hpp + opm/core/props/BlackoilPropertiesInterface.hpp + opm/core/props/IncompPropertiesBasic.hpp + opm/core/props/IncompPropertiesFromDeck.hpp + opm/core/props/IncompPropertiesInterface.hpp + opm/core/props/IncompPropertiesShadow.hpp + opm/core/props/IncompPropertiesShadow_impl.hpp + opm/core/props/IncompPropertiesSinglePhase.hpp + opm/core/props/phaseUsageFromDeck.hpp + opm/core/props/pvt/PvtPropertiesBasic.hpp + opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp + opm/core/props/pvt/ThermalGasPvtWrapper.hpp + opm/core/props/pvt/ThermalOilPvtWrapper.hpp + opm/core/props/pvt/ThermalWaterPvtWrapper.hpp + opm/core/props/rock/RockBasic.hpp + opm/core/props/rock/RockCompressibility.hpp + opm/core/props/rock/RockFromDeck.hpp + opm/core/props/satfunc/RelpermDiagnostics.hpp + opm/core/props/satfunc/SaturationPropsBasic.hpp + opm/core/props/satfunc/SaturationPropsFromDeck.hpp + opm/core/props/satfunc/SaturationPropsInterface.hpp opm/core/props/satfunc/RelpermDiagnostics_impl.hpp - opm/core/simulator/AdaptiveSimulatorTimer.hpp - opm/core/simulator/AdaptiveTimeStepping.hpp - opm/core/simulator/AdaptiveTimeStepping_impl.hpp - opm/core/simulator/BlackoilState.hpp - opm/core/simulator/BlackoilStateToFluidState.hpp - opm/core/simulator/EquilibrationHelpers.hpp - opm/core/simulator/ExplicitArraysFluidState.hpp - opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp - opm/core/simulator/TimeStepControl.hpp - opm/core/simulator/SimulatorCompressibleTwophase.hpp - opm/core/simulator/SimulatorIncompTwophase.hpp - opm/core/simulator/SimulatorOutput.hpp - opm/core/simulator/SimulatorReport.hpp - opm/core/simulator/SimulatorState.hpp - opm/core/simulator/SimulatorTimerInterface.hpp - opm/core/simulator/SimulatorTimer.hpp - opm/core/simulator/TimeStepControlInterface.hpp - opm/core/simulator/TwophaseState.hpp - opm/core/simulator/TwophaseState_impl.hpp - opm/core/simulator/WellState.hpp - opm/core/simulator/initState.hpp - opm/core/simulator/initState_impl.hpp - opm/core/simulator/initStateEquil.hpp - opm/core/simulator/initStateEquil_impl.hpp - opm/core/flowdiagnostics/AnisotropicEikonal.hpp - opm/core/flowdiagnostics/DGBasis.hpp - opm/core/flowdiagnostics/FlowDiagnostics.hpp - opm/core/flowdiagnostics/TofReorder.hpp - opm/core/flowdiagnostics/TofDiscGalReorder.hpp - opm/core/transport/TransportSolverTwophaseInterface.hpp - opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp - opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp - opm/core/transport/implicit/ImplicitAssembly.hpp - opm/core/transport/implicit/ImplicitTransport.hpp - opm/core/transport/implicit/JacobianSystem.hpp - opm/core/transport/implicit/NormSupport.hpp - opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp - opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp - opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp - opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp - opm/core/transport/implicit/transport_source.h - opm/core/transport/minimal/spu_explicit.h - opm/core/transport/minimal/spu_implicit.h - opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp - opm/core/transport/reorder/ReorderSolverInterface.hpp - opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp - opm/core/transport/reorder/reordersequence.h - opm/core/transport/reorder/tarjan.h - opm/core/utility/Average.hpp - opm/core/utility/CompressedPropertyAccess.hpp - opm/core/utility/compressedToCartesian.hpp - opm/core/utility/DataMap.hpp - opm/core/utility/Event.hpp - opm/core/utility/Event_impl.hpp - opm/core/utility/Factory.hpp - opm/core/utility/MonotCubicInterpolator.hpp - opm/core/utility/NonuniformTableLinear.hpp - opm/core/utility/NullStream.hpp - opm/core/utility/RegionMapping.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/VelocityInterpolation.hpp - opm/core/utility/WachspressCoord.hpp - opm/core/utility/buildUniformMonotoneTable.hpp - opm/core/utility/have_boost_redef.hpp - opm/core/utility/linearInterpolation.hpp - opm/core/utility/miscUtilities.hpp - opm/core/utility/miscUtilities_impl.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/share_obj.hpp - opm/core/utility/thresholdPressures.hpp - 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 - opm/core/wells/WellsManager_impl.hpp + opm/core/simulator/AdaptiveSimulatorTimer.hpp + opm/core/simulator/AdaptiveTimeStepping.hpp + opm/core/simulator/AdaptiveTimeStepping_impl.hpp + opm/core/simulator/BlackoilState.hpp + opm/core/simulator/BlackoilStateToFluidState.hpp + opm/core/simulator/EquilibrationHelpers.hpp + opm/core/simulator/ExplicitArraysFluidState.hpp + opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp + opm/core/simulator/SimulatorCompressibleTwophase.hpp + opm/core/simulator/SimulatorIncompTwophase.hpp + opm/core/simulator/SimulatorOutput.hpp + opm/core/simulator/SimulatorReport.hpp + opm/core/simulator/SimulatorState.hpp + opm/core/simulator/SimulatorTimer.hpp + opm/core/simulator/SimulatorTimerInterface.hpp + opm/core/simulator/TimeStepControl.hpp + opm/core/simulator/TimeStepControlInterface.hpp + opm/core/simulator/TwophaseState.hpp + opm/core/simulator/TwophaseState_impl.hpp + opm/core/simulator/WellState.hpp + opm/core/simulator/initState.hpp + opm/core/simulator/initStateEquil.hpp + opm/core/simulator/initStateEquil_impl.hpp + opm/core/simulator/initState_impl.hpp + opm/core/transport/TransportSolverTwophaseInterface.hpp + opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp + opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp + opm/core/transport/implicit/ImplicitAssembly.hpp + opm/core/transport/implicit/ImplicitTransport.hpp + opm/core/transport/implicit/JacobianSystem.hpp + opm/core/transport/implicit/NormSupport.hpp + opm/core/transport/implicit/SimpleFluid2pWrappingProps.hpp + opm/core/transport/implicit/SimpleFluid2pWrappingProps_impl.hpp + opm/core/transport/implicit/SinglePointUpwindTwoPhase.hpp + opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp + opm/core/transport/implicit/transport_source.h + opm/core/transport/minimal/spu_explicit.h + opm/core/transport/minimal/spu_implicit.h + opm/core/transport/reorder/ReorderSolverInterface.hpp + opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp + opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp + opm/core/transport/reorder/reordersequence.h + opm/core/transport/reorder/tarjan.h + opm/core/utility/Average.hpp + opm/core/utility/CompressedPropertyAccess.hpp + opm/core/utility/DataMap.hpp + opm/core/utility/Event.hpp + opm/core/utility/Event_impl.hpp + opm/core/utility/Factory.hpp + opm/core/utility/MonotCubicInterpolator.hpp + opm/core/utility/NonuniformTableLinear.hpp + opm/core/utility/NullStream.hpp + opm/core/utility/RegionMapping.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/VelocityInterpolation.hpp + opm/core/utility/WachspressCoord.hpp + opm/core/utility/buildUniformMonotoneTable.hpp + opm/core/utility/compressedToCartesian.hpp + opm/core/utility/extractPvtTableIndex.hpp + opm/core/utility/have_boost_redef.hpp + opm/core/utility/linearInterpolation.hpp + opm/core/utility/miscUtilities.hpp + opm/core/utility/miscUtilitiesBlackoil.hpp + opm/core/utility/miscUtilities_impl.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/share_obj.hpp + opm/core/utility/thresholdPressures.hpp + opm/core/well_controls.h + opm/core/wells.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 + opm/core/wells/WellsManager_impl.hpp ) diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index 8807650a..4ca597a4 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -20,14 +20,15 @@ #include "config.h" #include #include +#include #include #include +#include #include #include namespace Opm { - BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState, const UnstructuredGrid& grid, @@ -133,16 +134,15 @@ namespace Opm if (init_rock){ rock_.init(eclState, number_of_cells, global_cell, cart_dims); } - pvt_.init(deck, eclState, /*numSamples=*/0); + phaseUsage_ = phaseUsageFromDeck(deck); + initSurfaceDensities_(deck); + oilPvt_.initFromDeck(deck, eclState); + gasPvt_.initFromDeck(deck, eclState); + waterPvt_.initFromDeck(deck, eclState); SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); ptr->init(phaseUsageFromDeck(deck), materialLawManager); satprops_.reset(ptr); - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } } inline void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr deck, @@ -162,8 +162,11 @@ namespace Opm rock_.init(eclState, number_of_cells, global_cell, cart_dims); } - const int pvt_samples = param.getDefault("pvt_tab_size", -1); - pvt_.init(deck, eclState, pvt_samples); + phaseUsage_ = phaseUsageFromDeck(deck); + initSurfaceDensities_(deck); + oilPvt_.initFromDeck(deck, eclState); + gasPvt_.initFromDeck(deck, eclState); + waterPvt_.initFromDeck(deck, eclState); // Unfortunate lack of pointer smartness here... std::string threephase_model = param.getDefault("threephase_model", "gwseg"); @@ -175,18 +178,12 @@ namespace Opm = new SaturationPropsFromDeck(); ptr->init(phaseUsageFromDeck(deck), materialLawManager); satprops_.reset(ptr); - - if (pvt_.numPhases() != satprops_->numPhases()) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" - << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); - } } BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck() { } - /// \return D, the number of spatial dimensions. int BlackoilPropertiesFromDeck::numDimensions() const { @@ -219,13 +216,13 @@ namespace Opm /// \return P, the number of phases (also the number of components). int BlackoilPropertiesFromDeck::numPhases() const { - return pvt_.numPhases(); + return phaseUsage_.num_phases; } /// \return Object describing the active phases. PhaseUsage BlackoilPropertiesFromDeck::phaseUsage() const { - return pvt_.phaseUsage(); + return phaseUsage_; } /// \param[in] n Number of data points. @@ -244,16 +241,45 @@ namespace Opm double* mu, double* dmudp) const { - if (dmudp) { - OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented."); - } else { - const int *cellPvtTableIdx = cellPvtRegionIndex(); - assert(cellPvtTableIdx != 0); - std::vector pvtTableIdx(n); - for (int i = 0; i < n; ++ i) - pvtTableIdx[i] = cellPvtTableIdx[cells[i]]; + const auto& pu = phaseUsage(); - pvt_.mu(n, &pvtTableIdx[0], p, T, z, mu); + enum PressureEvalTag {}; + typedef Opm::LocalAd::Evaluation LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + LadEval RsLad = 0.0; + LadEval RvLad = 0.0; + LadEval muLad = 0.0; + + pLad.derivatives[0] = 1.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad.value = p[i]; + TLad.value = T[i]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + muLad = waterPvt_.viscosity(pvtRegionIdx, TLad, pLad); + int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]; + mu[offset] = muLad.value; + dmudp[offset] = muLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + muLad = oilPvt_.viscosity(pvtRegionIdx, TLad, pLad, RsLad); + int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]; + mu[offset] = muLad.value; + dmudp[offset] = muLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + muLad = gasPvt_.viscosity(pvtRegionIdx, TLad, pLad, RvLad); + int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]; + mu[offset] = muLad.value; + dmudp[offset] = muLad.derivatives[0]; + } } } @@ -278,27 +304,23 @@ namespace Opm { const int np = numPhases(); - const int *cellPvtTableIdx = cellPvtRegionIndex(); - std::vector pvtTableIdx(n); - for (int i = 0; i < n; ++ i) - pvtTableIdx[i] = cellPvtTableIdx[cells[i]]; - B_.resize(n*np); R_.resize(n*np); if (dAdp) { dB_.resize(n*np); dR_.resize(n*np); - pvt_.dBdp(n, &pvtTableIdx[0], p, T, z, &B_[0], &dB_[0]); - pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]); + + this->compute_dBdp_(n, p, T, z, cells, &B_[0], &dB_[0]); + this->compute_dRdp_(n, p, T, z, cells, &R_[0], &dR_[0]); } else { - pvt_.B(n, &pvtTableIdx[0], p, T, z, &B_[0]); - pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]); + this->compute_B_(n, p, T, z, cells, &B_[0]); + this->compute_R_(n, p, T, z, cells, &R_[0]); } - const int* phase_pos = pvt_.phasePosition(); - bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] && - pvt_.phaseUsed()[BlackoilPhases::Vapour]; - const int o = phase_pos[BlackoilPhases::Liquid]; - const int g = phase_pos[BlackoilPhases::Vapour]; + const auto& pu = phaseUsage(); + bool oil_and_gas = pu.phase_pos[BlackoilPhases::Liquid] && + pu.phase_pos[BlackoilPhases::Vapour]; + const int o = pu.phase_pos[BlackoilPhases::Liquid]; + const int g = pu.phase_pos[BlackoilPhases::Vapour]; // Compute A matrix // #pragma omp parallel for @@ -360,6 +382,276 @@ namespace Opm } } + void BlackoilPropertiesFromDeck::compute_B_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B) const + { + const auto& pu = phaseUsage(); + + typedef double LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + LadEval RsLad = 0.0; + LadEval RvLad = 0.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad = p[i]; + TLad = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + + B[waterOffset] = BLad; + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + double currentRs = 0.0; + double maxRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad); + } + LadEval BLad; + if (currentRs >= maxRs) { + BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RsLad = currentRs; + BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad); + } + + B[oilOffset] = BLad; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + double currentRv = 0.0; + double maxRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad); + } + LadEval BLad; + if (currentRv >= maxRv) { + BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RvLad = currentRv; + BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad); + } + + B[gasOffset] = BLad; + } + } + } + + void BlackoilPropertiesFromDeck::compute_dBdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B, + double* dBdp) const + { + const auto& pu = phaseUsage(); + + enum PressureEvalTag {}; + typedef Opm::LocalAd::Evaluation LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + LadEval RsLad = 0.0; + LadEval RvLad = 0.0; + + pLad.derivatives[0] = 1.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad.value = p[i]; + TLad.value = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + + B[waterOffset] = BLad.value; + dBdp[waterOffset] = BLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + double currentRs = 0.0; + double maxRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad.value, pLad.value); + } + LadEval BLad; + if (currentRs >= maxRs) { + BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RsLad.value = currentRs; + BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad); + } + + B[oilOffset] = BLad.value; + dBdp[oilOffset] = BLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + double currentRv = 0.0; + double maxRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad.value, pLad.value); + } + LadEval BLad; + if (currentRv >= maxRv) { + BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad); + } + else { + RvLad.value = currentRv; + BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad); + } + + B[gasOffset] = BLad.value; + dBdp[gasOffset] = BLad.derivatives[0]; + } + } + } + + void BlackoilPropertiesFromDeck::compute_R_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R) const + { + const auto& pu = phaseUsage(); + + typedef double LadEval; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad = p[i]; + TLad = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + R[waterOffset] = 0.0; // water is always immiscible! + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad); + + double currentRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + } + + RsSatLad = std::min(RsSatLad, currentRs); + + R[oilOffset] = RsSatLad; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad); + + double currentRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + } + + RvSatLad = std::min(RvSatLad, currentRv); + + R[gasOffset] = RvSatLad; + } + } + } + + void BlackoilPropertiesFromDeck::compute_dRdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R, + double* dRdp) const + { + const auto& pu = phaseUsage(); + + enum PressureEvalTag {}; + typedef Opm::LocalAd::Evaluation LadEval; + typedef Opm::MathToolbox Toolbox; + + LadEval pLad = 0.0; + LadEval TLad = 0.0; + + pLad.derivatives[0] = 1.0; + + for (int i = 0; i < n; ++ i) { + int cellIdx = cells[i]; + int pvtRegionIdx = cellPvtRegionIdx_[cellIdx]; + pLad.value = p[i]; + TLad.value = T[i]; + + int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid]; + int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour]; + int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + R[waterOffset] = 0.0; // water is always immiscible! + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad); + + LadEval currentRs = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) { + currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset]; + } + + RsSatLad = Toolbox::min(RsSatLad, currentRs); + + R[oilOffset] = RsSatLad.value; + dRdp[oilOffset] = RsSatLad.derivatives[0]; + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad); + + LadEval currentRv = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) { + currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset]; + } + + RvSatLad = Toolbox::min(RvSatLad, currentRv); + + R[gasOffset] = RvSatLad.value; + dRdp[gasOffset] = RvSatLad.derivatives[0]; + } + } + } + /// \param[in] n Number of data points. /// \param[in] A Array of nP^2 values, where the P^2 values for a cell give the /// matrix A = RB^{-1} which relates z to u by z = Au. The matrices @@ -376,8 +668,7 @@ namespace Opm // #pragma omp parallel for for (int i = 0; i < n; ++i) { int cellIdx = cells?cells[i]:i; - int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx); - const double* sdens = pvt_.surfaceDensities(pvtRegionIdx); + const double *sdens = surfaceDensity(cellIdx); for (int phase = 0; phase < np; ++phase) { rho[np*i + phase] = 0.0; for (int comp = 0; comp < np; ++comp) { @@ -391,8 +682,37 @@ namespace Opm /// \return Array of P density values. const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const { + const auto& pu = phaseUsage(); int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx); - return pvt_.surfaceDensities(pvtRegionIdx); + return &surfaceDensities_[pvtRegionIdx*pu.num_phases]; + } + + void BlackoilPropertiesFromDeck::initSurfaceDensities_(Opm::DeckConstPtr deck) + { + const auto& pu = phaseUsage(); + int np = pu.num_phases; + int numPvtRegions = 1; + if (deck->hasKeyword("TABDIMS")) { + const auto& tabdimsKeyword = deck->getKeyword("TABDIMS"); + numPvtRegions = tabdimsKeyword.getRecord(0).getItem("NTPVT").template get(0); + } + + const auto& densityKeyword = deck->getKeyword("DENSITY"); + + surfaceDensities_.resize(np*numPvtRegions); + for (int pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) { + if (pu.phase_used[BlackoilPhases::Aqua]) + surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Aqua]] = + densityKeyword.getRecord(pvtRegionIdx).getItem("WATER").getSIDouble(0); + + if (pu.phase_used[BlackoilPhases::Liquid]) + surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Liquid]] = + densityKeyword.getRecord(pvtRegionIdx).getItem("OIL").getSIDouble(0); + + if (pu.phase_used[BlackoilPhases::Vapour]) + surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Vapour]] = + densityKeyword.getRecord(pvtRegionIdx).getItem("GAS").getSIDouble(0); + } } /// \param[in] n Number of data points. diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index fb37a7f8..e5c18388 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -23,9 +23,11 @@ #include #include -#include #include #include +#include +#include +#include #include @@ -233,10 +235,19 @@ namespace Opm const double pcow, double & swat); - // return a reference to the "raw" PVT fluid object for a phase. - const PvtInterface& pvt(int phaseIdx) const + const OilPvtMultiplexer& oilPvt() const { - return pvt_.pvt(phaseIdx); + return oilPvt_; + } + + const GasPvtMultiplexer& gasPvt() const + { + return gasPvt_; + } + + const WaterPvtMultiplexer& waterPvt() const + { + return waterPvt_; } private: @@ -247,6 +258,38 @@ namespace Opm return pvtTableIdx[cellIdx]; } + void initSurfaceDensities_(Opm::DeckConstPtr deck); + + void compute_B_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B) const; + + void compute_dBdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* B, + double* dBdp) const; + + void compute_R_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R) const; + + void compute_dRdp_(const int n, + const double* p, + const double* T, + const double* z, + const int* cells, + double* R, + double* dRdp) const; + void init(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState, std::shared_ptr materialLawManager, @@ -265,10 +308,14 @@ namespace Opm bool init_rock); RockFromDeck rock_; + PhaseUsage phaseUsage_; std::vector cellPvtRegionIdx_; - BlackoilPvtProperties pvt_; + OilPvtMultiplexer oilPvt_; + GasPvtMultiplexer gasPvt_; + WaterPvtMultiplexer waterPvt_; std::shared_ptr materialLawManager_; std::shared_ptr satprops_; + std::vector surfaceDensities_; mutable std::vector B_; mutable std::vector dB_; mutable std::vector R_; diff --git a/opm/core/props/IncompPropertiesFromDeck.hpp b/opm/core/props/IncompPropertiesFromDeck.hpp index e71ecb78..4932ed80 100644 --- a/opm/core/props/IncompPropertiesFromDeck.hpp +++ b/opm/core/props/IncompPropertiesFromDeck.hpp @@ -23,8 +23,8 @@ #include #include #include -#include #include +#include #include struct UnstructuredGrid; diff --git a/opm/core/props/pvt/BlackoilPvtProperties.cpp b/opm/core/props/pvt/BlackoilPvtProperties.cpp deleted file mode 100644 index 490370f5..00000000 --- a/opm/core/props/pvt/BlackoilPvtProperties.cpp +++ /dev/null @@ -1,255 +0,0 @@ -/* - 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 "config.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -namespace Opm -{ - - BlackoilPvtProperties::BlackoilPvtProperties() - { - } - - void BlackoilPvtProperties::init(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, - int numSamples) - { - phase_usage_ = phaseUsageFromDeck(deck); - - // Surface densities. Accounting for different orders in eclipse and our code. - const auto& densityKeyword = deck->getKeyword("DENSITY"); - int numRegions = densityKeyword.size(); - - densities_.resize(numRegions); - for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - if (phase_usage_.phase_used[Liquid]) { - densities_[regionIdx][phase_usage_.phase_pos[Liquid]] - = densityKeyword.getRecord(regionIdx).getItem("OIL").getSIDouble(0); - } - if (phase_usage_.phase_used[Aqua]) { - densities_[regionIdx][phase_usage_.phase_pos[Aqua]] - = densityKeyword.getRecord(regionIdx).getItem("WATER").getSIDouble(0); - } - if (phase_usage_.phase_used[Vapour]) { - densities_[regionIdx][phase_usage_.phase_pos[Vapour]] - = densityKeyword.getRecord(regionIdx).getItem("GAS").getSIDouble(0); - } - } - - // Resize the property objects container - props_.resize(phase_usage_.num_phases); - - // Water PVT - if (phase_usage_.phase_used[Aqua]) { - // if water is used, we require the presence of the "PVTW" - // keyword for now... - std::shared_ptr pvtw(new PvtConstCompr); - pvtw->initFromWater(deck->getKeyword("PVTW")); - - props_[phase_usage_.phase_pos[Aqua]] = pvtw; - } - - { - auto tables = eclipseState->getTableManager(); - // Oil PVT - if (phase_usage_.phase_used[Liquid]) { - // for oil, we support the "PVDO", "PVTO" and "PVCDO" - // keywords... - const auto &pvdoTables = tables->getPvdoTables(); - const auto &pvtoTables = tables->getPvtoTables(); - if (pvdoTables.size() > 0) { - if (numSamples > 0) { - auto splinePvt = std::shared_ptr(new PvtDeadSpline); - splinePvt->initFromOil(pvdoTables, numSamples); - props_[phase_usage_.phase_pos[Liquid]] = splinePvt; - } else { - auto deadPvt = std::shared_ptr(new PvtDead); - deadPvt->initFromOil(pvdoTables); - props_[phase_usage_.phase_pos[Liquid]] = deadPvt; - } - } else if (pvtoTables.size() > 0) { - props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(pvtoTables)); - } else if (deck->hasKeyword("PVCDO")) { - std::shared_ptr pvcdo(new PvtConstCompr); - pvcdo->initFromOil(deck->getKeyword("PVCDO")); - - props_[phase_usage_.phase_pos[Liquid]] = pvcdo; - } else { - OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n"); - } - } - // Gas PVT - if (phase_usage_.phase_used[Vapour]) { - // gas can be specified using the "PVDG" or "PVTG" keywords... - const auto &pvdgTables = tables->getPvdgTables(); - const auto &pvtgTables = tables->getPvtgTables(); - if (pvdgTables.size() > 0) { - if (numSamples > 0) { - std::shared_ptr splinePvt(new PvtDeadSpline); - splinePvt->initFromGas(pvdgTables, numSamples); - - props_[phase_usage_.phase_pos[Vapour]] = splinePvt; - } else { - std::shared_ptr deadPvt(new PvtDead); - deadPvt->initFromGas(pvdgTables); - - props_[phase_usage_.phase_pos[Vapour]] = deadPvt; - } - } else if (pvtgTables.size() > 0) { - props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(pvtgTables)); - } else { - OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); - } - } - } - } - - const double* BlackoilPvtProperties::surfaceDensities(int regionIdx) const - { - return &densities_[regionIdx][0]; - } - - - PhaseUsage BlackoilPvtProperties::phaseUsage() const - { - return phase_usage_; - } - - int BlackoilPvtProperties::numPhases() const - { - return phase_usage_.num_phases; - } - - const int* BlackoilPvtProperties::phaseUsed() const - { - return phase_usage_.phase_used; - } - - const int* BlackoilPvtProperties::phasePosition() const - { - return phase_usage_.phase_pos; - } - - - void BlackoilPvtProperties::mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_mu) const - { - data1_.resize(n); - for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->mu(n, pvtTableIdx, p, T, z, &data1_[0]); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_mu[phase_usage_.num_phases*i + phase] = data1_[i]; - } - } - } - - void BlackoilPvtProperties::B(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B) const - { - data1_.resize(n); - for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->B(n, pvtTableIdx, p, T, z, &data1_[0]); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_B[phase_usage_.num_phases*i + phase] = data1_[i]; - } - } - } - - void BlackoilPvtProperties::dBdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B, - double* output_dBdp) const - { - data1_.resize(n); - data2_.resize(n); - for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->dBdp(n, pvtTableIdx, p, T, z, &data1_[0], &data2_[0]); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_B[phase_usage_.num_phases*i + phase] = data1_[i]; - output_dBdp[phase_usage_.num_phases*i + phase] = data2_[i]; - } - } - } - - - void BlackoilPvtProperties::R(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R) const - { - data1_.resize(n); - for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->R(n, pvtTableIdx, p, z, &data1_[0]); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_R[phase_usage_.num_phases*i + phase] = data1_[i]; - } - } - } - - void BlackoilPvtProperties::dRdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const - { - data1_.resize(n); - data2_.resize(n); - for (int phase = 0; phase < phase_usage_.num_phases; ++phase) { - props_[phase]->dRdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_R[phase_usage_.num_phases*i + phase] = data1_[i]; - output_dRdp[phase_usage_.num_phases*i + phase] = data2_[i]; - } - } - } -} // namespace Opm diff --git a/opm/core/props/pvt/BlackoilPvtProperties.hpp b/opm/core/props/pvt/BlackoilPvtProperties.hpp deleted file mode 100644 index c1778094..00000000 --- a/opm/core/props/pvt/BlackoilPvtProperties.hpp +++ /dev/null @@ -1,146 +0,0 @@ -/* - 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_BLACKOILPVTPROPERTIES_HEADER_INCLUDED -#define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED - -#include -#include - -#include -#include - -#include -#include -#include - -namespace Opm -{ - - /// Class collecting the pvt properties for all active phases. - /// For all the methods, the following apply: - /// - p and z are expected to be of size n and n*num_phases, respectively. - /// - pvtTableIdx specifies the PVT table to be used for each data - /// point and is thus expected to be an array of size n - /// - Output arrays shall be of size n*num_phases, and must be valid - /// before calling the method. - /// NOTE: The difference between this interface and the one defined - /// by PvtInterface is that this collects all phases' properties, - /// and therefore the output arrays are of size n*num_phases as opposed - /// to size n in PvtInterface. - class BlackoilPvtProperties : public BlackoilPhases - { - public: - /// Default constructor. - BlackoilPvtProperties(); - - /// Initialize from deck. - /// - /// \param deck An input deck from the opm-parser module. - void init(Opm::DeckConstPtr deck, - Opm::EclipseStateConstPtr eclipseState, - int samples); - - /// \return Object describing the active phases. - PhaseUsage phaseUsage() const; - - /// Number of active phases. - int numPhases() const; - - /// For each canonical phase, indicates if it is - /// active or not (boolean usage of int). - /// \return Array of size MaxNumPhases - const int* phaseUsed() const; - - /// Positions of canonical phases in arrays of phase - /// properties (saturations etc.). - /// \return Array of size MaxNumPhases - const int* phasePosition() const; - - /// Densities of stock components at surface conditions. - /// \return Array of size numPhases(). - const double* surfaceDensities(int regionIdx = 0) const; - - /// Viscosity as a function of p, T and z. - void mu(const int n, - const int *pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_mu) const; - - /// Formation volume factor as a function of p, T and z. - void B(const int n, - const int *pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B) const; - - /// Formation volume factor and p-derivative as functions of p, T and z. - void dBdp(const int n, - const int *pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B, - double* output_dBdp) const; - - /// Solution factor as a function of p and z. - void R(const int n, - const int *pvtTableIdx, - const double* p, - const double* z, - double* output_R) const; - - /// Solution factor and p-derivative as functions of p and z. - void dRdp(const int n, - const int *pvtTableIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const; - - // return a reference to the "raw" PVT fluid object for a phase. - const PvtInterface& pvt(int phaseIdx) const - { - return *props_[phaseIdx]; - } - - private: - // Disabling copying (just to avoid surprises, since we use shared_ptr). - BlackoilPvtProperties(const BlackoilPvtProperties&); - BlackoilPvtProperties& operator=(const BlackoilPvtProperties&); - - PhaseUsage phase_usage_; - - // The PVT properties. We need to store one object per PVT - // region per active fluid phase. - std::vector > props_; - std::vector > densities_; - - mutable std::vector data1_; - mutable std::vector data2_; - }; - -} - - - -#endif // OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED diff --git a/opm/core/props/pvt/PvtConstCompr.hpp b/opm/core/props/pvt/PvtConstCompr.hpp deleted file mode 100644 index 653f811f..00000000 --- a/opm/core/props/pvt/PvtConstCompr.hpp +++ /dev/null @@ -1,314 +0,0 @@ -/* - Copyright 2010, 2011, 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_PVTCONSTCOMPR_HEADER_INCLUDED -#define OPM_PVTCONSTCOMPR_HEADER_INCLUDED - -#include -#include -#include -#include -#include - -#include -#include - - -namespace Opm -{ - - /// Class for constant compressible phases (PVTW or PVCDO). - /// The PVT properties can either be given as a function of - /// pressure (p) and surface volume (z) or pressure (p) and gas - /// resolution factor (r). Also, since this class supports - /// multiple PVT regions, the concrete table to be used for each - /// data point needs to be specified via the pvtTableIdx argument - /// of the respective method. For all the virtual methods, the - /// following apply: pvtTableIdx, p, r and z are expected to be of - /// size n, size n, size n and n*num_phases, respectively. Output - /// arrays shall be of size n, and must be valid before calling - /// the method. - class PvtConstCompr : public PvtInterface - { - public: - PvtConstCompr() - {} - - void initFromWater(const Opm::DeckKeyword& pvtwKeyword) - { - int numRegions = pvtwKeyword.size(); - - ref_press_.resize(numRegions); - ref_B_.resize(numRegions); - comp_.resize(numRegions); - viscosity_.resize(numRegions); - visc_comp_.resize(numRegions); - - for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { - const auto& pvtwRecord = pvtwKeyword.getRecord(regionIdx); - - ref_press_[regionIdx] = pvtwRecord.getItem("P_REF").getSIDouble(0); - ref_B_[regionIdx] = pvtwRecord.getItem("WATER_VOL_FACTOR").getSIDouble(0); - comp_[regionIdx] = pvtwRecord.getItem("WATER_COMPRESSIBILITY").getSIDouble(0); - viscosity_[regionIdx] = pvtwRecord.getItem("WATER_VISCOSITY").getSIDouble(0); - visc_comp_[regionIdx] = pvtwRecord.getItem("WATER_VISCOSIBILITY").getSIDouble(0); - } - } - - void initFromOil(const DeckKeyword& pvcdoKeyword) - { - int numRegions = pvcdoKeyword.size(); - - ref_press_.resize(numRegions); - ref_B_.resize(numRegions); - comp_.resize(numRegions); - viscosity_.resize(numRegions); - visc_comp_.resize(numRegions); - - for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { - const auto& pvcdoRecord = pvcdoKeyword.getRecord(regionIdx); - - ref_press_[regionIdx] = pvcdoRecord.getItem("P_REF").getSIDouble(0); - ref_B_[regionIdx] = pvcdoRecord.getItem("OIL_VOL_FACTOR").getSIDouble(0); - comp_[regionIdx] = pvcdoRecord.getItem("OIL_COMPRESSIBILITY").getSIDouble(0); - viscosity_[regionIdx] = pvcdoRecord.getItem("OIL_VISCOSITY").getSIDouble(0); - visc_comp_[regionIdx] = pvcdoRecord.getItem("OIL_VISCOSIBILITY").getSIDouble(0); - } - } - - /*! - * \brief Create a PVT object with a given viscosity that - * assumes all fluid phases to be incompressible. - */ - explicit PvtConstCompr(double visc) - : ref_press_(1, 0.0), - ref_B_(1, 1.0), - comp_(1, 0.0), - viscosity_(1, visc), - visc_comp_(1, 0.0) - { - } - - virtual ~PvtConstCompr() - { - } - - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* /*z*/, - double* output_mu) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(pvtRegionIdx, i); - double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); - output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x); - } - } - - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(pvtRegionIdx, i); - double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); - double d = (1.0 + x + 0.5*x*x); - output_mu[i] = viscosity_[tableIdx]/d; - output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx]; - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - } - - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(pvtRegionIdx, i); - double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); - double d = (1.0 + x + 0.5*x*x); - output_mu[i] = viscosity_[tableIdx]/d; - output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx]; - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - } - - virtual void B(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* /*z*/, - double* output_B) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(pvtRegionIdx, i); - double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); - output_B[i] = ref_B_[tableIdx]/(1.0 + x + 0.5*x*x); - } - } - - virtual void dBdp(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* /*z*/, - double* output_B, - double* output_dBdp) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int tableIdx = getTableIndex_(pvtRegionIdx, i); - double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); - double d = (1.0 + x + 0.5*x*x); - output_B[i] = ref_B_[tableIdx]/d; - output_dBdp[i] = (-ref_B_[tableIdx]/(d*d))*(1 + x) * comp_[tableIdx]; - } - } - - virtual void b(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(pvtRegionIdx, i); - double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); - double d = (1.0 + x + 0.5*x*x); - - // b = 1/B = d/ref_B_B; - output_b[i] = d/ref_B_[tableIdx]; - output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx]; - } - - std::fill(output_dbdr, output_dbdr + n, 0.0); - } - - virtual void b(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - // Computing a polynomial approximation to the exponential. - int tableIdx = getTableIndex_(pvtRegionIdx, i); - double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]); - double d = (1.0 + x + 0.5*x*x); - - // b = 1/B = d/ref_B_[tableIdx]B; - output_b[i] = d/ref_B_[tableIdx]; - output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx]; - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - } - - virtual void rsSat(const int n, - const int* /*pvtRegionIdx*/, - const double* /*p*/, - double* output_rsSat, - double* output_drsSatdp) const - { - std::fill(output_rsSat, output_rsSat + n, 0.0); - std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); - } - - virtual void rvSat(const int n, - const int* /*pvtRegionIdx*/, - const double* /*p*/, - double* output_rvSat, - double* output_drvSatdp) const - { - std::fill(output_rvSat, output_rvSat + n, 0.0); - std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); - } - - virtual void R(const int n, - const int* /*pvtRegionIdx*/, - const double* /*p*/, - const double* /*z*/, - double* output_R) const - { - std::fill(output_R, output_R + n, 0.0); - } - - virtual void dRdp(const int n, - const int* /*pvtRegionIdx*/, - const double* /*p*/, - const double* /*z*/, - double* output_R, - double* output_dRdp) const - { - std::fill(output_R, output_R + n, 0.0); - std::fill(output_dRdp, output_dRdp + n, 0.0); - } - - private: - int getTableIndex_(const int* pvtTableIdx, int cellIdx) const - { - if (!pvtTableIdx) - return 0; - return pvtTableIdx[cellIdx]; - } - - // The PVT properties. We need to store one value per PVT - // region. - std::vector ref_press_; - std::vector ref_B_; - std::vector comp_; - std::vector viscosity_; - std::vector visc_comp_; - }; - -} - - -#endif // OPM_PVTCONSTCOMPR_HEADER_INCLUDED - diff --git a/opm/core/props/pvt/PvtDead.cpp b/opm/core/props/pvt/PvtDead.cpp deleted file mode 100644 index ca2750df..00000000 --- a/opm/core/props/pvt/PvtDead.cpp +++ /dev/null @@ -1,294 +0,0 @@ -/* - 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 "config.h" -#include -#include -#include -#include - -// Extra includes for debug dumping of tables. -// #include -// #include -// #include - -namespace Opm -{ - - //------------------------------------------------------------------------ - // Member functions - //------------------------------------------------------------------------- - /// Constructor - void PvtDead::initFromOil(const TableContainer& pvdoTables) - { - int numRegions = pvdoTables.size(); - - // resize the attributes of the object - b_.resize(numRegions); - viscosity_.resize(numRegions); - inverseBmu_.resize(numRegions); - - for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - const Opm::PvdoTable& pvdoTable = pvdoTables.getTable(regionIdx); - - // Copy data - const auto& press = pvdoTable.getColumn("P"); - const auto& b_var = pvdoTable.getColumn("BO"); - const auto& visc = pvdoTable.getColumn("MUO"); - - const int sz = b_var.size(); - std::vector inverseB(sz); - for (int i = 0; i < sz; ++i) { - inverseB[i] = 1.0 / b_var[i]; - } - - std::vector inverseBmu(sz); - for (int i = 0; i < sz; ++i) { - inverseBmu[i] = 1.0 / (b_var[i] * visc[i]); - } - - b_[regionIdx] = NonuniformTableLinear(press, inverseB); - viscosity_[regionIdx] = NonuniformTableLinear(press, visc); - inverseBmu_[regionIdx] = NonuniformTableLinear(press, inverseBmu); - } - } - - - void PvtDead::initFromGas(const TableContainer& pvdgTables) - { - int numRegions = pvdgTables.size(); - - // resize the attributes of the object - b_.resize(numRegions); - viscosity_.resize(numRegions); - inverseBmu_.resize(numRegions); - - for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - const Opm::PvdgTable& pvdgTable = pvdgTables.getTable(regionIdx); - - // Copy data - const auto& press = pvdgTable.getColumn("P"); - const auto& b = pvdgTable.getColumn("BG"); - const auto& visc = pvdgTable.getColumn("MUG"); - - const int sz = b.size(); - std::vector inverseB(sz); - for (int i = 0; i < sz; ++i) { - inverseB[i] = 1.0 / b[i]; - } - - std::vector inverseBmu(sz); - for (int i = 0; i < sz; ++i) { - inverseBmu[i] = 1.0 / (b[i] * visc[i]); - } - - b_[regionIdx] = NonuniformTableLinear(press, inverseB); - viscosity_[regionIdx] = NonuniformTableLinear(press, visc); - inverseBmu_[regionIdx] = NonuniformTableLinear(press, inverseBmu); - } - } - - // Destructor - PvtDead::~PvtDead() - { - } - - - - void PvtDead::mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*z*/, - double* output_mu) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - double tempInvB = b_[regionIdx](p[i]); - double tempInvBmu = inverseBmu_[regionIdx](p[i]); - output_mu[i] = tempInvB / tempInvBmu; - } - } - - void PvtDead::mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - double tempInvB = b_[regionIdx](p[i]); - double tempInvBmu = inverseBmu_[regionIdx](p[i]); - output_mu[i] = tempInvB / tempInvBmu; - output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) - - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - - } - - void PvtDead::mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - double tempInvB = b_[regionIdx](p[i]); - double tempInvBmu = inverseBmu_[regionIdx](p[i]); - output_mu[i] = tempInvB / tempInvBmu; - output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i]) - - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) - / (tempInvBmu * tempInvBmu); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - - } - - void PvtDead::B(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*z*/, - double* output_B) const - { -// #pragma omp parallel for - // B = 1/b - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - output_B[i] = 1.0/b_[regionIdx](p[i]); - } - } - - void PvtDead::dBdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* /*z*/, - double* output_B, - double* output_dBdp) const - { - B(n, pvtTableIdx, p, T, 0, output_B); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - double Bg = output_B[i]; - output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]); - } - } - - void PvtDead::b(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - - output_b[i] = b_[regionIdx](p[i]); - output_dbdp[i] = b_[regionIdx].derivative(p[i]); - - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - - } - - void PvtDead::b(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - - output_b[i] = b_[regionIdx](p[i]); - output_dbdp[i] = b_[regionIdx].derivative(p[i]); - - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - - } - - void PvtDead::rsSat(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - double* output_rsSat, - double* output_drsSatdp) const - { - std::fill(output_rsSat, output_rsSat + n, 0.0); - std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); - } - - void PvtDead::rvSat(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - double* output_rvSat, - double* output_drvSatdp) const - { - std::fill(output_rvSat, output_rvSat + n, 0.0); - std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); - } - - void PvtDead::R(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - const double* /*z*/, - double* output_R) const - { - std::fill(output_R, output_R + n, 0.0); - } - - void PvtDead::dRdp(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - const double* /*z*/, - double* output_R, - double* output_dRdp) const - { - std::fill(output_R, output_R + n, 0.0); - std::fill(output_dRdp, output_dRdp + n, 0.0); - } - -} diff --git a/opm/core/props/pvt/PvtDead.hpp b/opm/core/props/pvt/PvtDead.hpp deleted file mode 100644 index b238909a..00000000 --- a/opm/core/props/pvt/PvtDead.hpp +++ /dev/null @@ -1,171 +0,0 @@ -/* - 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_PVTDEAD_HEADER_INCLUDED -#define OPM_PVTDEAD_HEADER_INCLUDED - -#include -#include - -#include -#include - -#include - -namespace Opm -{ - - /// Class for immiscible dead oil and dry gas. - /// The PVT properties can either be given as a function of - /// pressure (p) and surface volume (z) or pressure (p) and gas - /// resolution factor (r). Also, since this class supports - /// multiple PVT regions, the concrete table to be used for each - /// data point needs to be specified via the pvtTableIdx argument - /// of the respective method. For all the virtual methods, the - /// following apply: pvtTableIdx, p, r and z are expected to be of - /// size n, size n, size n and n*num_phases, respectively. Output - /// arrays shall be of size n, and must be valid before calling - /// the method. - class PvtDead : public PvtInterface - { - public: - PvtDead() {}; - - void initFromOil(const TableContainer& pvdoTables); - void initFromGas(const TableContainer& pvdgTables); - virtual ~PvtDead(); - - /// Viscosity as a function of p, T and z. - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_mu) const; - - /// Viscosity and its p and r derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Viscosity as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Formation volume factor as a function of p, T and z. - virtual void B(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B) const; - - /// Formation volume factor and p-derivative as functions of p, T and z. - virtual void dBdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B, - double* output_dBdp) const; - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void b(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void b(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - - /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. - virtual void rsSat(const int n, - const int* pvtTableIdx, - const double* p, - double* output_rsSat, - double* output_drsSatdp) const; - - /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. - virtual void rvSat(const int n, - const int* pvtTableIdx, - const double* p, - double* output_rvSat, - double* output_drvSatdp) const; - - /// Solution factor as a function of p and z. - virtual void R(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R) const; - - /// Solution factor and p-derivative as functions of p and z. - virtual void dRdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const; - private: - int getTableIndex_(const int* pvtTableIdx, int cellIdx) const - { - if (!pvtTableIdx) - return 0; - return pvtTableIdx[cellIdx]; - } - - // PVT properties of dry gas or dead oil. We need to store one - // table per PVT region. - std::vector > b_; - std::vector > viscosity_; - std::vector > inverseBmu_; - }; -} - - -#endif // OPM_PVTDEAD_HEADER_INCLUDED diff --git a/opm/core/props/pvt/PvtDeadSpline.cpp b/opm/core/props/pvt/PvtDeadSpline.cpp deleted file mode 100644 index b7452c52..00000000 --- a/opm/core/props/pvt/PvtDeadSpline.cpp +++ /dev/null @@ -1,268 +0,0 @@ -/* - Copyright 2010, 2011, 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 "config.h" -#include -#include -#include -#include - - -#include - -// Extra includes for debug dumping of tables. -// #include -// #include -// #include - -namespace Opm -{ - //------------------------------------------------------------------------ - // Member functions - //------------------------------------------------------------------------- - - PvtDeadSpline::PvtDeadSpline() - {} - - void PvtDeadSpline::initFromOil(const TableContainer& pvdoTables, - int numSamples) - { - int numRegions = pvdoTables.size(); - - // resize the attributes of the object - b_.resize(numRegions); - viscosity_.resize(numRegions); - - for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - const Opm::PvdoTable& pvdoTable = pvdoTables.getTable(regionIdx); - - int numRows = pvdoTable.numRows(); - - // Copy data - std::vector press = pvdoTable.getColumn("P").vectorCopy(); - std::vector B_var = pvdoTable.getColumn("BO").vectorCopy(); - std::vector visc = pvdoTable.getColumn("MUO").vectorCopy(); - - std::vector B_inv(numRows); - for (int i = 0; i < numRows; ++i) { - B_inv[i] = 1.0 / B_var[i]; - } - - buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]); - buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]); - } - } - - void PvtDeadSpline::initFromGas(const TableContainer& pvdgTables, - int numSamples) - { - int numRegions = pvdgTables.size(); - - // resize the attributes of the object - b_.resize(numRegions); - viscosity_.resize(numRegions); - - for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - const Opm::PvdgTable& pvdgTable = pvdgTables.getTable(regionIdx); - - int numRows = pvdgTable.numRows(); - - // Copy data - std::vector press = pvdgTable.getColumn("P").vectorCopy(); - std::vector B = pvdgTable.getColumn("BG").vectorCopy(); - std::vector visc = pvdgTable.getColumn("MUG").vectorCopy(); - - std::vector B_inv(numRows); - for (int i = 0; i < numRows; ++i) { - B_inv[i] = 1.0 / B[i]; - } - - buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]); - buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]); - } - } - - // Destructor - PvtDeadSpline::~PvtDeadSpline() - { - } - - - - void PvtDeadSpline::mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*z*/, - double* output_mu) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = viscosity_[regionIdx](p[i]); - } - } - - void PvtDeadSpline::mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = viscosity_[regionIdx](p[i]); - output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - } - - void PvtDeadSpline::mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - // #pragma omp parallel for - - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - output_mu[i] = viscosity_[regionIdx](p[i]); - output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]); - } - std::fill(output_dmudr, output_dmudr + n, 0.0); - } - - void PvtDeadSpline::B(const int n, - const int* pvtTableIdx, - const double* p, - const double* /*T*/, - const double* /*z*/, - double* output_B) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - output_B[i] = 1.0/b_[regionIdx](p[i]); - } - } - - void PvtDeadSpline::dBdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* /*z*/, - double* output_B, - double* output_dBdp) const - { - B(n, pvtTableIdx, p, T, 0, output_B); -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - double Bg = output_B[i]; - output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]); - } - } - - void PvtDeadSpline::b(const int n, - const int* pvtTableIdx, - const double* p, - const double* /* T */, - const double* /*r*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - output_b[i] = b_[regionIdx](p[i]); - output_dbdp[i] = b_[regionIdx].derivative(p[i]); - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - } - - void PvtDeadSpline::b(const int n, - const int* pvtTableIdx, - const double* p, - const double* /* T */, - const double* /*r*/, - const PhasePresence* /*cond*/, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - int regionIdx = getTableIndex_(pvtTableIdx, i); - output_b[i] = b_[regionIdx](p[i]); - output_dbdp[i] = b_[regionIdx].derivative(p[i]); - } - std::fill(output_dbdr, output_dbdr + n, 0.0); - } - - void PvtDeadSpline::rsSat(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - double* output_rsSat, - double* output_drsSatdp) const - { - std::fill(output_rsSat, output_rsSat + n, 0.0); - std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); - } - - void PvtDeadSpline::rvSat(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - double* output_rvSat, - double* output_drvSatdp) const - { - std::fill(output_rvSat, output_rvSat + n, 0.0); - std::fill(output_drvSatdp, output_drvSatdp + n, 0.0); - } - - void PvtDeadSpline::R(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - const double* /*z*/, - double* output_R) const - { - std::fill(output_R, output_R + n, 0.0); - } - - void PvtDeadSpline::dRdp(const int n, - const int* /*pvtTableIdx*/, - const double* /*p*/, - const double* /*z*/, - double* output_R, - double* output_dRdp) const - { - std::fill(output_R, output_R + n, 0.0); - std::fill(output_dRdp, output_dRdp + n, 0.0); - } - -} diff --git a/opm/core/props/pvt/PvtDeadSpline.hpp b/opm/core/props/pvt/PvtDeadSpline.hpp deleted file mode 100644 index 9c087127..00000000 --- a/opm/core/props/pvt/PvtDeadSpline.hpp +++ /dev/null @@ -1,169 +0,0 @@ -/* - Copyright 2010, 2011, 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_PVTDEADSPLINE_HEADER_INCLUDED -#define OPM_PVTDEADSPLINE_HEADER_INCLUDED - -#include -#include - -#include -#include - -#include - -namespace Opm -{ - - /// Class for immiscible dead oil and dry gas. - /// The PVT properties can either be given as a function of pressure (p) and surface volume (z) - /// or pressure (p) and gas resolution factor (r). - /// For all the virtual methods, the following apply: p, r and z - /// are expected to be of size n, size n and n*num_phases, respectively. - /// Output arrays shall be of size n, and must be valid before - /// calling the method. - class PvtDeadSpline : public PvtInterface - { - public: - PvtDeadSpline(); - - void initFromOil(const TableContainer& pvdoTables, - int numSamples); - void initFromGas(const TableContainer& pvdgTables, - int numSamples); - - virtual ~PvtDeadSpline(); - - /// Viscosity as a function of p, T and z. - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_mu) const; - - /// Viscosity and its p and r derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Viscosity as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Formation volume factor as a function of p, T and z. - virtual void B(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B) const; - - /// Formation volume factor and p-derivative as functions of p, T and z. - virtual void dBdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B, - double* output_dBdp) const; - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void b(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void b(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. - virtual void rsSat(const int n, - const int* pvtTableIdx, - const double* p, - double* output_rsSat, - double* output_drsSatdp) const; - - /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. - virtual void rvSat(const int n, - const int* pvtTableIdx, - const double* p, - double* output_rvSat, - double* output_drvSatdp) const; - - /// Solution factor as a function of p and z. - virtual void R(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R) const; - - /// Solution factor and p-derivative as functions of p and z. - virtual void dRdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const; - private: - int getTableIndex_(const int* pvtTableIdx, int cellIdx) const - { - if (!pvtTableIdx) - return 0; - return pvtTableIdx[cellIdx]; - } - - // PVT properties of dry gas or dead oil. We need to store one - // table per PVT region. - std::vector > b_; - std::vector > viscosity_; - }; - -} - -#endif // OPM_PVTDEADSPLINE_HEADER_INCLUDED - diff --git a/opm/core/props/pvt/PvtInterface.cpp b/opm/core/props/pvt/PvtInterface.cpp deleted file mode 100644 index aa31a576..00000000 --- a/opm/core/props/pvt/PvtInterface.cpp +++ /dev/null @@ -1,67 +0,0 @@ -/* - Copyright 2010, 2011, 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 "config.h" -#include -#include - -#include - -namespace Opm -{ - - PvtInterface::PvtInterface() - : num_phases_(MaxNumPhases) - { - for (int i = 0; i < MaxNumPhases; ++i) { - phase_pos_[i] = i; - } - } - - PvtInterface::~PvtInterface() - { - } - - void PvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos) - { - num_phases_ = num_phases; - for (int i = 0; i < MaxNumPhases; ++i) { - phase_pos_[i] = phase_pos[i]; - } - } - - void extractPvtTableIndex(std::vector &pvtTableIdx, - Opm::EclipseStateConstPtr eclState, - size_t numCompressed, - const int *compressedToCartesianCellIdx) - { - //Get the PVTNUM data - const std::vector& pvtnumData = eclState->getIntGridProperty("PVTNUM")->getData(); - // Convert this into an array of compressed cells - // Eclipse uses Fortran-style indices which start at 1 - // instead of 0, we subtract 1. - pvtTableIdx.resize(numCompressed); - for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) { - size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx; - assert(cartesianCellIdx < pvtnumData.size()); - pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1; - } - } - -} // namespace Opm diff --git a/opm/core/props/pvt/PvtInterface.hpp b/opm/core/props/pvt/PvtInterface.hpp deleted file mode 100644 index 7c9c24cb..00000000 --- a/opm/core/props/pvt/PvtInterface.hpp +++ /dev/null @@ -1,187 +0,0 @@ -/* - Copyright 2010, 2011, 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_PVTINTERFACE_HEADER_INCLUDED -#define OPM_PVTINTERFACE_HEADER_INCLUDED - -#include - -#include -#include - - -namespace Opm -{ - - class PvtInterface : public BlackoilPhases - { - public: - PvtInterface(); - - virtual ~PvtInterface(); - - /// \param[in] num_phases The number of active phases. - /// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases - /// integers, giving the relative - /// positions of the three canonical - /// phases A, L, V in order to handle - /// arbitrary two-phase and three-phase situations. - void setPhaseConfiguration(const int num_phases, const int* phase_pos); - - /// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z) - /// or pressure (p), temperature (T) and gas resolution factor (r). - /// For all the virtual methods, the following apply: - /// - pvtRegionIdx is an array of size n and represents the - /// index of the PVT table which should be used to calculate - /// the output. NULL can also be passed and is interpreted - /// such that the first table should be used for the output - /// - p, r and z are expected to be of size n, size n and - /// n*num_phases, respectively. - /// - Output arrays shall be of size n, and must be valid before - /// calling the method. - - /// Viscosity as a function of p, T and z. - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* z, - double* output_mu) const = 0; - - /// Viscosity as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const = 0; - - /// Viscosity as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const = 0; - - /// Formation volume factor as a function of p, T and z. - virtual void B(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* z, - double* output_B) const = 0; - - /// Formation volume factor and p-derivative as functions of p and z. - virtual void dBdp(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* z, - double* output_B, - double* output_dBdp) const = 0; - - /// The inverse of the volume factor b = 1 / B as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void b(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - double* output_b, - double* output_dbdp, - double* output_dpdr) const = 0; - - /// The inverse of the volume factor b = 1 / B as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void b(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_b, - double* output_dbdp, - double* output_dpdr) const = 0; - - /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. - virtual void rsSat(const int n, - const int* pvtRegionIdx, - const double* p, - double* output_rsSat, - double* output_drsSatdp) const = 0; - - /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. - virtual void rvSat(const int n, - const int* pvtRegionIdx, - const double* p, - double* output_rvSat, - double* output_drvSatdp) const = 0; - - - /// Solution factor as a function of p and z. - virtual void R(const int n, - const int* pvtRegionIdx, - const double* p, - const double* z, - double* output_R) const = 0; - - /// Solution factor and p-derivative as functions of p and z. - virtual void dRdp(const int n, - const int* pvtRegionIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const = 0; - - - protected: - int num_phases_; - int phase_pos_[MaxNumPhases]; - }; - - /*! - * \brief Helper function to create an array containing the (C-Style) - * PVT table index for each compressed cell from an Eclipse deck. - * - * This function assumes that the degrees of freedom where PVT - * properties need to be calculated are grid cells. The main point - * of this function is to avoid code duplication because the - * Eclipse deck only contains Fortran-style PVT table indices - * which start at 1 instead of 0 and -- more relevantly -- it uses - * logically cartesian cell indices to specify the table index of - * a cell while the classes which use the PvtInterface - * implementations usually use compressed cells. - */ - void extractPvtTableIndex(std::vector& pvtTableIdx, - Opm::EclipseStateConstPtr eclState, - size_t numCompressed, - const int* compressedToCartesianIdx); - -} // namespace Opm - -#endif // OPM_PVTINTERFACE_HEADER_INCLUDED - diff --git a/opm/core/props/pvt/PvtLiveGas.cpp b/opm/core/props/pvt/PvtLiveGas.cpp deleted file mode 100644 index 51a22003..00000000 --- a/opm/core/props/pvt/PvtLiveGas.cpp +++ /dev/null @@ -1,613 +0,0 @@ -//=========================================================================== -// -// File: MiscibilityLiveGas.cpp -// -// Created: Wed Feb 10 09:21:53 2010 -// -// Author: Bjørn Spjelkavik -// -// Revision: $Id$ -// -//=========================================================================== -/* - Copyright 2010 SINTEF ICT, Applied Mathematics. - Copyright 2015 IRIS AS - - 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 "config.h" -#include -#include -#include -#include -#include - -#include -#include - -namespace Opm -{ - - using Opm::linearInterpolation; - using Opm::linearInterpolationDerivative; - - - //------------------------------------------------------------------------ - // Member functions - //------------------------------------------------------------------------- - PvtLiveGas::PvtLiveGas(const std::vector& pvtgTables) - { - int numTables = pvtgTables.size(); - saturated_gas_table_.resize(numTables); - undersat_gas_tables_.resize(numTables); - - for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) { - const Opm::PvtgTable& pvtgTable = pvtgTables[pvtTableIdx]; - const auto& saturatedPvtg = pvtgTable.getSaturatedTable( ); - // GAS, P - // saturated_gas_table_[pvtTableIdx].resize(4); - // Adding one extra line to PVTG to store 1./(Bg*mu_g) - saturated_gas_table_[pvtTableIdx].resize(5); - for (int k=0; k<5; ++k) { - saturated_gas_table_[pvtTableIdx][k].resize(saturatedPvtg.numRows()); - } - - for (size_t row=0; row < saturatedPvtg.numRows(); row++) { - saturated_gas_table_[pvtTableIdx][0][row] = saturatedPvtg.get("PG" , row); // Pressure - saturated_gas_table_[pvtTableIdx][1][row] = saturatedPvtg.get("BG" , row); // Bg - saturated_gas_table_[pvtTableIdx][2][row] = saturatedPvtg.get("MUG" , row); // mu_g - // 1/Bg will go in [3] - saturated_gas_table_[pvtTableIdx][4][row] = saturatedPvtg.get("RV" , row); // Rv - } - - - int sz = saturatedPvtg.numRows(); - undersat_gas_tables_[pvtTableIdx].resize(sz); - for (int i=0; i 1/Bg - for (int i=0; i 1) { - continue; - } - // Look ahead for next record containing undersaturated data - if (iNext < i) { - iNext = i+1; - while (iNext < sz && undersat_gas_tables_[pvtTableIdx][iNext][0].size() < 2) { - ++iNext; - } - if (iNext == sz) { - OPM_THROW(std::runtime_error,"Unable to complete undersaturated table."); - } - } - // Undersaturated data is added to the current record in the same way as for liveoil. - // It is unclear whether this expantion maintains the compressibility and viscosibility, - // but it seems like it reproduces eclipse results for spe3. - typedef std::vector > >::size_type sz_t; - auto& from_table = undersat_gas_tables_[pvtTableIdx][iNext]; - auto& to_table = undersat_gas_tables_[pvtTableIdx][i]; - enum {RV = 0, BG = 1, MUG = 2, BGMUG = 3}; - - for (sz_t j = 1; j < from_table[0].size(); ++j) { - double diffSolubility = from_table[RV][j] - from_table[RV][j-1]; - double solubility = to_table[RV].back() + diffSolubility; - to_table[RV].push_back(solubility); - double compr = (1.0/from_table[BG][j] - 1.0/from_table[BG][j-1]) - / (0.5*(1.0/from_table[BG][j] + 1.0/from_table[BG][j-1])); - double B_var = (1.0/to_table[BG].back()) * (1.0+0.5*compr) / (1.0-0.5*compr); - to_table[BG].push_back(1.0/B_var); - double visc = (from_table[MUG][j] - from_table[MUG][j-1]) - / (0.5*(from_table[MUG][j] + from_table[MUG][j-1])); - double mu_var = (to_table[MUG].back()) * (1.0+0.5*visc) / (1.0-0.5*visc); - to_table[MUG].push_back(mu_var); - - // A try to expolate the 1/BMu with the expolated mu and B - double inverseBMu = 1.0 / (B_var*mu_var); - to_table[BGMUG].push_back(inverseBMu); - } - } - - } - } - - // Destructor - PvtLiveGas::~PvtLiveGas() - { - } - - - void PvtLiveGas::mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* z, - double* output_mu) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - double inverseB = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 1, false); - double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false); - - output_mu[i] = inverseB / inverseBMu; - } - } - - /// Viscosity and its p and r derivatives as a function of p, T and r. - void PvtLiveGas::mu(const int /*n*/, - const int* /*pvtRegionIdx*/, - const double* /*p*/, - const double* /*T*/, - const double* /*r*/, - double* /*output_mu*/, - double* /*output_dmudp*/, - double* /*output_dmudr*/) const - { - OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented"); - } - - /// Viscosity and its p and r derivatives as a function of p, T and r. - void PvtLiveGas::mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* r, - const PhasePresence* cond, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const - { - for (int i = 0; i < n; ++i) { - const PhasePresence& cnd = cond[i]; - int tableIdx = getTableIndex_(pvtRegionIdx, i); - - double inverseBMu = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 0); - double inverseB = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0); - - output_mu[i] = inverseB / inverseBMu; - - double dinverseBmudp = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 1); - double dinverseBdp = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1); - - output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp) - / (inverseBMu * inverseBMu); - - double dinverseBmudr = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 2); - double dinverseBdr = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2); - - output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr) - / (inverseBMu * inverseBMu); - - } - - } - - - /// Formation volume factor as a function of p, T and z. - void PvtLiveGas::B(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* z, - double* output_B) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_B[i] = evalB(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i)); - } - - } - - - /// Formation volume factor and p-derivative as functions of p and z. - void PvtLiveGas::dBdp(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* z, - double* output_B, - double* output_dBdp) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - evalBDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_B[i], output_dBdp[i]); - } - } - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - void PvtLiveGas::b(const int /*n*/, - const int* /*pvtRegionIdx*/, - const double* /*p*/, - const double* /*T*/, - const double* /*r*/, - double* /*output_b*/, - double* /*output_dbdp*/, - double* /*output_dbdr*/) const - - { - OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented"); - } - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - void PvtLiveGas::b(const int n, - const int* pvtRegionIdx, - const double* p, - const double* /*T*/, - const double* r, - const PhasePresence* cond, - double* output_b, - double* output_dbdp, - double* output_dbdr) const - - { - // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - const PhasePresence& cnd = cond[i]; - - int tableIdx = getTableIndex_(pvtRegionIdx, i); - output_b[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0); - output_dbdp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1); - output_dbdr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2); - - } - } - - /// Gas resolution and its derivatives at bublepoint as a function of p. - void PvtLiveGas::rvSat(const int n, - const int* pvtRegionIdx, - const double* p, - double* output_rvSat, - double* output_drvSatdp) const - { - for (int i = 0; i < n; ++i) { - int pvtTableIdx = getTableIndex_(pvtRegionIdx, i); - output_rvSat[i] = linearInterpolation(saturated_gas_table_[pvtTableIdx][0], - saturated_gas_table_[pvtTableIdx][4],p[i]); - output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0], - saturated_gas_table_[pvtTableIdx][4],p[i]); - - } - } - - void PvtLiveGas::rsSat(const int n, - const int* /*pvtRegionIdx*/, - const double* /*p*/, - double* output_rsSat, - double* output_drsSatdp) const - { - std::fill(output_rsSat, output_rsSat + n, 0.0); - std::fill(output_drsSatdp, output_drsSatdp + n, 0.0); - } - - /// Solution factor as a function of p and z. - void PvtLiveGas::R(const int n, - const int* pvtRegionIdx, - const double* p, - const double* z, - double* output_R) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - output_R[i] = evalR(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i)); - } - - } - - - /// Solution factor and p-derivative as functions of p and z. - void PvtLiveGas::dRdp(const int n, - const int* pvtRegionIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const - { -// #pragma omp parallel for - for (int i = 0; i < n; ++i) { - evalRDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_R[i], output_dRdp[i]); - } - } - - - // ---- Private methods ---- - - double PvtLiveGas::evalB(const double press, const double* surfvol, int pvtTableIdx) const - { - if (surfvol[phase_pos_[Vapour]] == 0.0) { - // To handle no-gas case. - return 1.0; - } - return 1.0/miscible_gas(press, surfvol, pvtTableIdx, 1, false); - } - - void PvtLiveGas::evalBDeriv(const double press, const double* surfvol, int pvtTableIdx, - double& Bval, double& dBdpval) const - { - if (surfvol[phase_pos_[Vapour]] == 0.0) { - // To handle no-gas case. - Bval = 1.0; - dBdpval = 0.0; - return; - } - Bval = evalB(press, surfvol, pvtTableIdx); - dBdpval = -Bval*Bval*miscible_gas(press, surfvol, pvtTableIdx, 1, true); - } - - double PvtLiveGas::evalR(const double press, const double* surfvol, int pvtTableIdx) const - { - if (surfvol[phase_pos_[Liquid]] == 0.0) { - // To handle no-gas case. - return 0.0; - } - double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0], - saturated_gas_table_[pvtTableIdx][4], press); - double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]]; - if (satR < maxR ) { - // Saturated case - return satR; - } else { - // Undersaturated case - return maxR; - } - } - - void PvtLiveGas::evalRDeriv(const double press, const double* surfvol, int pvtTableIdx, - double& Rval, double& dRdpval) const - { - if (surfvol[phase_pos_[Liquid]] == 0.0) { - // To handle no-gas case. - Rval = 0.0; - dRdpval = 0.0; - return; - } - double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0], - saturated_gas_table_[pvtTableIdx][4], press); - double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]]; - if (satR < maxR ) { - // Saturated case - Rval = satR; - dRdpval = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0], - saturated_gas_table_[pvtTableIdx][4], - press); - } else { - // Undersaturated case - Rval = maxR; - dRdpval = 0.0; - } - } - - double PvtLiveGas::miscible_gas(const double press, - const double* surfvol, - const int pvtTableIdx, - const int item, - const bool deriv) const - { - const std::vector > &saturatedGasTable = - saturated_gas_table_[pvtTableIdx]; - const std::vector > > &undersatGasTables = - undersat_gas_tables_[pvtTableIdx]; - - int section; - double Rval = linearInterpolation(saturatedGasTable[0], - saturatedGasTable[4], press, - section); - double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]]; - if (deriv) { - if (Rval < maxR ) { // Saturated case - return linearInterpolationDerivative(saturatedGasTable[0], - saturatedGasTable[item], - press); - } else { // Undersaturated case - int is = section; - if (undersatGasTables[is][0].size() < 2) { - double val = (saturatedGasTable[item][is+1] - - saturatedGasTable[item][is]) / - (saturatedGasTable[0][is+1] - - saturatedGasTable[0][is]); - return val; - } - double val1 = - linearInterpolation(undersatGasTables[is][0], - undersatGasTables[is][item], - maxR); - double val2 = - linearInterpolation(undersatGasTables[is+1][0], - undersatGasTables[is+1][item], - maxR); - double val = (val2 - val1)/ - (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]); - return val; - } - } else { - if (Rval < maxR ) { // Saturated case - return linearInterpolation(saturatedGasTable[0], - saturatedGasTable[item], - press); - } else { // Undersaturated case - int is = section; - // Extrapolate from first table section - if (is == 0 && press < saturatedGasTable[0][0]) { - return linearInterpolation(undersatGasTables[0][0], - undersatGasTables[0][item], - maxR); - } - - // Extrapolate from last table section - int ltp = saturatedGasTable[0].size() - 1; - if (is+1 == ltp && press > saturatedGasTable[0][ltp]) { - return linearInterpolation(undersatGasTables[ltp][0], - undersatGasTables[ltp][item], - maxR); - } - - // Interpolate between table sections - double w = (press - saturatedGasTable[0][is]) / - (saturatedGasTable[0][is+1] - - saturatedGasTable[0][is]); - if (undersatGasTables[is][0].size() < 2) { - double val = saturatedGasTable[item][is] + - w*(saturatedGasTable[item][is+1] - - saturatedGasTable[item][is]); - return val; - } - double val1 = - linearInterpolation(undersatGasTables[is][0], - undersatGasTables[is][item], - maxR); - double val2 = - linearInterpolation(undersatGasTables[is+1][0], - undersatGasTables[is+1][item], - maxR); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - double PvtLiveGas::miscible_gas(const double press, - const double r, - const PhasePresence& cond, - const int pvtTableIdx, - const int item, - const int deriv) const - { - const std::vector > &saturatedGasTable = - saturated_gas_table_[pvtTableIdx]; - const std::vector > > &undersatGasTables = - undersat_gas_tables_[pvtTableIdx]; - - const bool isSat = cond.hasFreeOil(); - - // Derivative w.r.t p - if (deriv == 1) { - if (isSat) { // Saturated case - return linearInterpolationDerivative(saturatedGasTable[0], - saturatedGasTable[item], - press); - } else { // Undersaturated case - int is = tableIndex(saturatedGasTable[0], press); - if (undersatGasTables[is][0].size() < 2) { - double val = (saturatedGasTable[item][is+1] - - saturatedGasTable[item][is]) / - (saturatedGasTable[0][is+1] - - saturatedGasTable[0][is]); - return val; - } - double val1 = - linearInterpolation(undersatGasTables[is][0], - undersatGasTables[is][item], - r); - double val2 = - linearInterpolation(undersatGasTables[is+1][0], - undersatGasTables[is+1][item], - r); - double val = (val2 - val1)/ - (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]); - return val; - } - } else if (deriv == 2){ - if (isSat) { - return 0; - } else { - int is = tableIndex(saturatedGasTable[0], press); - double w = (press - saturatedGasTable[0][is]) / - (saturatedGasTable[0][is+1] - saturatedGasTable[0][is]); - assert(undersatGasTables[is][0].size() >= 2); - assert(undersatGasTables[is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersatGasTables[is][0], - undersatGasTables[is][item], - r); - double val2 = - linearInterpolationDerivative(undersatGasTables[is+1][0], - undersatGasTables[is+1][item], - r); - - double val = val1 + w * (val2 - val1); - return val; - - } - } else { - if (isSat) { // Saturated case - return linearInterpolation(saturatedGasTable[0], - saturatedGasTable[item], - press); - } else { // Undersaturated case - int is = tableIndex(saturatedGasTable[0], press); - - // Extrapolate from last table section - //int ltp = saturatedGasTable[0].size() - 1; - //if (is+1 == ltp && press > saturatedGasTable[0][ltp]) { - // return linearInterpolation(undersatGasTables[ltp][0], - // undersatGasTables[ltp][item], - // r); - //} - - // Interpolate between table sections - double w = (press - saturatedGasTable[0][is]) / - (saturatedGasTable[0][is+1] - - saturatedGasTable[0][is]); - if (undersatGasTables[is][0].size() < 2) { - double val = saturatedGasTable[item][is] + - w*(saturatedGasTable[item][is+1] - - saturatedGasTable[item][is]); - return val; - } - double val1 = - linearInterpolation(undersatGasTables[is][0], - undersatGasTables[is][item], - r); - double val2 = - linearInterpolation(undersatGasTables[is+1][0], - undersatGasTables[is+1][item], - r); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - - - -} // namespace Opm diff --git a/opm/core/props/pvt/PvtLiveGas.hpp b/opm/core/props/pvt/PvtLiveGas.hpp deleted file mode 100644 index 93833de9..00000000 --- a/opm/core/props/pvt/PvtLiveGas.hpp +++ /dev/null @@ -1,182 +0,0 @@ -/* - Copyright 2010, 2011, 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_PVTLIVEGAS_HEADER_INCLUDED -#define OPM_PVTLIVEGAS_HEADER_INCLUDED - -#include - -#include - -#include - -namespace Opm -{ - class PvtgTable; - - /// Class for miscible wet gas (with vaporized oil in vapour phase). - /// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z) - /// or pressure (p), temperature (T) and gas resolution factor (r). - /// For all the virtual methods, the following apply: p, r and z - /// are expected to be of size n, size n and n*num_phases, respectively. - /// Output arrays shall be of size n, and must be valid before - /// calling the method. - class PvtLiveGas : public PvtInterface - { - public: - PvtLiveGas(const std::vector& pvtgTables); - virtual ~PvtLiveGas(); - - /// Viscosity as a function of p, T and z. - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* z, - double* output_mu) const; - - /// Viscosity and its p and r derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Viscosity as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void mu(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Formation volume factor as a function of p, T and z. - virtual void B(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* z, - double* output_B) const; - - /// Formation volume factor and p-derivative as functions of p, T and z. - virtual void dBdp(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* z, - double* output_B, - double* output_dBdp) const; - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void b(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - /// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void b(const int n, - const int* pvtRegionIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - - - /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p, T. - virtual void rsSat(const int n, - const int* pvtRegionIdx, - const double* p, - double* output_rsSat, - double* output_drsSatdp) const; - - /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. - virtual void rvSat(const int n, - const int* pvtRegionIdx, - const double* p, - double* output_rvSat, - double* output_drvSatdp) const; - - /// Solution factor as a function of p and z. - virtual void R(const int n, - const int* pvtRegionIdx, - const double* p, - const double* z, - double* output_R) const; - - /// Solution factor and p-derivative as functions of p and z. - virtual void dRdp(const int n, - const int* pvtRegionIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const; - - protected: - int getTableIndex_(const int* pvtTableIdx, int cellIdx) const - { - if (!pvtTableIdx) - return 0; - return pvtTableIdx[cellIdx]; - } - - double evalB(double press, const double* surfvol, int pvtTableIdx) const; - void evalBDeriv(double press, const double* surfvol, int pvtTableIdx, double& B, double& dBdp) const; - double evalR(double press, const double* surfvol, int pvtTableIdx) const; - void evalRDeriv(double press, const double* surfvol, int pvtTableIdx, double& R, double& dRdp) const; - - // item: 1=>B 2=>mu; - double miscible_gas(const double press, - const double* surfvol, - const int pvtTableIdx, - const int item, - const bool deriv = false) const; - double miscible_gas(const double press, - const double r, - const PhasePresence& cond, - const int pvtTableIdx, - const int item, - const int deriv = 0) const; - // PVT properties of wet gas (with vaporised oil). We need to - // store one table per PVT region. - std::vector< std::vector > > saturated_gas_table_; - std::vector< std::vector > > > undersat_gas_tables_; - }; - -} - -#endif // OPM_PVTLIVEGAS_HEADER_INCLUDED - diff --git a/opm/core/props/pvt/PvtLiveOil.cpp b/opm/core/props/pvt/PvtLiveOil.cpp deleted file mode 100644 index cfa25d43..00000000 --- a/opm/core/props/pvt/PvtLiveOil.cpp +++ /dev/null @@ -1,653 +0,0 @@ -/* - Copyright 2010, 2011, 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 "config.h" -#include -#include -#include -#include -#include - -#include - -namespace Opm -{ - - using Opm::linearInterpolation; - using Opm::linearInterpolationDerivative; - using Opm::tableIndex; - - //------------------------------------------------------------------------ - // Member functions - //------------------------------------------------------------------------- - PvtLiveOil::PvtLiveOil(const std::vector& pvtoTables) - { - int numTables = pvtoTables.size(); - saturated_oil_table_.resize(numTables); - undersat_oil_tables_.resize(numTables); - - for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) { - const Opm::PvtoTable& pvtoTable = pvtoTables[pvtTableIdx]; - - const auto& saturatedPvto = pvtoTable.getSaturatedTable(); - - // OIL, PVTO - // saturated_oil_table_[pvtTableIdx].resize(4); - // adding a extra colummn to the PVTO to store 1/(B*mu) - saturated_oil_table_[pvtTableIdx].resize(5); - const int sz = saturatedPvto.numRows(); - // for (int k=0; k<4; ++k) { - for (int k=0; k<5; ++k) { - saturated_oil_table_[pvtTableIdx][k].resize(sz); - } - for (int i=0; i 1) { - continue; - } - // Look ahead for next record containing undersaturated data - if (iNext < i) { - iNext = i+1; - while (iNext > >::size_type sz_t; - for (sz_t j=1; j= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } else { - if (Rval < maxR ) { // Saturated case - return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][item], - press); - } else { // Undersaturated case - // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR); - double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) / - (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - double PvtLiveOil::miscible_oil(const double press, - const double r, - const int pvtTableIdx, - const int item, - const int deriv) const - { - int section; - double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][4], - press, section); - // derivative with respect to frist component (pressure) - if (deriv == 1) { - if (Rval <= r ) { // Saturated case - return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][item], - press); - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); - double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / - (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - // derivative with respect to second component (r) - } else if (deriv == 2) { - if (Rval <= r ) { // Saturated case - return 0; - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - - double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]); - return val; - } - - - } else { - if (Rval <= r ) { // Saturated case - return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][item], - press); - } else { // Undersaturated case - // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); - double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / - (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - - double PvtLiveOil::miscible_oil(const double press, - const double r, - const PhasePresence& cond, - const int pvtTableIdx, - const int item, - const int deriv) const - { - const bool isSat = cond.hasFreeGas(); - - // derivative with respect to frist component (pressure) - if (deriv == 1) { - if (isSat) { // Saturated case - return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][item], - press); - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); - double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / - (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - // derivative with respect to second component (r) - } else if (deriv == 2) { - if (isSat) { // Saturated case - return 0; - } else { // Undersaturated case - int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - - double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]); - return val; - } - - - } - // no derivative - else { - if (isSat) { // Saturated case - return linearInterpolation(saturated_oil_table_[pvtTableIdx][0], - saturated_oil_table_[pvtTableIdx][item], - press); - } else { // Undersaturated case - // Interpolate between table sections - int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r); - double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) / - (saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]); - assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2); - assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2); - double val1 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0], - undersat_oil_tables_[pvtTableIdx][is][item], - press); - double val2 = - linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0], - undersat_oil_tables_[pvtTableIdx][is+1][item], - press); - double val = val1 + w*(val2 - val1); - return val; - } - } - } - -} // namespace Opm diff --git a/opm/core/props/pvt/PvtLiveOil.hpp b/opm/core/props/pvt/PvtLiveOil.hpp deleted file mode 100644 index 24ede096..00000000 --- a/opm/core/props/pvt/PvtLiveOil.hpp +++ /dev/null @@ -1,190 +0,0 @@ -/* - Copyright 2010, 2011, 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_PVTLIVEOIL_HEADER_INCLUDED -#define OPM_PVTLIVEOIL_HEADER_INCLUDED - -#include - -#include -#include - -#include - -namespace Opm -{ - - class PvtoTable; - - /// Class for miscible live oil (with dissolved gas in liquid phase). - /// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z) - /// or pressure (p), temperature (T) and gas resolution factor (r). - /// For all the virtual methods, the following apply: p, r and z - /// are expected to be of size n, size n and n*num_phases, respectively. - /// Output arrays shall be of size n, and must be valid before - /// calling the method. - class PvtLiveOil : public PvtInterface - { - public: - PvtLiveOil(const std::vector& pvtoTables); - virtual ~PvtLiveOil(); - - /// Viscosity as a function of p, T and z. - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_mu) const; - - /// Viscosity and its p and r derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Viscosity as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void mu(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_mu, - double* output_dmudp, - double* output_dmudr) const; - - /// Formation volume factor as a function of p, T and z. - virtual void B(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B) const; - - /// Formation volume factor and p-derivative as functions of p, T and z. - virtual void dBdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* z, - double* output_B, - double* output_dBdp) const; - - /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p, T and r. - /// The fluid is considered saturated if r >= rsSat(p). - virtual void b(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - /// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p, T and r. - /// State condition determined by 'cond'. - virtual void b(const int n, - const int* pvtTableIdx, - const double* p, - const double* T, - const double* r, - const PhasePresence* cond, - double* output_b, - double* output_dbdp, - double* output_dbdr) const; - - /// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p. - virtual void rsSat(const int n, - const int* pvtTableIdx, - const double* p, - double* output_rsSat, - double* output_drsSatdp) const; - - /// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p. - virtual void rvSat(const int n, - const int* pvtTableIdx, - const double* p, - double* output_rvSat, - double* output_drvSatdp) const; - - /// Solution factor as a function of p and z. - virtual void R(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R) const; - - /// Solution factor and p-derivative as functions of p and z. - virtual void dRdp(const int n, - const int* pvtTableIdx, - const double* p, - const double* z, - double* output_R, - double* output_dRdp) const; - - private: - int getTableIndex_(const int* pvtTableIdx, int cellIdx) const - { - if (!pvtTableIdx) - return 0; - return pvtTableIdx[cellIdx]; - } - - double evalB(size_t pvtTableIdx, double press, const double* surfvol) const; - void evalBDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& B, double& dBdp) const; - double evalR(size_t pvtTableIdx, double press, const double* surfvol) const; - void evalRDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& R, double& dRdp) const; - - // item: 1=>1/B 2=>mu; - double miscible_oil(const double press, - const double* surfvol, - const int pvtTableIdx, - const int item, - const bool deriv = false) const; - - double miscible_oil(const double press, - const double r, - const int pvtTableIdx, - const int item, - const int deriv = 0) const; - - double miscible_oil(const double press, - const double r, - const PhasePresence& cond, - const int pvtTableIdx, - const int item, - const int deriv = 0) const; - - // PVT properties of live oil (with dissolved gas). We need to - // store one table per PVT region. - std::vector > > saturated_oil_table_; - std::vector > > > undersat_oil_tables_; - }; - -} - -#endif // OPM_PVTLIVEOIL_HEADER_INCLUDED - diff --git a/opm/core/utility/extractPvtTableIndex.cpp b/opm/core/utility/extractPvtTableIndex.cpp new file mode 100644 index 00000000..ec1d4dc9 --- /dev/null +++ b/opm/core/utility/extractPvtTableIndex.cpp @@ -0,0 +1,46 @@ +/* + 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 "extractPvtTableIndex.hpp" + +#include +#include + +#include + +namespace Opm { + +void extractPvtTableIndex(std::vector &pvtTableIdx, + Opm::EclipseStateConstPtr eclState, + size_t numCompressed, + const int *compressedToCartesianCellIdx) +{ + //Get the PVTNUM data + const std::vector& pvtnumData = eclState->getIntGridProperty("PVTNUM")->getData(); + // Convert this into an array of compressed cells + // Eclipse uses Fortran-style indices which start at 1 + // instead of 0, we subtract 1. + pvtTableIdx.resize(numCompressed); + for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) { + size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx; + assert(cartesianCellIdx < pvtnumData.size()); + pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1; + } +} + +} diff --git a/opm/core/utility/extractPvtTableIndex.hpp b/opm/core/utility/extractPvtTableIndex.hpp new file mode 100644 index 00000000..550e258f --- /dev/null +++ b/opm/core/utility/extractPvtTableIndex.hpp @@ -0,0 +1,33 @@ +/* + 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_EXTRACT_PVT_TABLE_INDEX_HPP +#define OPM_EXTRACT_PVT_TABLE_INDEX_HPP + +#include + +#include + +namespace Opm { + +void extractPvtTableIndex(std::vector &pvtTableIdx, + Opm::EclipseStateConstPtr eclState, + size_t numCompressed, + const int *compressedToCartesianCellIdx); + +} + +#endif // OPM_EXTRACT_PVT_TABLE_INDEX_HPP diff --git a/opm/core/utility/thresholdPressures.hpp b/opm/core/utility/thresholdPressures.hpp index 5a2be77a..31a91658 100644 --- a/opm/core/utility/thresholdPressures.hpp +++ b/opm/core/utility/thresholdPressures.hpp @@ -128,8 +128,6 @@ void computeMaxDp(std::map, double>& maxDp, } // calculate the initial fluid densities for the gravity correction. - std::vector b(numCells); - std::vector> rho(numPhases); for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { rho[phaseIdx].resize(numCells); @@ -176,69 +174,55 @@ void computeMaxDp(std::map, double>& maxDp, } } - // calculate the inverse formation volume factors for the active phases and each cell + // calculate the densities of the active phases for each cell if (pu.phase_used[BlackoilPhases::Aqua]) { - std::vector dummy(numCells*BlackoilPhases::MaxNumPhases); const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; - const PvtInterface& pvtw = props.pvt(wpos); - pvtw.b(numCells, - pvtRegion.data(), - phasePressure[wpos].data(), - initialState.temperature().data(), - initialState.gasoilratio().data(), - cond.data(), - b.data(), - dummy.data(), - dummy.data()); + const auto& pvtw = props.waterPvt(); for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { - rho[wpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][wpos]*b[cellIdx]; + int pvtRegionIdx = pvtRegion[cellIdx]; + + double T = initialState.temperature()[cellIdx]; + double p = initialState.pressure()[cellIdx]; + double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p); + + rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b; } } if (pu.phase_used[BlackoilPhases::Liquid]) { - std::vector dummy(numCells*BlackoilPhases::MaxNumPhases); const int opos = pu.phase_pos[BlackoilPhases::Liquid]; - const PvtInterface& pvto = props.pvt(opos); - pvto.b(numCells, - pvtRegion.data(), - phasePressure[opos].data(), - initialState.temperature().data(), - initialState.gasoilratio().data(), - cond.data(), - b.data(), - dummy.data(), - dummy.data()); + const auto& pvto = props.oilPvt(); for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { - rho[opos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][opos]*b[cellIdx]; + int pvtRegionIdx = pvtRegion[cellIdx]; + double T = initialState.temperature()[cellIdx]; + double p = initialState.pressure()[cellIdx]; + double Rs = initialState.gasoilratio()[cellIdx]; + double b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); + + rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b; if (pu.phase_used[BlackoilPhases::Vapour]) { int gpos = pu.phase_pos[BlackoilPhases::Vapour]; - rho[opos][cellIdx] += - surfaceDensity[pvtRegion[cellIdx]][gpos]*initialState.gasoilratio()[cellIdx]*b[cellIdx]; + rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b; } } } if (pu.phase_used[BlackoilPhases::Vapour]) { - std::vector dummy(numCells*BlackoilPhases::MaxNumPhases); - const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; - const PvtInterface& pvtg = props.pvt(gpos); - pvtg.b(numCells, - pvtRegion.data(), - phasePressure[gpos].data(), - initialState.temperature().data(), - initialState.rv().data(), - cond.data(), - b.data(), - dummy.data(), - dummy.data()); + const int gpos = pu.phase_pos[BlackoilPhases::Liquid]; + const auto& pvtg = props.gasPvt(); for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { - rho[gpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][gpos]*b[cellIdx]; + int pvtRegionIdx = pvtRegion[cellIdx]; + double T = initialState.temperature()[cellIdx]; + double p = initialState.pressure()[cellIdx]; + double Rv = initialState.rv()[cellIdx]; + double b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv); + + rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b; if (pu.phase_used[BlackoilPhases::Liquid]) { - const int opos = pu.phase_pos[BlackoilPhases::Liquid]; - rho[gpos][cellIdx] += - surfaceDensity[pvtRegion[cellIdx]][opos]*initialState.rv()[cellIdx]*b[cellIdx]; + int opos = pu.phase_pos[BlackoilPhases::Liquid]; + rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b; } } } diff --git a/tests/norne_pvt.data b/tests/norne_pvt.data index 10de7806..bdb7af92 100644 --- a/tests/norne_pvt.data +++ b/tests/norne_pvt.data @@ -5,11 +5,22 @@ -- Copyright (C) 2015 Statoil +RUNSPEC TABDIMS --ntsfun ntpvt nssfun nppvt ntfip nrpvt ntendp 2 2 33 60 16 60 / +DIMENS +1 1 1 / + +GRID + +PROPS + +DENSITY + 859.5 1033.0 0.854 / Justert 22/7 + 860.04 1033.0 0.853 / Justert 22/7 PVTG diff --git a/tests/test_blackoilfluid.cpp b/tests/test_blackoilfluid.cpp deleted file mode 100644 index 862872ad..00000000 --- a/tests/test_blackoilfluid.cpp +++ /dev/null @@ -1,388 +0,0 @@ -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#if HAVE_DYNAMIC_BOOST_TEST -#define BOOST_TEST_DYN_LINK -#endif -#define NVERBOSE // to suppress our messages when throwing -#define BOOST_TEST_MODULE BlackoilFluidTest -#define BOOST_TEST_MAIN -#include -#include -#include -#include -#include -#include - -using namespace Opm; -using namespace std; - - -std::vector > getProps(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclipseState, PhaseUsage phase_usage_){ - enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 }; - int samples = 0; - std::shared_ptr tables = eclipseState->getTableManager(); - std::vector > props_; - // Set the properties. - props_.resize(phase_usage_.num_phases); - - // Water PVT - if (phase_usage_.phase_used[Aqua]) { - if (deck->hasKeyword("PVTW")) { - std::shared_ptr pvtw(new PvtConstCompr); - pvtw->initFromWater(deck->getKeyword("PVTW")); - - props_[phase_usage_.phase_pos[Aqua]] = pvtw; - } else { - // Eclipse 100 default. - props_[phase_usage_.phase_pos[Aqua]].reset(new PvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise)); - } - } - - // Oil PVT - if (phase_usage_.phase_used[Liquid]) { - const auto& pvdoTables = tables->getPvdoTables(); - const auto& pvtoTables = tables->getPvtoTables(); - if (!pvdoTables.empty()) { - if (samples > 0) { - std::shared_ptr splinePvt(new PvtDeadSpline); - splinePvt->initFromOil(pvdoTables, samples); - props_[phase_usage_.phase_pos[Liquid]] = splinePvt; - } else { - std::shared_ptr deadPvt(new PvtDead); - deadPvt->initFromOil(pvdoTables); - props_[phase_usage_.phase_pos[Liquid]] = deadPvt; - } - } else if (!pvtoTables.empty()) { - props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(pvtoTables)); - } else if (deck->hasKeyword("PVCDO")) { - std::shared_ptr pvcdo(new PvtConstCompr); - pvcdo->initFromOil(deck->getKeyword("PVCDO")); - - props_[phase_usage_.phase_pos[Liquid]] = pvcdo; - } else { - OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n"); - } - } - // Gas PVT - if (phase_usage_.phase_used[Vapour]) { - const auto& pvdgTables = tables->getPvdgTables(); - const auto& pvtgTables = tables->getPvtgTables(); - if (!pvdgTables.empty()) { - if (samples > 0) { - std::shared_ptr splinePvt(new PvtDeadSpline); - splinePvt->initFromGas(pvdgTables, samples); - props_[phase_usage_.phase_pos[Vapour]] = splinePvt; - } else { - std::shared_ptr deadPvt(new PvtDead); - deadPvt->initFromGas(pvdgTables); - props_[phase_usage_.phase_pos[Vapour]] = deadPvt; - } - } else if (!pvtgTables.empty()) { - props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(pvtgTables)); - } else { - OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); - } - } - - return props_; -} - -void testmu(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector T, std::vector r,std::vector z, - std::vector > props_, std::vector condition) -{ - std::vector mu(n); - std::vector dmudp(n); - std::vector dmudr(n); - std::vector mu_new(n); - double dmudp_diff; - double dmudr_diff; - double dmudp_diff_u; - double dmudr_diff_u; - - // test mu - for (int phase = 0; phase < np; ++phase) { - props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &T[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]); - props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &T[0], &z[0], &mu[0]); - dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]); - dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]); - dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]); - dmudr_diff_u = (mu_new[5]-mu_new[3])/(r[5]-r[3]); - - for (int i = 0; i < n; ++i){ - BOOST_CHECK_CLOSE(mu_new[i],mu[i],reltol); - } - - // saturated case - BOOST_CHECK_CLOSE(dmudp_diff,dmudp[0],reltol); - BOOST_CHECK_CLOSE(dmudr_diff,dmudr[0],reltol); - - // unsaturated case - BOOST_CHECK_CLOSE(dmudp_diff_u,dmudp[3],reltol); - BOOST_CHECK_CLOSE(dmudr_diff_u,dmudr[3],reltol); - - } -} - -void testb(const double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector T, std::vector r,std::vector z, - std::vector > props_, std::vector condition) -{ - // test b - std::vector b(n); - std::vector B(n); - std::vector invB(n); - std::vector dinvBdp(n); - std::vector dBdp(n); - std::vector dbdr(n); - std::vector dbdp(n); - double dbdp_diff; - double dbdr_diff; - double dbdp_diff_u; - double dbdr_diff_u; - - for (int phase = 0; phase < np; ++phase) { - props_[phase]->b(n, &pvtTableIdx[0], &p[0], &T[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]); - props_[phase]->dBdp(n, &pvtTableIdx[0], &p[0], &T[0], &z[0], &B[0], &dBdp[0]); - dbdp_diff = (b[1]-b[0])/(p[1]-p[0]); - dbdr_diff = (b[2]-b[0])/(r[2]-r[0]); - dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]); - dbdr_diff_u = (b[5]-b[3])/(r[5]-r[3]); - for (int i = 0; i < n; ++i){ - invB[i] = 1/B[i]; - dinvBdp[i] = -1/pow(B[i],2) * dBdp[i]; - } - - for (int i = 0; i < n; ++i){ - BOOST_CHECK_CLOSE(invB[i],b[i] , reltol); - BOOST_CHECK_CLOSE(dinvBdp[i],dbdp[i] , reltol); - - } - // saturated case - BOOST_CHECK_CLOSE(dbdp_diff,dbdp[0], reltol); - BOOST_CHECK_CLOSE(dbdr_diff,dbdr[0], reltol); - - // unsaturated case - BOOST_CHECK_CLOSE(dbdp_diff_u,dbdp[3], reltol); - BOOST_CHECK_CLOSE(dbdr_diff_u,dbdr[3], reltol); - } -} - -void testrsSat(double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector > props_){ - // test bublepoint pressure - std::vector rs(n); - std::vector drsdp(n); - double drsdp_diff; - double drsdp_diff_u; - - for (int phase = 0; phase < np; ++phase) { - props_[phase] ->rsSat(n, &pvtTableIdx[0], &p[0], &rs[0], &drsdp[0]); - - drsdp_diff = (rs[1]-rs[0])/(p[1]-p[0]); - drsdp_diff_u = (rs[4]-rs[3])/(p[4]-p[3]); - - // saturated case - BOOST_CHECK_CLOSE(drsdp_diff,drsdp[0], reltol); - - // unsaturad case - BOOST_CHECK_CLOSE(drsdp_diff_u,drsdp[3], reltol); - - } -} - -void testrvSat(double reltol, int n, int np, const std::vector &pvtTableIdx, std::vector p, std::vector > props_){ - // test rv saturated - std::vector rv(n); - std::vector drvdp(n); - double drvdp_diff; - double drvdp_diff_u; - - for (int phase = 0; phase < np; ++phase) { - props_[phase] ->rvSat(n, &pvtTableIdx[0], &p[0], &rv[0], &drvdp[0]); - - drvdp_diff = (rv[1]-rv[0])/(p[1]-p[0]); - drvdp_diff_u = (rv[4]-rv[3])/(p[4]-p[3]); - - // saturated case - BOOST_CHECK_CLOSE(drvdp_diff,drvdp[0], reltol); - - // unsaturad case - BOOST_CHECK_CLOSE(drvdp_diff_u,drvdp[3], reltol); - - } -} - -BOOST_AUTO_TEST_CASE(test_liveoil) -{ - - - // read eclipse deck - const std::string filename = "liveoil.DATA"; - cout << "Reading deck: " << filename << endl; - Opm::ParserPtr parser(new Opm::Parser()); - Opm::ParseMode parseMode; - Opm::DeckConstPtr deck(parser->parseFile(filename , parseMode)); - Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck, parseMode)); - - // setup pvt interface - PhaseUsage phase_usage_ = phaseUsageFromDeck(deck); - std::vector > props_ = getProps(deck, eclipseState, phase_usage_); - - - // setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference - // approximation of the derivatives. - const int n = 6; - const int np = phase_usage_.num_phases; - std::vector pvtRegionIdx(n, 0); - - // the relative tolerance in percentage for acceptable difference in values - const double reltolper = 1e-9; - // the relative tolerance in percentage for acceptable difference in values for viscosity - const double reltolpermu = 1e-1; - - std::vector p(n); - std::vector T(n, 273.15 + 20); - std::vector r(n); - std::vector condition(n); - std::vector z(n * np); - - - // Used for forward difference calculations - const double h_p = 1e4; - const double h_r = 1; - - // saturated - p[0] = 1e7; - p[1] = p[0] + h_p; - p[2] = p[0]; - - r[0] = 200; - r[1] = 200; - r[2] = 200 + h_r; - - condition[0].setFreeGas(); - condition[1].setFreeGas(); - condition[2].setFreeGas(); - - - // undersaturated - p[3] = p[0]; - p[4] = p[1]; - p[5] = p[2]; - - r[3] = 50; - r[4] = 50; - r[5] = 50 +h_r; - - // Corresponing z factors, used to compare with the [p,z] interface - for (int i = 0; i < n; ++i) { - z[0+i*np] = 0; z[1+i*np] = 1; - z[2+i*np] = r[i]; - - } - - testmu(reltolpermu, n, np, pvtRegionIdx, p, T, r,z, props_, condition); - - testb(reltolper,n,np,pvtRegionIdx,p,T,r,z,props_,condition); - - testrsSat(reltolper,n,np,pvtRegionIdx,p,props_); - - testrvSat(reltolper,n,np,pvtRegionIdx,p,props_); - - -} - -BOOST_AUTO_TEST_CASE(test_wetgas) -{ - - - // read eclipse deck - - const std::string filename = "wetgas.DATA"; - cout << "Reading deck: " << filename << endl; - Opm::ParserPtr parser(new Opm::Parser()); - Opm::ParseMode parseMode; - Opm::DeckConstPtr deck(parser->parseFile(filename , parseMode)); - Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck, parseMode)); - - // setup pvt interface - PhaseUsage phase_usage_ = phaseUsageFromDeck(deck); - std::vector > props_ = getProps(deck,eclipseState,phase_usage_); - - - // setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference - // approximation of the derivatives. - const int n = 6; - const int np = phase_usage_.num_phases; - std::vector pvtRegionIdx(n, 0); - - // the relative tolerance in percentage for acceptable difference in values - const double reltolper = 1e-9; - // the relative tolerance in percentage for acceptable difference in values for viscosity - const double reltolpermu = 1e-1; - - std::vector p(n); - std::vector T(n, 273.15+20); - std::vector r(n); - std::vector condition(n); - std::vector z(n * np); - - - // Used for forward difference calculations - const double h_p = 1e4; - const double h_r = 1e-7; - - // saturated - p[0] = 1e7; - p[1] = p[0] + h_p; - p[2] = p[0]; - - r[0] = 5e-5; - r[1] = 5e-5; - r[2] = 5e-5 + h_r; - - condition[0].setFreeOil(); - condition[1].setFreeOil(); - condition[2].setFreeOil(); - - - // undersaturated - p[3] = p[0]; - p[4] = p[1]; - p[5] = p[2]; - - r[3] = 1e-5; - r[4] = 1e-5; - r[5] = 1e-5 +h_r; - - // Corresponing z factors, used to compare with the [p,z] interface - for (int i = 0; i < n; ++i) { - z[0+i*np] = 0; z[1+i*np] = r[i]; - z[2+i*np] = 1; - - } - - testmu(reltolpermu, n, np, pvtRegionIdx, p,T, r,z, props_, condition); - - testb(reltolper,n,np,pvtRegionIdx,p,T,r,z,props_,condition); - - testrsSat(reltolper,n,np,pvtRegionIdx,p,props_); - - testrvSat(reltolper,n,np,pvtRegionIdx,p,props_); - -} diff --git a/tests/test_norne_pvt.cpp b/tests/test_norne_pvt.cpp index 39a9c868..57e28fcd 100644 --- a/tests/test_norne_pvt.cpp +++ b/tests/test_norne_pvt.cpp @@ -39,7 +39,7 @@ #include #include -#include +#include using namespace Opm; @@ -57,19 +57,9 @@ using namespace Opm; further semantic meaning. */ - -void check_vectors( const std::vector& v1 , const std::vector& v2) { - double tol = 1e-5; - for (decltype(v1.size()) i = 0; i < v1.size(); i++) { - BOOST_CHECK_CLOSE( v1[i] , v2[i] , tol ); - } -} - - - -void verify_norne_oil_pvt_region2(const TableManager& tableManager) { - auto pvtoTables = tableManager.getPvtoTables(); - PvtLiveOil oilPvt(pvtoTables); +void verify_norne_oil_pvt_region1(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState) { + Opm::LiveOilPvt oilPvt; + oilPvt.initFromDeck(deck, eclState); std::vector rs = {33, 33, 43, 43, @@ -113,33 +103,37 @@ void verify_norne_oil_pvt_region2(const TableManager& tableManager) { { std::vector tableIndex(P.size() , 0); - std::vector mu(P.size()); - std::vector dmudp(P.size()); - std::vector dmudr(P.size()); - - std::vector b(P.size()); - std::vector dbdp(P.size()); - std::vector dbdr(P.size()); - - + // convert the pressures to SI units (bar to Pascal) for (auto& value : P) value *= Metric::Pressure; + // convert the gas dissolution factors to SI units for (auto& value : rs) value *= Metric::GasDissolutionFactor; - oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data()); - oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data()); + for (unsigned i = 0; i < P.size(); ++i) { + double mu; + double b; + double RsSat = oilPvt.saturatedGasDissolutionFactor(/*tableIndex=*/0, /*T=*/273.15, P[i]); + if (rs[i] >= RsSat) { + mu = oilPvt.saturatedViscosity(/*tableIndex=*/0, /*T=*/273.15, P[i]); + b = oilPvt.saturatedInverseFormationVolumeFactor(/*tableIndex=*/0, /*T=*/273.15, P[i]); + } + else { + mu = oilPvt.viscosity(/*tableIndex=*/0, /*T=*/273.15, P[i], rs[i]); + b = oilPvt.inverseFormationVolumeFactor(/*tableIndex=*/0, /*T=*/273.15, P[i], rs[i]); + } - check_vectors( mu , mu_expected ); - check_vectors( b , b_expected ); + BOOST_CHECK_CLOSE( mu , mu_expected[i], 1e-5 ); + BOOST_CHECK_CLOSE( b , b_expected[i], 1e-5 ); + } } } -void verify_norne_oil_pvt_region1(const TableManager& tableManager) { - auto pvtoTables = tableManager.getPvtoTables(); - PvtLiveOil oilPvt(pvtoTables); +void verify_norne_oil_pvt_region2(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState) { + Opm::LiveOilPvt oilPvt; + oilPvt.initFromDeck(deck, eclState); std::vector rs = {21 , 21, 30 , 30, @@ -254,50 +248,41 @@ void verify_norne_oil_pvt_region1(const TableManager& tableManager) { 0.57289091037, 0.56019050084, 0.55474601877, 0.55809201119, 0.54526832277}; - { - std::vector tableIndex(P.size() , 1); + // convert the pressures to SI units (bar to Pascal) + for (auto& value : P) + value *= Metric::Pressure; - std::vector mu(P.size()); - std::vector dmudp(P.size()); - std::vector dmudr(P.size()); + // convert the gas dissolution factors to SI units + for (auto& value : rs) + value *= Metric::GasDissolutionFactor; - std::vector b(P.size()); - std::vector dbdp(P.size()); - std::vector dbdr(P.size()); + for (unsigned i = 0; i < P.size(); ++i) { + double mu; + double b; + double RsSat = oilPvt.saturatedGasDissolutionFactor(/*tableIndex=*/1, /*T=*/273.15, P[i]); + if (rs[i] >= RsSat) { + mu = oilPvt.saturatedViscosity(/*tableIndex=*/1, /*T=*/273.15, P[i]); + b = oilPvt.saturatedInverseFormationVolumeFactor(/*tableIndex=*/1, /*T=*/273.15, P[i]); + } + else { + mu = oilPvt.viscosity(/*tableIndex=*/1, /*T=*/273.15, P[i], rs[i]); + b = oilPvt.inverseFormationVolumeFactor(/*tableIndex=*/1, /*T=*/273.15, P[i], rs[i]); + } - - for (auto& value : P) - value *= Metric::Pressure; - - for (auto& value : rs) - value *= Metric::GasDissolutionFactor; - - oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data()); - oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data()); - - check_vectors( mu , mu_expected ); - check_vectors( b , b_expected ); + BOOST_CHECK_CLOSE( mu , mu_expected[i], 1e-5 ); + BOOST_CHECK_CLOSE( b , b_expected[i], 1e-5 ); } } - - - - -TableManager loadTables( const std::string& deck_file) { +BOOST_AUTO_TEST_CASE( Test_Norne_PVT) { Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); Opm::ParserPtr parser(new Parser()); + std::shared_ptr deck; + deck = parser->parseFile("norne_pvt.data", parseMode); - deck = parser->parseFile(deck_file, parseMode); - return TableManager(*deck); -} - - - -BOOST_AUTO_TEST_CASE( Test_Norne_PVT) { - TableManager tableManager = loadTables( "norne_pvt.data" ); - - verify_norne_oil_pvt_region1( tableManager ); - verify_norne_oil_pvt_region2( tableManager ); + Opm::EclipseStateConstPtr eclState(new EclipseState(deck, parseMode)); + + verify_norne_oil_pvt_region1( deck, eclState ); + verify_norne_oil_pvt_region2( deck, eclState ); }