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

Merge pull request #203 from atgeirr/combined

Browse Source 
Reorganisation and improvements in significant parts of opm-core
This commit is contained in:
Bård Skaflestad
2013-03-18 12:12:55 -07:00
parent e59db97128 b32674f3ea
commit 36b3f46e93
148 changed files with 1139 additions and 4040 deletions
Download Patch File Download Diff File Expand all files Collapse all files

353
CMakeLists_files.cmake
View File

@@ -26,56 +26,34 @@
# originally generated with the command:
# find opm -name '*.c*' -printf '\t%p\n' | sort
list (APPEND MAIN_SOURCE_FILES
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp
opm/core/fluid/BlackoilPropertiesBasic.cpp
opm/core/fluid/BlackoilPropertiesFromDeck.cpp
opm/core/fluid/blackoil/SinglePvtDead.cpp
opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp
opm/core/fluid/blackoil/SinglePvtInterface.cpp
opm/core/fluid/blackoil/SinglePvtLiveGas.cpp
opm/core/fluid/blackoil/SinglePvtLiveOil.cpp
opm/core/fluid/IncompPropertiesBasic.cpp
opm/core/fluid/IncompPropertiesFromDeck.cpp
opm/core/fluid/PvtPropertiesBasic.cpp
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp
opm/core/fluid/RockBasic.cpp
opm/core/fluid/RockCompressibility.cpp
opm/core/fluid/RockFromDeck.cpp
opm/core/fluid/SatFuncGwseg.cpp
opm/core/fluid/SatFuncSimple.cpp
opm/core/fluid/SatFuncStone2.cpp
opm/core/fluid/SaturationPropsBasic.cpp
opm/core/fluid/SaturationPropsFromDeck.cpp
opm/core/grid.c
opm/core/grid/GridManager.cpp
opm/core/grid/grid.c
opm/core/grid/cart_grid.c
opm/core/grid/cornerpoint_grid.c
opm/core/grid/cpgpreprocess/facetopology.c
opm/core/grid/cpgpreprocess/geometry.c
opm/core/grid/cpgpreprocess/preprocess.c
opm/core/grid/cpgpreprocess/processgrid.c
opm/core/grid/cpgpreprocess/uniquepoints.c
opm/core/GridManager.cpp
opm/core/io/eclipse/EclipseGridInspector.cpp
opm/core/io/eclipse/EclipseGridParser.cpp
opm/core/io/eclipse/writeECLData.cpp
opm/core/io/vag/vag.cpp
opm/core/io/vtk/writeVtkData.cpp
opm/core/linalg/call_umfpack.c
opm/core/linalg/LinearSolverAGMG.cpp
opm/core/linalg/LinearSolverFactory.cpp
opm/core/linalg/LinearSolverInterface.cpp
opm/core/linalg/LinearSolverIstl.cpp
opm/core/linalg/LinearSolverUmfpack.cpp
opm/core/linalg/call_umfpack.c
opm/core/linalg/sparse_sys.c
opm/core/newwells.c
opm/core/pressure/cfsh.c
opm/core/pressure/CompressibleTpfa.cpp
opm/core/pressure/flow_bc.c
opm/core/pressure/FlowBCManager.cpp
opm/core/pressure/IncompTpfa.cpp
opm/core/pressure/cfsh.c
opm/core/pressure/flow_bc.c
opm/core/pressure/fsh.c
opm/core/pressure/fsh_common_impl.c
opm/core/pressure/ifsh.c
opm/core/pressure/IncompTpfa.cpp
opm/core/pressure/mimetic/hybsys.c
opm/core/pressure/mimetic/hybsys_global.c
opm/core/pressure/mimetic/mimetic.c
@@ -93,26 +71,50 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/pressure/tpfa/compr_source.c
opm/core/pressure/tpfa/ifs_tpfa.c
opm/core/pressure/tpfa/trans_tpfa.c
opm/core/pressure/well.c
opm/core/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/pvt/BlackoilPvtProperties.cpp
opm/core/props/pvt/PvtPropertiesBasic.cpp
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
opm/core/props/pvt/SinglePvtDead.cpp
opm/core/props/pvt/SinglePvtDeadSpline.cpp
opm/core/props/pvt/SinglePvtInterface.cpp
opm/core/props/pvt/SinglePvtLiveGas.cpp
opm/core/props/pvt/SinglePvtLiveOil.cpp
opm/core/props/rock/RockBasic.cpp
opm/core/props/rock/RockCompressibility.cpp
opm/core/props/rock/RockFromDeck.cpp
opm/core/props/satfunc/SatFuncGwseg.cpp
opm/core/props/satfunc/SatFuncSimple.cpp
opm/core/props/satfunc/SatFuncStone2.cpp
opm/core/props/satfunc/SaturationPropsBasic.cpp
opm/core/props/satfunc/SaturationPropsFromDeck.cpp
opm/core/simulator/SimulatorCompressibleTwophase.cpp
opm/core/simulator/SimulatorIncompTwophase.cpp
opm/core/simulator/SimulatorReport.cpp
opm/core/simulator/SimulatorTimer.cpp
opm/core/transport/reorder/DGBasis.cpp
opm/core/transport/reorder/nlsolvers.c
opm/core/tof/DGBasis.cpp
opm/core/tof/TofReorder.cpp
opm/core/tof/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/transport/reorder/TransportModelCompressibleTwophase.cpp
opm/core/transport/reorder/TransportModelInterface.cpp
opm/core/transport/reorder/TransportModelTracerTof.cpp
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp
opm/core/transport/reorder/TransportModelTwophase.cpp
opm/core/transport/spu_explicit.c
opm/core/transport/spu_implicit.c
opm/core/transport/transport_source.c
opm/core/utility/miscUtilitiesBlackoil.cpp
opm/core/utility/miscUtilities.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/parameters/Parameter.cpp
opm/core/utility/parameters/ParameterGroup.cpp
opm/core/utility/parameters/ParameterTools.cpp
@@ -122,14 +124,12 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/utility/parameters/tinyxml/tinyxmlerror.cpp
opm/core/utility/parameters/tinyxml/tinyxmlparser.cpp
opm/core/utility/parameters/tinyxml/xmltest.cpp
opm/core/utility/StopWatch.cpp
opm/core/utility/VelocityInterpolation.cpp
opm/core/utility/WachspressCoord.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
)
# originally generated with the command:
@@ -160,12 +160,8 @@ list (APPEND EXAMPLE_SOURCE_FILES
examples/compute_tof.cpp
examples/compute_tof_from_files.cpp
examples/import_rewrite.cpp
examples/refine_wells.cpp
examples/scaneclipsedeck.c
examples/sim_2p_comp_reorder.cpp
examples/sim_2p_incomp_reorder.cpp
examples/sim_wateroil.cpp
examples/spu_2p.cpp
examples/sim_2p_incomp.cpp
examples/wells_example.cpp
tutorials/tutorial1.cpp
tutorials/tutorial2.cpp
@@ -176,157 +172,154 @@ list (APPEND EXAMPLE_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/transport/ImplicitTransport.hpp
opm/core/transport/SimpleFluid2pWrapper.hpp
opm/core/transport/CSRMatrixBlockAssembler.hpp
opm/core/transport/ImplicitAssembly.hpp
opm/core/transport/CSRMatrixUmfpackSolver.hpp
opm/core/transport/GravityColumnSolver_impl.hpp
opm/core/transport/spu_implicit.h
opm/core/transport/transport_source.h
opm/core/transport/SinglePointUpwindTwoPhase.hpp
opm/core/transport/JacobianSystem.hpp
opm/core/transport/reorder/TransportModelTwophase.hpp
opm/core/transport/reorder/TransportModelTracerTof.hpp
opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp
opm/core/transport/reorder/DGBasis.hpp
opm/core/transport/reorder/reordersequence.h
opm/core/transport/reorder/tarjan.h
opm/core/transport/reorder/TransportModelInterface.hpp
opm/core/transport/reorder/nlsolvers.h
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp
opm/core/transport/GravityColumnSolver.hpp
opm/core/transport/NormSupport.hpp
opm/core/transport/spu_explicit.h
opm/core/GridManager.hpp
opm/core/io/eclipse/SpecialEclipseFields.hpp
opm/core/io/eclipse/EclipseUnits.hpp
opm/core/io/eclipse/EclipseGridParserHelpers.hpp
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/GridManager.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/grdecl.h
opm/core/grid/cpgpreprocess/preprocess.h
opm/core/grid/cpgpreprocess/uniquepoints.h
opm/core/io/eclipse/CornerpointChopper.hpp
opm/core/io/eclipse/EclipseGridParser.hpp
opm/core/io/eclipse/EclipseGridInspector.hpp
opm/core/io/eclipse/EclipseGridParser.hpp
opm/core/io/eclipse/EclipseGridParserHelpers.hpp
opm/core/io/eclipse/EclipseUnits.hpp
opm/core/io/eclipse/SpecialEclipseFields.hpp
opm/core/io/eclipse/writeECLData.hpp
opm/core/io/vag/vag.hpp
opm/core/io/vtk/writeVtkData.hpp
opm/core/linalg/blas_lapack.h
opm/core/linalg/LinearSolverAGMG.hpp
opm/core/linalg/LinearSolverFactory.hpp
opm/core/linalg/LinearSolverUmfpack.hpp
opm/core/linalg/LinearSolverInterface.hpp
opm/core/linalg/LinearSolverIstl.hpp
opm/core/linalg/LinearSolverUmfpack.hpp
opm/core/linalg/blas_lapack.h
opm/core/linalg/call_umfpack.h
opm/core/linalg/sparse_sys.h
opm/core/linalg/LinearSolverInterface.hpp
opm/core/newwells.h
opm/core/doxygen_main.hpp
opm/core/simulator/SimulatorIncompTwophase.hpp
opm/core/simulator/WellState.hpp
opm/core/wells.h
opm/core/pressure/CompressibleTpfa.hpp
opm/core/pressure/FlowBCManager.hpp
opm/core/pressure/IncompTpfa.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/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/phaseUsageFromDeck.hpp
opm/core/props/pvt/BlackoilPvtProperties.hpp
opm/core/props/pvt/PvtPropertiesBasic.hpp
opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
opm/core/props/pvt/SinglePvtConstCompr.hpp
opm/core/props/pvt/SinglePvtDead.hpp
opm/core/props/pvt/SinglePvtDeadSpline.hpp
opm/core/props/pvt/SinglePvtInterface.hpp
opm/core/props/pvt/SinglePvtLiveGas.hpp
opm/core/props/pvt/SinglePvtLiveOil.hpp
opm/core/props/rock/RockBasic.hpp
opm/core/props/rock/RockCompressibility.hpp
opm/core/props/rock/RockFromDeck.hpp
opm/core/props/satfunc/SatFuncGwseg.hpp
opm/core/props/satfunc/SatFuncSimple.hpp
opm/core/props/satfunc/SatFuncStone2.hpp
opm/core/props/satfunc/SaturationPropsBasic.hpp
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
opm/core/props/satfunc/SaturationPropsInterface.hpp
opm/core/simulator/BlackoilState.hpp
opm/core/simulator/SimulatorCompressibleTwophase.hpp
opm/core/simulator/SimulatorIncompTwophase.hpp
opm/core/simulator/SimulatorReport.hpp
opm/core/simulator/SimulatorTimer.hpp
opm/core/simulator/TwophaseState.hpp
opm/core/simulator/BlackoilState.hpp
opm/core/utility/linearInterpolation.hpp
opm/core/utility/WachspressCoord.hpp
opm/core/utility/initState_impl.hpp
opm/core/utility/initState.hpp
opm/core/utility/NonuniformTableLinear.hpp
opm/core/utility/Units.hpp
opm/core/utility/DataMap.hpp
opm/core/utility/linInt.hpp
opm/core/utility/parameters/ParameterXML.hpp
opm/core/utility/parameters/ParameterGroup_impl.hpp
opm/core/utility/parameters/Parameter.hpp
opm/core/utility/parameters/ParameterTools.hpp
opm/core/utility/parameters/ParameterRequirement.hpp
opm/core/utility/parameters/ParameterMapItem.hpp
opm/core/utility/parameters/ParameterGroup.hpp
opm/core/utility/parameters/tinyxml/tinystr.h
opm/core/utility/parameters/tinyxml/tinyxml.h
opm/core/utility/parameters/ParameterStrings.hpp
opm/core/utility/miscUtilitiesBlackoil.hpp
opm/core/utility/MonotCubicInterpolator.hpp
opm/core/utility/UniformTableLinear.hpp
opm/core/simulator/WellState.hpp
opm/core/simulator/initState.hpp
opm/core/simulator/initState_impl.hpp
opm/core/tof/DGBasis.hpp
opm/core/tof/TofReorder.hpp
opm/core/tof/TofDiscGalReorder.hpp
opm/core/transport/TransportSolverTwophaseInterface.hpp
opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp
opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp
opm/core/transport/implicit/GravityColumnSolver.hpp
opm/core/transport/implicit/GravityColumnSolver_impl.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/Factory.hpp
opm/core/utility/VelocityInterpolation.hpp
opm/core/utility/SparseVector.hpp
opm/core/utility/StopWatch.hpp
opm/core/utility/buildUniformMonotoneTable.hpp
opm/core/utility/DataMap.hpp
opm/core/utility/ErrorMacros.hpp
opm/core/utility/miscUtilities.hpp
opm/core/utility/ColumnExtract.hpp
opm/core/utility/Factory.hpp
opm/core/utility/MonotCubicInterpolator.hpp
opm/core/utility/NonuniformTableLinear.hpp
opm/core/utility/RootFinders.hpp
opm/core/utility/SparseTable.hpp
opm/core/utility/SparseVector.hpp
opm/core/utility/StopWatch.hpp
opm/core/utility/UniformTableLinear.hpp
opm/core/utility/Units.hpp
opm/core/utility/VelocityInterpolation.hpp
opm/core/utility/WachspressCoord.hpp
opm/core/utility/buildUniformMonotoneTable.hpp
opm/core/utility/have_boost_redef.hpp
opm/core/grid/CellQuadrature.hpp
opm/core/grid/cornerpoint_grid.h
opm/core/grid/cpgpreprocess/geometry.h
opm/core/grid/cpgpreprocess/uniquepoints.h
opm/core/grid/cpgpreprocess/preprocess.h
opm/core/grid/cpgpreprocess/grdecl.h
opm/core/grid/cpgpreprocess/facetopology.h
opm/core/grid/cart_grid.h
opm/core/grid/FaceQuadrature.hpp
opm/core/pressure/TPFAPressureSolver.hpp
opm/core/pressure/fsh_common_impl.h
opm/core/pressure/msmfem/dfs.h
opm/core/pressure/msmfem/coarse_conn.h
opm/core/pressure/msmfem/hash_set.h
opm/core/pressure/msmfem/partition.h
opm/core/pressure/msmfem/ifsh_ms.h
opm/core/pressure/msmfem/coarse_sys.h
opm/core/pressure/IncompTpfa.hpp
opm/core/pressure/HybridPressureSolver.hpp
opm/core/pressure/fsh.h
opm/core/pressure/flow_bc.h
opm/core/pressure/tpfa/compr_quant_general.h
opm/core/pressure/tpfa/compr_bc.h
opm/core/pressure/tpfa/trans_tpfa.h
opm/core/pressure/tpfa/compr_quant.h
opm/core/pressure/tpfa/cfs_tpfa_residual.h
opm/core/pressure/tpfa/ifs_tpfa.h
opm/core/pressure/tpfa/compr_source.h
opm/core/pressure/tpfa/cfs_tpfa.h
opm/core/pressure/mimetic/hybsys_global.h
opm/core/pressure/mimetic/hybsys.h
opm/core/pressure/mimetic/mimetic.h
opm/core/pressure/FlowBCManager.hpp
opm/core/pressure/CompressibleTpfa.hpp
opm/core/pressure/TPFACompressiblePressureSolver.hpp
opm/core/utility/linInt.hpp
opm/core/utility/linearInterpolation.hpp
opm/core/utility/miscUtilities.hpp
opm/core/utility/miscUtilitiesBlackoil.hpp
opm/core/utility/parameters/Parameter.hpp
opm/core/utility/parameters/ParameterGroup.hpp
opm/core/utility/parameters/ParameterGroup_impl.hpp
opm/core/utility/parameters/ParameterMapItem.hpp
opm/core/utility/parameters/ParameterRequirement.hpp
opm/core/utility/parameters/ParameterStrings.hpp
opm/core/utility/parameters/ParameterTools.hpp
opm/core/utility/parameters/ParameterXML.hpp
opm/core/utility/parameters/tinyxml/tinystr.h
opm/core/utility/parameters/tinyxml/tinyxml.h
opm/core/wells/InjectionSpecification.hpp
opm/core/wells/WellsGroup.hpp
opm/core/wells/WellsManager.hpp
opm/core/wells/ProductionSpecification.hpp
opm/core/wells/WellCollection.hpp
opm/core/fluid/SaturationPropsInterface.hpp
opm/core/fluid/SaturationPropsFromDeck_impl.hpp
opm/core/fluid/RockFromDeck.hpp
opm/core/fluid/BlackoilPropertiesBasic.hpp
opm/core/fluid/RockCompressibility.hpp
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp
opm/core/fluid/IncompPropertiesBasic.hpp
opm/core/fluid/SaturationPropsFromDeck.hpp
opm/core/fluid/SaturationPropsBasic.hpp
opm/core/fluid/BlackoilPropertiesFromDeck.hpp
opm/core/fluid/BlackoilPropertiesInterface.hpp
opm/core/fluid/blackoil/SinglePvtInterface.hpp
opm/core/fluid/blackoil/BlackoilPhases.hpp
opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp
opm/core/fluid/blackoil/SinglePvtLiveOil.hpp
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp
opm/core/fluid/blackoil/SinglePvtDead.hpp
opm/core/fluid/blackoil/phaseUsageFromDeck.hpp
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp
opm/core/fluid/blackoil/SinglePvtLiveGas.hpp
opm/core/fluid/SimpleFluid2p.hpp
opm/core/fluid/IncompPropertiesInterface.hpp
opm/core/fluid/PvtPropertiesBasic.hpp
opm/core/fluid/SatFuncSimple.hpp
opm/core/fluid/SatFuncGwseg.hpp
opm/core/fluid/IncompPropertiesFromDeck.hpp
opm/core/fluid/RockBasic.hpp
opm/core/fluid/SatFuncStone2.hpp
opm/core/grid.h
opm/core/well.h
opm/core/GridAdapter.hpp
opm/core/wells/WellsGroup.hpp
opm/core/wells/WellsManager.hpp
)

20
examples/compute_tof.cpp
View File

@@ -26,25 +26,25 @@
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/tof/TofReorder.hpp>
#include <opm/core/tof/TofDiscGalReorder.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
@@ -170,9 +170,9 @@ main(int argc, char** argv)
bool use_dg = param.getDefault("use_dg", false);
bool use_multidim_upwind = false;
// Need to initialize dg solver here, since it uses parameters now.
boost::scoped_ptr<Opm::TransportModelTracerTofDiscGal> dg_solver;
boost::scoped_ptr<Opm::TofDiscGalReorder> dg_solver;
if (use_dg) {
dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(*grid->c_grid(), param));
dg_solver.reset(new Opm::TofDiscGalReorder(*grid->c_grid(), param));
} else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
}
@@ -237,7 +237,7 @@ main(int argc, char** argv)
if (use_dg) {
dg_solver->solveTof(&state.faceflux()[0], &porevol[0], &transport_src[0], tof);
} else {
Opm::TransportModelTracerTof tofsolver(*grid->c_grid(), use_multidim_upwind);
Opm::TofReorder tofsolver(*grid->c_grid(), use_multidim_upwind);
if (compute_tracer) {
tofsolver.solveTofTracer(&state.faceflux()[0], &porevol[0], &transport_src[0], tof, tracer);
} else {

20
examples/compute_tof_from_files.cpp
View File

@@ -25,25 +25,25 @@
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/tof/TofReorder.hpp>
#include <opm/core/tof/TofDiscGalReorder.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/filesystem.hpp>
@@ -124,9 +124,9 @@ main(int argc, char** argv)
bool use_dg = param.getDefault("use_dg", false);
bool use_multidim_upwind = false;
// Need to initialize dg solver here, since it uses parameters now.
boost::scoped_ptr<Opm::TransportModelTracerTofDiscGal> dg_solver;
boost::scoped_ptr<Opm::TofDiscGalReorder> dg_solver;
if (use_dg) {
dg_solver.reset(new Opm::TransportModelTracerTofDiscGal(grid, param));
dg_solver.reset(new Opm::TofDiscGalReorder(grid, param));
} else {
use_multidim_upwind = param.getDefault("use_multidim_upwind", false);
}
@@ -163,7 +163,7 @@ main(int argc, char** argv)
if (use_dg) {
dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof);
} else {
Opm::TransportModelTracerTof tofsolver(grid, use_multidim_upwind);
Opm::TofReorder tofsolver(grid, use_multidim_upwind);
if (compute_tracer) {
tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tof, tracer);
} else {

50
examples/refine_wells.cpp
View File

@@ -1,50 +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 <http://www.gnu.org/licenses/>.
*/
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
// Double I and J coordinates of wells and completions.
// Do not change any well productivity indices or any other data.
int main(int argc, char** argv)
{
using namespace Opm;
if (argc != 2) {
std::cerr << "Usage: " << argv[0] << " deck\n";
return 1;
}
EclipseGridParser deck(argv[1], false);
WELSPECS ws = deck.getWELSPECS();
const int nw = ws.welspecs.size();
for (int w = 0; w < nw; ++w) {
ws.welspecs[w].I_ *= 2;
ws.welspecs[w].J_ *= 2;
}
COMPDAT cd = deck.getCOMPDAT();
const int nc = cd.compdat.size();
for (int c = 0; c < nc; ++c) {
cd.compdat[c].grid_ind_[0] *= 2;
cd.compdat[c].grid_ind_[1] *= 2;
}
ws.write(std::cout);
std::cout << '\n';
cd.write(std::cout);
}

134
examples/scaneclipsedeck.c
View File

@@ -1,134 +0,0 @@
/* scaneclipsedec finds technically valid Eclipse keywords in an ascii file.
* Copyright (c) 2010 Jostein R. Natvig <jostein.natvig@gmail.com>
*
* The MIT License
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
* OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
* OTHER DEALINGS IN THE SOFTWARE.
*/
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
static char*
read_keyword(FILE *fp, char *buf)
{
int i, j, c;
/* Clear buf */
for (i=0; i<9; ++i) {
buf[i] = '\0';
}
/* Read first character and check if it is uppercase*/
buf[0] = fgetc(fp);
if ( !isupper( buf[0] ) ) {
return NULL; /* NOT VALID CHARACTER */
ungetc(buf[0], fp);
}
/* Scan as much as possible possible keyword, 8 characters long */
i = 1;
while ( (c = fgetc(fp)) &&
(c != EOF ) &&
(c != '\n' ) &&
(c != '/' ) &&
(i < 8 )) {
buf[i++] = c;
}
/* Skip rest of line */
if (c != '\n'){
while ( (c = fgetc(fp)) &&
(c != EOF ) &&
(c != '\n' )) {
;
}
}
if(c == '\n') {
ungetc(c, fp);
}
/* Find first non-uppercase or non-digit character */
for (i=0; i<8; ++i) {
if ( !(isupper(buf[i]) || isdigit(buf[i])) ) {
break;
}
}
/* Check if remaining characters are blank */
for (j = i; j<8; ++j) {
if(!isspace(buf[j]) && buf[j] != '\0') {
return NULL; /* CHARACTER AFTER SPACE OR INVALID CHARACTER */
}
buf[j] = '\0';
}
return buf;
}
int main(int argc, char *argv[])
{
int c, lineno, nkw;
FILE *fp;
char buf[10];
if (argc != 2)
{
fprintf(stderr, "Usage: <app> filename.grdecl\n");
exit(EXIT_FAILURE);
}
fp = fopen(argv[1], "ra");
if (fp == NULL)
{
fprintf(stderr, "No such file...\n");
exit(EXIT_FAILURE);
}
lineno = nkw = 0;
if (read_keyword(fp, buf) != NULL) {
++nkw;
fprintf(stderr, "%s\n", buf);
}
while ((c = getc(fp)) != EOF) { /* Eat large chunks */
if ( c == '\n') {
++lineno;
if (read_keyword(fp, buf) != NULL) {
fprintf(stderr, "%s\n", buf);
++nkw;
}
}
}
fprintf(stderr, "Scanned %d lines, found %d keywords.\n", lineno, nkw);
return 0;
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

12
examples/sim_2p_comp_reorder.cpp
View File

@@ -24,19 +24,19 @@
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/BlackoilPropertiesBasic.hpp>
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>

12
examples/sim_2p_incomp_reorder.cpp → examples/sim_2p_incomp.cpp
View File

@@ -25,19 +25,19 @@
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorReport.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>

477
examples/sim_wateroil.cpp
View File

@@ -1,477 +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 <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/pressure/CompressibleTpfa.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/BlackoilPropertiesBasic.hpp>
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/GravityColumnSolver.hpp>
#include <opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <cassert>
#include <cstddef>
#include <algorithm>
#include <tr1/array>
#include <functional>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iterator>
#include <vector>
#include <numeric>
template <class State>
static void outputState(const UnstructuredGrid& grid,
const State& state,
const int step,
const std::string& output_dir)
{
// Write data in VTK format.
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
{
// Write water cut curve.
std::string fname = output_dir + "/watercut.txt";
std::ofstream os(fname.c_str());
if (!os) {
THROW("Failed to open " << fname);
}
watercut.write(os);
}
static void outputWellReport(const Opm::WellReport& wellreport,
const std::string& output_dir)
{
// Write well report.
std::string fname = output_dir + "/wellreport.txt";
std::ofstream os(fname.c_str());
if (!os) {
THROW("Failed to open " << fname);
}
wellreport.write(os);
}
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
using namespace Opm;
std::cout << "\n================ Test program for weakly compressible two-phase flow ===============\n\n";
Opm::parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
// Reading various control parameters.
const bool output = param.getDefault("output", true);
std::string output_dir;
int output_interval = 1;
if (output) {
output_dir = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
output_interval = param.getDefault("output_interval", output_interval);
}
const int num_transport_substeps = param.getDefault("num_transport_substeps", 1);
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
boost::scoped_ptr<Opm::GridManager> grid;
boost::scoped_ptr<Opm::BlackoilPropertiesInterface> props;
boost::scoped_ptr<Opm::WellsManager> wells;
boost::scoped_ptr<Opm::RockCompressibility> rock_comp;
Opm::SimulatorTimer simtimer;
Opm::BlackoilState state;
bool check_well_controls = false;
int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename);
// Grid init
grid.reset(new Opm::GridManager(deck));
// Rock and fluid init
props.reset(new BlackoilPropertiesFromDeck(deck, *grid->c_grid(), param));
// Wells init.
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
check_well_controls = param.getDefault("check_well_controls", false);
max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Timer init.
if (deck.hasField("TSTEP")) {
simtimer.init(deck);
} else {
simtimer.init(param);
}
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(deck));
// Gravity.
gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
initBlackoilSurfvol(*grid->c_grid(), *props, state);
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0);
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new Opm::BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Wells init.
wells.reset(new Opm::WellsManager());
// Timer init.
simtimer.init(param);
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
}
// Warn if gravity but no density difference.
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
if (use_gravity) {
if (props->surfaceDensity()[0] == props->surfaceDensity()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
}
bool use_segregation_split = false;
if (use_gravity) {
use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
}
// Source-related variables init.
int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> totmob;
std::vector<double> omega; // Will remain empty if no gravity.
std::vector<double> rc; // Will remain empty if no rock compressibility.
// Extra rock init.
std::vector<double> porevol;
if (rock_comp->isActive()) {
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
} else {
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
}
std::vector<double> initial_porevol = porevol;
double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
// Initialising src
std::vector<double> src(num_cells, 0.0);
if (wells->c_wells()) {
// Do nothing, wells will be the driving force, not source terms.
// Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
} else {
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
*tot_porevol_init/Opm::unit::day;
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
}
std::vector<double> reorder_src = src;
// Solvers init.
// Linear solver.
Opm::LinearSolverFactory linsolver(param);
// Pressure solver.
const double nl_press_res_tol = param.getDefault("nl_press_res_tol", 1e-6);
const double nl_press_change_tol = param.getDefault("nl_press_change_tol", 10.0);
const int nl_press_maxiter = param.getDefault("nl_press_maxiter", 20);
const double *grav = use_gravity ? &gravity[0] : 0;
Opm::CompressibleTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
nl_press_res_tol, nl_press_change_tol, nl_press_maxiter,
grav, wells->c_wells());
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportModelCompressibleTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
if (use_segregation_split) {
reorder_model.initGravity(grav);
}
// Column-based gravity segregation solver.
std::vector<std::vector<int> > columns;
if (use_segregation_split) {
Opm::extractColumn(*grid->c_grid(), columns);
}
// The allcells vector is used in calls to computeTotalMobility()
// and computeTotalMobilityOmega().
std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell;
}
// Warn if any parameters are unused.
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
// Write parameters used for later reference.
if (output) {
param.writeParam(output_dir + "/spu_2p.param");
}
// Main simulation loop.
Opm::time::StopWatch pressure_timer;
double ptime = 0.0;
Opm::time::StopWatch transport_timer;
double ttime = 0.0;
Opm::time::StopWatch total_timer;
total_timer.start();
std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl;
double init_satvol[2] = { 0.0 };
double satvol[2] = { 0.0 };
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init
<< " " << init_satvol[1]/tot_porevol_init << std::endl;
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport;
Opm::WellState well_state;
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
int num_wells = 0;
if (wells->c_wells()) {
num_wells = wells->c_wells()->number_of_wells;
well_state.init(wells->c_wells(), state);
well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(num_wells), 0.0);
wellreport.push(*props, *wells->c_wells(),
state.pressure(), state.surfacevol(), state.saturation(),
0.0, well_state.bhp(), well_state.perfRates());
}
for (; !simtimer.done(); ++simtimer) {
// Report timestep and (optionally) write state to disk.
simtimer.report(std::cout);
if (output && (simtimer.currentStepNum() % output_interval == 0)) {
outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir);
}
// Solve pressure.
if (check_well_controls) {
computeFractionalFlow(*props, allcells, state.pressure(), state.surfacevol(), state.saturation(), fractional_flows);
}
if (check_well_controls) {
wells->applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
bool well_control_passed = !check_well_controls;
int well_control_iteration = 0;
do { // Well control outer loop.
pressure_timer.start();
psolver.solve(simtimer.currentStepLength(), state, well_state);
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
if (check_well_controls) {
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
fractional_flows,
well_state.perfRates(),
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells->conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_state.perfRates(), reorder_src);
// Compute new porevolumes after pressure solve, if necessary.
if (rock_comp->isActive()) {
initial_porevol = porevol;
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
}
// Solve transport.
transport_timer.start();
double stepsize = simtimer.currentStepLength();
if (num_transport_substeps != 1) {
stepsize /= double(num_transport_substeps);
std::cout << "Making " << num_transport_substeps << " transport substeps." << std::endl;
}
for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
// Note that for now we do not handle rock compressibility,
// although the transport solver should be able to.
reorder_model.solve(&state.faceflux()[0], &state.pressure()[0],
&porevol[0], &initial_porevol[0], &reorder_src[0], stepsize,
state.saturation(), state.surfacevol());
// Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
if (use_segregation_split) {
reorder_model.solveGravity(columns,
stepsize, state.saturation(), state.surfacevol());
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
std::cout.precision(5);
const int width = 18;
std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n";
std::cout << " Saturated volumes: "
<< std::setw(width) << satvol[0]/tot_porevol_init
<< std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
std::cout << " Injected volumes: "
<< std::setw(width) << injected[0]/tot_porevol_init
<< std::setw(width) << injected[1]/tot_porevol_init << std::endl;
std::cout << " Produced volumes: "
<< std::setw(width) << produced[0]/tot_porevol_init
<< std::setw(width) << produced[1]/tot_porevol_init << std::endl;
std::cout << " Total inj volumes: "
<< std::setw(width) << tot_injected[0]/tot_porevol_init
<< std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
std::cout << " Total prod volumes: "
<< std::setw(width) << tot_produced[0]/tot_porevol_init
<< std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
std::cout << " In-place + prod - inj: "
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
std::cout << " Init - now - pr + inj: "
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
<< std::endl;
std::cout.precision(8);
watercut.push(simtimer.currentTime() + simtimer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
if (wells->c_wells()) {
wellreport.push(*props, *wells->c_wells(),
state.pressure(), state.surfacevol(), state.saturation(),
simtimer.currentTime() + simtimer.currentStepLength(),
well_state.bhp(), well_state.perfRates());
}
}
total_timer.stop();
std::cout << "\n\n================ End of simulation ===============\n"
<< "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime
<< "\n Transport time: " << ttime << std::endl;
if (output) {
outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir);
outputWaterCut(watercut, output_dir);
if (wells->c_wells()) {
outputWellReport(wellreport, output_dir);
}
}
}

718
examples/spu_2p.cpp
View File

@@ -1,718 +0,0 @@
/*===========================================================================
//
// File: spu_2p.cpp
//
// Created: 2011-10-05 10:29:01+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011, 2012 SINTEF ICT, Applied Mathematics.
Copyright 2011, 2012 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/StopWatch.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/transport/transport_source.h>
#include <opm/core/transport/CSRMatrixUmfpackSolver.hpp>
#include <opm/core/transport/NormSupport.hpp>
#include <opm/core/transport/ImplicitAssembly.hpp>
#include <opm/core/transport/ImplicitTransport.hpp>
#include <opm/core/transport/JacobianSystem.hpp>
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/GravityColumnSolver.hpp>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
#include <cassert>
#include <cstddef>
#include <algorithm>
#include <tr1/array>
#include <functional>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <iterator>
#include <vector>
#include <numeric>
#ifdef HAVE_ERT
#include <opm/core/io/eclipse/writeECLData.hpp>
#endif
static void outputState(const UnstructuredGrid& grid,
const Opm::TwophaseState& state,
const Opm::SimulatorTimer& simtimer,
const std::string& output_dir)
{
// Write data in VTK format.
int step = simtimer.currentStepNum();
std::ostringstream vtkfilename;
vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str());
if (!vtkfile) {
THROW("Failed to open " << vtkfilename.str());
}
Opm::DataMap dm;
dm["saturation"] = &state.saturation();
dm["pressure"] = &state.pressure();
std::vector<double> cell_velocity;
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile);
#ifdef HAVE_ERT
Opm::writeECLData(grid, dm, simtimer.currentStepNum(), simtimer.currentTime(), simtimer.currentDateTime(), output_dir, "OPM" );
#endif
// Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
std::ostringstream fname;
fname << output_dir << "/" << it->first << "-" << std::setw(3) << std::setfill('0') << step << ".dat";
std::ofstream file(fname.str().c_str());
if (!file) {
THROW("Failed to open " << fname.str());
}
const std::vector<double>& d = *(it->second);
std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
}
}
static void outputWaterCut(const Opm::Watercut& watercut,
const std::string& output_dir)
{
// Write water cut curve.
std::string fname = output_dir + "/watercut.txt";
std::ofstream os(fname.c_str());
if (!os) {
THROW("Failed to open " << fname);
}
watercut.write(os);
}
static void outputWellReport(const Opm::WellReport& wellreport,
const std::string& output_dir)
{
// Write well report.
std::string fname = output_dir + "/wellreport.txt";
std::ofstream os(fname.c_str());
if (!os) {
THROW("Failed to open " << fname);
}
wellreport.write(os);
}
// --------------- Types needed to define transport solver ---------------
class SimpleFluid2pWrappingProps
{
public:
SimpleFluid2pWrappingProps(const Opm::IncompPropertiesInterface& props)
: props_(props),
smin_(props.numCells()*props.numPhases()),
smax_(props.numCells()*props.numPhases())
{
if (props.numPhases() != 2) {
THROW("SimpleFluid2pWrapper requires 2 phases.");
}
const int num_cells = props.numCells();
std::vector<int> cells(num_cells);
for (int c = 0; c < num_cells; ++c) {
cells[c] = c;
}
props.satRange(num_cells, &cells[0], &smin_[0], &smax_[0]);
}
double density(int phase) const
{
return props_.density()[phase];
}
template <class Sat,
class Mob,
class DMob>
void mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const
{
props_.relperm(1, &s[0], &c, &mob[0], &dmob[0]);
const double* mu = props_.viscosity();
mob[0] /= mu[0];
mob[1] /= mu[1];
// Recall that we use Fortran ordering for kr derivatives,
// therefore dmob[i*2 + j] is row j and column i of the
// matrix.
// Each row corresponds to a kr function, so which mu to
// divide by also depends on the row, j.
dmob[0*2 + 0] /= mu[0];
dmob[0*2 + 1] /= mu[1];
dmob[1*2 + 0] /= mu[0];
dmob[1*2 + 1] /= mu[1];
}
template <class Sat,
class Pcap,
class DPcap>
void pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const
{
double pcow[2];
double dpcow[4];
props_.capPress(1, &s[0], &c, pcow, dpcow);
pcap = pcow[0];
ASSERT(pcow[1] == 0.0);
dpcap = dpcow[0];
ASSERT(dpcow[1] == 0.0);
ASSERT(dpcow[2] == 0.0);
ASSERT(dpcow[3] == 0.0);
}
double s_min(int c) const
{
return smin_[2*c + 0];
}
double s_max(int c) const
{
return smax_[2*c + 0];
}
private:
const Opm::IncompPropertiesInterface& props_;
std::vector<double> smin_;
std::vector<double> smax_;
};
typedef SimpleFluid2pWrappingProps TwophaseFluid;
typedef Opm::SinglePointUpwindTwoPhase<TwophaseFluid> TransportModel;
using namespace Opm::ImplicitTransportDefault;
typedef NewtonVectorCollection< ::std::vector<double> > NVecColl;
typedef JacobianSystem < struct CSRMatrix, NVecColl > JacSys;
template <class Vector>
class MaxNorm {
public:
static double
norm(const Vector& v) {
return AccumulationNorm <Vector, MaxAbs>::norm(v);
}
};
typedef Opm::ImplicitTransport<TransportModel,
JacSys ,
MaxNorm ,
VectorNegater ,
VectorZero ,
MatrixZero ,
VectorAssign > TransportSolver;
// ----------------- Main program -----------------
int
main(int argc, char** argv)
{
using namespace Opm;
std::cout << "\n================ Test program for incompressible two-phase flow ===============\n\n";
Opm::parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
// Reading various control parameters.
const bool guess_old_solution = param.getDefault("guess_old_solution", false);
const bool use_reorder = param.getDefault("use_reorder", true);
const bool output = param.getDefault("output", true);
std::string output_dir;
int output_interval = 1;
if (output) {
output_dir = param.getDefault("output_dir", std::string("output"));
// Ensure that output dir exists
boost::filesystem::path fpath(output_dir);
try {
create_directories(fpath);
}
catch (...) {
THROW("Creating directories failed: " << fpath);
}
output_interval = param.getDefault("output_interval", output_interval);
}
const int num_transport_substeps = param.getDefault("num_transport_substeps", 1);
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
boost::scoped_ptr<Opm::GridManager> grid;
boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
boost::scoped_ptr<Opm::WellsManager> wells;
boost::scoped_ptr<Opm::RockCompressibility> rock_comp;
Opm::SimulatorTimer simtimer;
Opm::TwophaseState state;
bool check_well_controls = false;
int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::EclipseGridParser deck(deck_filename);
// Grid init
grid.reset(new Opm::GridManager(deck));
// Rock and fluid init
props.reset(new Opm::IncompPropertiesFromDeck(deck, *grid->c_grid()));
// Wells init.
wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability()));
check_well_controls = param.getDefault("check_well_controls", false);
max_well_control_iterations = param.getDefault("max_well_control_iterations", 10);
// Timer init.
if (deck.hasField("TSTEP")) {
simtimer.init(deck);
} else {
simtimer.init(param);
}
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(deck));
// Gravity.
gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
// Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
} else {
// Grid init.
const int nx = param.getDefault("nx", 100);
const int ny = param.getDefault("ny", 100);
const int nz = param.getDefault("nz", 1);
const double dx = param.getDefault("dx", 1.0);
const double dy = param.getDefault("dy", 1.0);
const double dz = param.getDefault("dz", 1.0);
grid.reset(new Opm::GridManager(nx, ny, nz, dx, dy, dz));
// Rock and fluid init.
props.reset(new Opm::IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells));
// Wells init.
wells.reset(new Opm::WellsManager());
// Timer init.
simtimer.init(param);
// Rock compressibility.
rock_comp.reset(new Opm::RockCompressibility(param));
// Gravity.
gravity[2] = param.getDefault("gravity", 0.0);
// Init state variables (saturation and pressure).
initStateBasic(*grid->c_grid(), *props, param, gravity[2], state);
}
// Extra fluid init for transport solver.
TwophaseFluid fluid(*props);
// Warn if gravity but no density difference.
bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0);
if (use_gravity) {
if (props->density()[0] == props->density()[1]) {
std::cout << "**** Warning: nonzero gravity, but zero density difference." << std::endl;
}
}
bool use_segregation_split = false;
bool use_column_solver = false;
bool use_gauss_seidel_gravity = false;
if (use_gravity && use_reorder) {
use_segregation_split = param.getDefault("use_segregation_split", use_segregation_split);
if (use_segregation_split) {
use_column_solver = param.getDefault("use_column_solver", use_column_solver);
if (use_column_solver) {
use_gauss_seidel_gravity = param.getDefault("use_gauss_seidel_gravity", use_gauss_seidel_gravity);
}
}
}
// Check that rock compressibility is not used with solvers that do not handle it.
int nl_pressure_maxiter = 0;
double nl_pressure_residual_tolerance = 0.0;
double nl_pressure_change_tolerance = 0.0;
if (rock_comp->isActive()) {
if (!use_reorder) {
THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
}
nl_pressure_residual_tolerance = param.getDefault("nl_pressure_residual_tolerance", 0.0);
nl_pressure_change_tolerance = param.getDefault("nl_pressure_change_tolerance", 1.0); // In Pascal.
nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10);
}
// Source-related variables init.
int num_cells = grid->c_grid()->number_of_cells;
// Extra rock init.
std::vector<double> porevol;
if (rock_comp->isActive()) {
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
} else {
computePorevolume(*grid->c_grid(), props->porosity(), porevol);
}
double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
// Initialising src
std::vector<double> src(num_cells, 0.0);
if (wells->c_wells()) {
// Do nothing, wells will be the driving force, not source terms.
// Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
} else {
const double default_injection = use_gravity ? 0.0 : 0.1;
const double flow_per_sec = param.getDefault<double>("injected_porevolumes_per_day", default_injection)
*tot_porevol_init/Opm::unit::day;
src[0] = flow_per_sec;
src[num_cells - 1] = -flow_per_sec;
}
TransportSource* tsrc = create_transport_source(2, 2);
double ssrc[] = { 1.0, 0.0 };
double ssink[] = { 0.0, 1.0 };
double zdummy[] = { 0.0, 0.0 };
for (int cell = 0; cell < num_cells; ++cell) {
if (src[cell] > 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssrc, zdummy, tsrc);
} else if (src[cell] < 0.0) {
append_transport_source(cell, 2, 0, src[cell], ssink, zdummy, tsrc);
}
}
std::vector<double> reorder_src = src;
// Boundary conditions.
Opm::FlowBCManager bcs;
if (param.getDefault("use_pside", false)) {
int pside = param.get<int>("pside");
double pside_pressure = param.get<double>("pside_pressure");
bcs.pressureSide(*grid->c_grid(), Opm::FlowBCManager::Side(pside), pside_pressure);
}
// Solvers init.
// Linear solver.
Opm::LinearSolverFactory linsolver(param);
// Pressure solver.
const double *grav = use_gravity ? &gravity[0] : 0;
Opm::IncompTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
nl_pressure_residual_tolerance, nl_pressure_change_tolerance,
nl_pressure_maxiter,
grav, wells->c_wells(), src, bcs.c_bcs());
// Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30);
Opm::TransportModelTwophase reorder_model(*grid->c_grid(), *props, nl_tolerance, nl_maxiter);
if (use_gauss_seidel_gravity) {
reorder_model.initGravity(grav);
}
// Non-reordering solver.
TransportModel model (fluid, *grid->c_grid(), porevol, grav, guess_old_solution);
if (use_gravity) {
model.initGravityTrans(*grid->c_grid(), psolver.getHalfTrans());
}
TransportSolver tsolver(model);
// Column-based gravity segregation solver.
std::vector<std::vector<int> > columns;
if (use_column_solver) {
Opm::extractColumn(*grid->c_grid(), columns);
}
Opm::GravityColumnSolver<TransportModel> colsolver(model, *grid->c_grid(), nl_tolerance, nl_maxiter);
// Control init.
Opm::ImplicitTransportDetails::NRReport rpt;
Opm::ImplicitTransportDetails::NRControl ctrl;
if (!use_reorder || use_segregation_split) {
ctrl.max_it = param.getDefault("max_it", 20);
ctrl.verbosity = param.getDefault("verbosity", 0);
ctrl.max_it_ls = param.getDefault("max_it_ls", 5);
}
// Linear solver init.
using Opm::ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
// The allcells vector is used in calls to computeTotalMobility()
// and computeTotalMobilityOmega().
std::vector<int> allcells(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
allcells[cell] = cell;
}
// Warn if any parameters are unused.
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
// Write parameters used for later reference.
if (output) {
param.writeParam(output_dir + "/simulation.param");
}
// Main simulation loop.
Opm::time::StopWatch pressure_timer;
double ptime = 0.0;
Opm::time::StopWatch transport_timer;
double ttime = 0.0;
Opm::time::StopWatch total_timer;
total_timer.start();
std::cout << "\n\n================ Starting main simulation loop ===============" << std::endl;
double init_satvol[2] = { 0.0 };
double satvol[2] = { 0.0 };
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
double tot_injected[2] = { 0.0 };
double tot_produced[2] = { 0.0 };
Opm::computeSaturatedVol(porevol, state.saturation(), init_satvol);
std::cout << "\nInitial saturations are " << init_satvol[0]/tot_porevol_init
<< " " << init_satvol[1]/tot_porevol_init << std::endl;
Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport;
Opm::WellState well_state;
well_state.init(wells->c_wells(), state);
std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase;
int num_wells = 0;
if (wells->c_wells()) {
num_wells = wells->c_wells()->number_of_wells;
well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0);
wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
}
for (; !simtimer.done(); ++simtimer) {
// Report timestep and (optionally) write state to disk.
simtimer.report(std::cout);
if (output && (simtimer.currentStepNum() % output_interval == 0)) {
outputState(*grid->c_grid(), state, simtimer , output_dir);
}
// Solve pressure.
if (check_well_controls) {
computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows);
}
if (check_well_controls) {
wells->applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
}
bool well_control_passed = !check_well_controls;
int well_control_iteration = 0;
do {
pressure_timer.start();
std::vector<double> initial_pressure = state.pressure();
psolver.solve(simtimer.currentStepLength(), state, well_state);
if (!rock_comp->isActive()) {
// Compute average pressures of previous and last
// step, and total volume.
double av_prev_press = 0.;
double av_press = 0.;
double tot_vol = 0.;
for (int cell = 0; cell < num_cells; ++cell) {
av_prev_press += initial_pressure[cell]*grid->c_grid()->cell_volumes[cell];
av_press += state.pressure()[cell]*grid->c_grid()->cell_volumes[cell];
tot_vol += grid->c_grid()->cell_volumes[cell];
}
// Renormalization constant
const double ren_const = (av_prev_press - av_press)/tot_vol;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += ren_const;
}
for (int well = 0; well < num_wells; ++well) {
well_state.bhp()[well] += ren_const;
}
}
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
ptime += pt;
if (check_well_controls) {
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
fractional_flows,
well_state.perfRates(),
well_resflows_phase);
std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir
well_control_passed = wells->conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
}
if (!well_control_passed) {
std::cout << "Well controls not passed, solving again." << std::endl;
} else {
std::cout << "Well conditions met." << std::endl;
}
}
} while (!well_control_passed);
// Update pore volumes if rock is compressible.
if (rock_comp->isActive()) {
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol);
}
// Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_state.perfRates(), reorder_src);
if (!use_reorder) {
clear_transport_source(tsrc);
for (int cell = 0; cell < num_cells; ++cell) {
if (reorder_src[cell] > 0.0) {
append_transport_source(cell, 2, 0, reorder_src[cell], ssrc, zdummy, tsrc);
} else if (reorder_src[cell] < 0.0) {
append_transport_source(cell, 2, 0, reorder_src[cell], ssink, zdummy, tsrc);
}
}
}
// Solve transport.
transport_timer.start();
double stepsize = simtimer.currentStepLength();
if (num_transport_substeps != 1) {
stepsize /= double(num_transport_substeps);
std::cout << "Making " << num_transport_substeps << " transport substeps." << std::endl;
}
for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
if (use_reorder) {
reorder_model.solve(&state.faceflux()[0], &porevol[0], &reorder_src[0],
stepsize, state.saturation());
Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
if (use_segregation_split) {
if (use_column_solver) {
if (use_gauss_seidel_gravity) {
reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation());
} else {
colsolver.solve(columns, stepsize, state.saturation());
}
} else {
std::vector<double> fluxes = state.faceflux();
std::fill(state.faceflux().begin(), state.faceflux().end(), 0.0);
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
std::cout << rpt;
state.faceflux() = fluxes;
}
}
} else {
tsolver.solve(*grid->c_grid(), tsrc, stepsize, ctrl, state, linsolve, rpt);
std::cout << rpt;
Opm::computeInjectedProduced(*props, state.saturation(), reorder_src, stepsize, injected, produced);
}
}
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
ttime += tt;
// Report volume balances.
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
tot_injected[0] += injected[0];
tot_injected[1] += injected[1];
tot_produced[0] += produced[0];
tot_produced[1] += produced[1];
std::cout.precision(5);
const int width = 18;
std::cout << "\nVolume balance report (all numbers relative to total pore volume).\n";
std::cout << " Saturated volumes: "
<< std::setw(width) << satvol[0]/tot_porevol_init
<< std::setw(width) << satvol[1]/tot_porevol_init << std::endl;
std::cout << " Injected volumes: "
<< std::setw(width) << injected[0]/tot_porevol_init
<< std::setw(width) << injected[1]/tot_porevol_init << std::endl;
std::cout << " Produced volumes: "
<< std::setw(width) << produced[0]/tot_porevol_init
<< std::setw(width) << produced[1]/tot_porevol_init << std::endl;
std::cout << " Total inj volumes: "
<< std::setw(width) << tot_injected[0]/tot_porevol_init
<< std::setw(width) << tot_injected[1]/tot_porevol_init << std::endl;
std::cout << " Total prod volumes: "
<< std::setw(width) << tot_produced[0]/tot_porevol_init
<< std::setw(width) << tot_produced[1]/tot_porevol_init << std::endl;
std::cout << " In-place + prod - inj: "
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init << std::endl;
std::cout << " Init - now - pr + inj: "
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
<< std::endl;
std::cout.precision(8);
watercut.push(simtimer.currentTime() + simtimer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
if (wells->c_wells()) {
wellreport.push(*props, *wells->c_wells(), state.saturation(),
simtimer.currentTime() + simtimer.currentStepLength(),
well_state.bhp(), well_state.perfRates());
}
}
total_timer.stop();
std::cout << "\n\n================ End of simulation ===============\n"
<< "Total time taken: " << total_timer.secsSinceStart()
<< "\n Pressure time: " << ptime
<< "\n Transport time: " << ttime << std::endl;
if (output) {
outputState(*grid->c_grid(), state, simtimer, output_dir);
outputWaterCut(watercut, output_dir);
if (wells->c_wells()) {
outputWellReport(wellreport, output_dir);
}
}
destroy_transport_source(tsrc);
}

10
examples/wells_example.cpp
View File

@@ -3,20 +3,20 @@
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/initState.hpp>
#include <opm/core/simulator/initState.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/GridManager.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/newwells.h>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/wells.h>
#include <opm/core/grid.h>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/pressure/FlowBCManager.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
int main(int argc, char** argv)
{

297
opm/core/GridAdapter.hpp
View File

@@ -1,297 +0,0 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_GRIDADAPTER_HEADER_INCLUDED
#define OPM_GRIDADAPTER_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <stdexcept>
class GridAdapter
{
public:
/// @brief
/// Initialize the grid.
/// @tparam Grid This must conform to the SimpleGrid concept.
/// @param grid The grid object.
template <class Grid>
void init(const Grid& grid)
{
buildTopology(grid);
buildGeometry(grid);
}
UnstructuredGrid* c_grid()
{
return &g_;
}
/// Access the underlying C grid.
const UnstructuredGrid* c_grid() const
{
return &g_;
}
// ------ Forwarding the same interface that init() expects ------
//
// This is only done in order to verify that init() works correctly.
//
enum { dimension = 3 }; // This is actually a hack used for testing (dim is a runtime parameter).
struct Vector
{
explicit Vector(const double* source)
{
for (int i = 0; i < dimension; ++i) {
data[i] = source[i];
}
}
double& operator[] (const int ix)
{
return data[ix];
}
double operator[] (const int ix) const
{
return data[ix];
}
double data[dimension];
};
// Topology
int numCells() const
{
return g_.number_of_cells;
}
int numFaces() const
{
return g_.number_of_faces;
}
int numVertices() const
{
return g_.number_of_nodes;
}
int numCellFaces(int cell) const
{
return cell_facepos_[cell + 1] - cell_facepos_[cell];
}
int cellFace(int cell, int local_index) const
{
return cell_faces_[cell_facepos_[cell] + local_index];
}
int faceCell(int face, int local_index) const
{
return face_cells_[2*face + local_index];
}
int numFaceVertices(int face) const
{
return face_nodepos_[face + 1] - face_nodepos_[face];
}
int faceVertex(int face, int local_index) const
{
return face_nodes_[face_nodepos_[face] + local_index];
}
// Geometry
Vector vertexPosition(int vertex) const
{
return Vector(&node_coordinates_[g_.dimensions*vertex]);
}
double faceArea(int face) const
{
return face_areas_[face];
}
Vector faceCentroid(int face) const
{
return Vector(&face_centroids_[g_.dimensions*face]);
}
Vector faceNormal(int face) const
{
Vector fn(&face_normals_[g_.dimensions*face]);
// We must renormalize since the stored normals are
// 'unit normal * face area'.
double invfa = 1.0 / faceArea(face);
for (int i = 0; i < dimension; ++i) {
fn[i] *= invfa;
}
return fn;
}
double cellVolume(int cell) const
{
return cell_volumes_[cell];
}
Vector cellCentroid(int cell) const
{
return Vector(&cell_centroids_[g_.dimensions*cell]);
}
bool operator==(const GridAdapter& other)
{
return face_nodes_ == other.face_nodes_
&& face_nodepos_ == other.face_nodepos_
&& face_cells_ == other.face_cells_
&& cell_faces_ == other.cell_faces_
&& cell_facepos_ == other.cell_facepos_
&& node_coordinates_ == other.node_coordinates_
&& face_centroids_ == other.face_centroids_
&& face_areas_ == other.face_areas_
&& face_normals_ == other.face_normals_
&& cell_centroids_ == other.cell_centroids_
&& cell_volumes_ == other.cell_volumes_;
}
// make a grid which looks periodic but do not have 2 half faces for each
// periodic boundary
void makeQPeriodic(const std::vector<int>& hf_ind,const std::vector<int>& periodic_cells){
for(int i=0; i<int(hf_ind.size());++i){
//std::array<int,2> cells;
int& cell0=face_cells_[2*cell_faces_[ hf_ind[i] ]+0];
int& cell1=face_cells_[2*cell_faces_[ hf_ind[i] ]+1];
assert(periodic_cells[2*i+1]>=0);
if(periodic_cells[2*i+0] == cell0){
assert(cell1==-1);
cell1=periodic_cells[2*i+1];
}else{
assert(periodic_cells[2*i+0] == cell1);
assert(cell0==-1);
cell0=periodic_cells[2*i+1];
}
}
}
private:
UnstructuredGrid g_;
// Topology storage.
std::vector<int> face_nodes_;
std::vector<int> face_nodepos_;
std::vector<int> face_cells_;
std::vector<int> cell_faces_;
std::vector<int> cell_facepos_;
// Geometry storage.
std::vector<double> node_coordinates_;
std::vector<double> face_centroids_;
std::vector<double> face_areas_;
std::vector<double> face_normals_;
std::vector<double> cell_centroids_;
std::vector<double> cell_volumes_;
/// Build (copy of) topological structure from grid.
template <class Grid>
void buildTopology(const Grid& grid)
{
// Face topology.
int num_cells = grid.numCells();
int num_faces = grid.numFaces();
face_nodepos_.resize(num_faces + 1);
int facenodecount = 0;
for (int f = 0; f < num_faces; ++f) {
face_nodepos_[f] = facenodecount;
facenodecount += grid.numFaceVertices(f);
}
face_nodepos_.back() = facenodecount;
face_nodes_.resize(facenodecount);
for (int f = 0; f < num_faces; ++f) {
for (int local = 0; local < grid.numFaceVertices(f); ++local) {
face_nodes_[face_nodepos_[f] + local] = grid.faceVertex(f, local);
}
}
face_cells_.resize(2*num_faces);
for (int f = 0; f < num_faces; ++f) {
face_cells_[2*f] = grid.faceCell(f, 0);
face_cells_[2*f + 1] = grid.faceCell(f, 1);
}
// Cell topology.
int cellfacecount = 0;
cell_facepos_.resize(num_cells + 1);
for (int c = 0; c < num_cells; ++c) {
cell_facepos_[c] = cellfacecount;
cellfacecount += grid.numCellFaces(c);
}
cell_facepos_.back() = cellfacecount;
cell_faces_.resize(cellfacecount);
for (int c = 0; c < num_cells; ++c) {
for (int local = 0; local < grid.numCellFaces(c); ++local) {
cell_faces_[cell_facepos_[c] + local] = grid.cellFace(c, local);
}
}
// Set C grid members.
g_.dimensions = Grid::dimension;
g_.number_of_cells = grid.numCells();
g_.number_of_faces = grid.numFaces();
g_.number_of_nodes = grid.numVertices();
g_.face_nodes = &face_nodes_[0];
g_.face_nodepos = &face_nodepos_[0];
g_.face_cells = &face_cells_[0];
g_.cell_faces = &cell_faces_[0];
g_.cell_facepos = &cell_facepos_[0];
}
/// Build (copy of) geometric properties of grid.
/// Assumes that buildTopology() has been called.
template <class Grid>
void buildGeometry(const Grid& grid)
{
// Node geometry.
int num_cells = grid.numCells();
int num_nodes = grid.numVertices();
int num_faces = grid.numFaces();
int dim = Grid::dimension;
node_coordinates_.resize(dim*num_nodes);
for (int n = 0; n < num_nodes; ++n) {
for (int dd = 0; dd < dim; ++dd) {
node_coordinates_[dim*n + dd] = grid.vertexPosition(n)[dd];
}
}
// Face geometry.
face_centroids_.resize(dim*num_faces);
face_areas_.resize(num_faces);
face_normals_.resize(dim*num_faces);
for (int f = 0; f < num_faces; ++f) {
face_areas_[f] = grid.faceArea(f);
for (int dd = 0; dd < dim; ++dd) {
face_centroids_[dim*f + dd] = grid.faceCentroid(f)[dd];
face_normals_[dim*f + dd] = grid.faceNormal(f)[dd]*face_areas_[f];
}
}
// Cell geometry.
cell_centroids_.resize(dim*num_cells);
cell_volumes_.resize(num_cells);
for (int c = 0; c < num_cells; ++c) {
cell_volumes_[c] = grid.cellVolume(c);
for (int dd = 0; dd < dim; ++dd) {
cell_centroids_[dim*c + dd] = grid.cellCentroid(c)[dd];
}
}
// Set C grid members.
g_.node_coordinates = &node_coordinates_[0];
g_.face_centroids = &face_centroids_[0];
g_.face_areas = &face_areas_[0];
g_.face_normals = &face_normals_[0];
g_.cell_centroids = &cell_centroids_[0];
g_.cell_volumes = &cell_volumes_[0];
}
};
#endif // OPM_GRIDADAPTER_HEADER_INCLUDED

107
opm/core/fluid/SimpleFluid2p.hpp
View File

@@ -1,107 +0,0 @@
/*===========================================================================
//
// File: SimpleFluid2p.hpp
//
// Created: 2011-09-30 11:38:28+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SIMPLEFLUID2P_HPP_HEADER
#define OPM_SIMPLEFLUID2P_HPP_HEADER
#include <tr1/array>
#include <cmath>
namespace Opm {
template <int n = 2>
class SimpleFluid2p {
public:
SimpleFluid2p(const std::tr1::array<double, 2>& mu,
const std::tr1::array<double, 2>& rho)
: mu_(mu), rho_(rho)
{
}
double density(int p) const { return rho_[p]; }
template <class Sat ,
class Mob ,
class DMob>
void
mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const {
(void) c; // Unused
const double s1 = s[0];
const double s2 = 1 - s1;
if (n == 1) {
mob [ 0] = s1; mob [ 1] = s2;
dmob[0*2 + 0] = 1 ; dmob[1*2 + 1] = 1 ;
} else if (n == 2) {
mob [ 0] = s1 * s1; mob [ 1] = s2 * s2;
dmob[0*2 + 0] = 2 * s1 ; dmob[1*2 + 1] = 2 * s2 ;
} else {
mob [ 0] = pow(s1, double(n));
mob [ 1] = pow(s2, double(n));
dmob[0*2 + 0] = double(n) * pow(s1, double(n) - 1);
dmob[1*2 + 1] = double(n) * pow(s2, double(n) - 1);
}
mob[0] /= mu_[0]; dmob[0*2 + 0] /= mu_[0];
mob[1] /= mu_[1]; dmob[1*2 + 1] /= mu_[1];
dmob[0*2 + 1] = dmob[1*2 + 0] = 0;
}
template <class Sat ,
class Pcap ,
class DPcap>
void
pc(int c, const Sat& s, Pcap& pcap, DPcap& dpcap) const {
(void) c; // Unused
(void) s; // Unused
pcap = 0.0;
dpcap = 0.0;
}
double s_min(int c) const { (void) c; return 0.0; }
double s_max(int c) const { (void) c; return 1.0; }
private:
std::tr1::array<double, 2> mu_ ;
std::tr1::array<double, 2> rho_;
};
}
#endif /* OPM_SIMPLEFLUID2P_HPP_HEADER */

0
opm/core/utility/ColumnExtract.hpp → opm/core/grid/ColumnExtract.hpp
View File

12
opm/core/GridManager.cpp → opm/core/grid/GridManager.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/GridManager.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
@@ -63,13 +63,19 @@ namespace Opm
/// Construct a 2d cartesian grid with cells of unit size.
GridManager::GridManager(int nx, int ny)
{
ug_ = create_grid_cart2d(nx, ny);
ug_ = create_grid_cart2d(nx, ny, 1.0, 1.0);
if (!ug_) {
THROW("Failed to construct grid.");
}
}
GridManager::GridManager(int nx, int ny,double dx, double dy)
{
ug_ = create_grid_cart2d(nx, ny, dx, dy);
if (!ug_) {
THROW("Failed to construct grid.");
}
}
/// Construct a 3d cartesian grid with cells of unit size.

3
opm/core/GridManager.hpp → opm/core/grid/GridManager.hpp
View File

@@ -47,6 +47,9 @@ namespace Opm
/// Construct a 2d cartesian grid with cells of unit size.
GridManager(int nx, int ny);
/// Construct a 2d cartesian grid with cells of size [dx, dy].
GridManager(int nx, int ny, double dx, double dy);
/// Construct a 3d cartesian grid with cells of unit size.
GridManager(int nx, int ny, int nz);

6
opm/core/grid/cart_grid.c
View File

@@ -91,7 +91,7 @@ static void fill_cart_geometry_2d(struct UnstructuredGrid *G,
const double *y);
struct UnstructuredGrid*
create_grid_cart2d(int nx, int ny)
create_grid_cart2d(int nx, int ny, double dx, double dy)
{
int i;
double *x, *y;
@@ -104,8 +104,8 @@ create_grid_cart2d(int nx, int ny)
G = NULL;
} else {
for (i = 0; i < nx + 1; i++) { x[i] = i; }
for (i = 0; i < ny + 1; i++) { y[i] = i; }
for (i = 0; i < nx + 1; i++) { x[i] = i*dx; }
for (i = 0; i < ny + 1; i++) { y[i] = i*dy; }
G = create_grid_tensor2d(nx, ny, x, y);
}

6
opm/core/grid/cart_grid.h
View File

@@ -44,7 +44,7 @@
extern "C" {
#endif
struct UnstructuredGrid;
/**
* Form geometrically Cartesian grid in two space dimensions with unit-sized
@@ -52,12 +52,14 @@ struct UnstructuredGrid;
*
* @param[in] nx Number of cells in @c x direction.
* @param[in] ny Number of cells in @c y direction.
* @param[in] dx Length, in meters, of each cell's @c x extent.
* @param[in] dy Length, in meters, of each cell's @c y extent.
*
* @return Fully formed grid structure containing valid geometric primitives.
* Must be destroyed using function destroy_grid().
*/
struct UnstructuredGrid *
create_grid_cart2d(int nx, int ny);
create_grid_cart2d(int nx, int ny,double dx, double dy);
/**

92
opm/core/grid/cpgpreprocess/geometry.c
View File

@@ -2,6 +2,7 @@
* Copyright 2010 (c) SINTEF ICT, Applied Mathematics.
* Jostein R. Natvig <Jostein.R.Natvig at sintef.no>
*/
#include <omp.h>
#include <math.h>
#include <stdio.h>
#include "geometry.h"
@@ -47,11 +48,20 @@ compute_face_geometry_3d(double *coords, int nfaces,
double cface[3] = {0};
double n[3] = {0};
const double twothirds = 0.666666666666666666666666666667;
double twothirds = 0.666666666666666666666666666667;
double a;
int num_face_nodes;
double area;
/*#pragma omp parallel for */
/*#pragma omp parallel for shared(fnormals,fcentroids,fareas)*/
#pragma omp parallel for default(none) \
private(f,x,u,v,w,i,k,node,cface,n,a,num_face_nodes,area) \
shared(fnormals,fcentroids,fareas \
,coords, nfaces, nodepos, facenodes) \
firstprivate(ndims, twothirds)
for (f=0; f<nfaces; ++f)
{
int num_face_nodes;
double area = 0.0;
for(i=0; i<ndims; ++i) x[i] = 0.0;
for(i=0; i<ndims; ++i) n[i] = 0.0;
for(i=0; i<ndims; ++i) cface[i] = 0.0;
@@ -71,11 +81,11 @@ compute_face_geometry_3d(double *coords, int nfaces,
node = facenodes[nodepos[f+1]-1];
for(i=0; i<ndims; ++i) u[i] = coords[3*node+i] - x[i];
area=0.0;
/* Compute triangular contrib. to face normal and face centroid*/
for(k=nodepos[f]; k<nodepos[f+1]; ++k)
{
double a;
node = facenodes[k];
for (i=0; i<ndims; ++i) v[i] = coords[3*node+i] - x[i];
@@ -83,11 +93,11 @@ compute_face_geometry_3d(double *coords, int nfaces,
cross(u,v,w);
a = 0.5*norm(w);
area += a;
if(!(a>0))
/* if(!(a>0))
{
fprintf(stderr, "Internal error in compute_face_geometry.");
}
*/
/* face normal */
for (i=0; i<ndims; ++i) n[i] += w[i];
@@ -221,17 +231,19 @@ compute_cell_geometry_3d(double *coords,
double xcell[3];
double ccell[3];
double cface[3] = {0};
const double twothirds = 0.666666666666666666666666666667;
int ndigits;
ndigits = ((int) (log(ncells) / log(10.0))) + 1;
int num_faces;
double volume;
double tet_volume, subnormal_sign;
double twothirds = 0.666666666666666666666666666667;
#pragma omp parallel for default(none) \
private(i,k,f,c,face,node,x,u,v,w,xcell \
,ccell ,cface,num_faces,volume, tet_volume, subnormal_sign) \
shared(coords,nodepos,facenodes,neighbors, \
fnormals,fcentroids,facepos,cellfaces,ccentroids,cvolumes) \
firstprivate(ncells,ndims,twothirds)
for (c=0; c<ncells; ++c)
{
int num_faces;
double volume = 0.0;
for(i=0; i<ndims; ++i) xcell[i] = 0.0;
for(i=0; i<ndims; ++i) ccell[i] = 0.0;
@@ -255,6 +267,7 @@ compute_cell_geometry_3d(double *coords,
* For all faces, add tetrahedron's volume and centroid to
* 'cvolume' and 'ccentroid'.
*/
volume=0.0;
for(f=facepos[c]; f<facepos[c+1]; ++f)
{
int num_face_nodes;
@@ -279,38 +292,38 @@ compute_cell_geometry_3d(double *coords,
node = facenodes[nodepos[face+1]-1];
for(i=0; i<ndims; ++i) u[i] = coords[3*node+i] - x[i];
/* Compute triangular contributions to face normal and face centroid */
for(k=nodepos[face]; k<nodepos[face+1]; ++k)
{
double tet_volume, subnormal_sign;
node = facenodes[k];
for (i=0; i<ndims; ++i) v[i] = coords[3*node+i] - x[i];
cross(u,v,w);
tet_volume = 0;
tet_volume = 0.0;
for(i=0; i<ndims; ++i){
tet_volume += w[i] * (x[i] - xcell[i]);
}
tet_volume += w[i]*(x[i]-xcell[i]);
}
tet_volume *= 0.5 / 3;
subnormal_sign=0.0;
for(i=0; i<ndims; ++i){
subnormal_sign += w[i]*fnormals[3*face+i];
}
subnormal_sign=0.0;
for(i=0; i<ndims; ++i){
subnormal_sign += w[i]*fnormals[3*face+i];
}
if (subnormal_sign < 0.0) {
tet_volume = - tet_volume;
}
if (neighbors[2*face + 0] != c) {
tet_volume = - tet_volume;
}
volume += tet_volume;
if(subnormal_sign < 0.0){
tet_volume =- tet_volume;
}
if(!(neighbors[2*face+0]==c)){
tet_volume = -tet_volume;
}
volume += tet_volume;
/* face centroid of triangle */
for (i=0; i<ndims; ++i) cface[i] = (x[i]+twothirds*0.5*(u[i]+v[i]));
for (i=0; i<ndims; ++i) cface[i] = (x[i]+(twothirds)*0.5*(u[i]+v[i]));
/* Cell centroid */
for (i=0; i<ndims; ++i) ccell[i] += tet_volume * 3/4.0*(cface[i] - xcell[i]);
@@ -320,13 +333,6 @@ compute_cell_geometry_3d(double *coords,
for (i=0; i<ndims; ++i) u[i] = v[i];
}
}
if (! (volume > 0.0)) {
fprintf(stderr,
"Internal error in mex_compute_geometry(%*d): "
"negative volume\n", ndigits, c);
}
for (i=0; i<ndims; ++i) ccentroids[3*c+i] = xcell[i] + ccell[i]/volume;
cvolumes[c] = volume;
}

0
opm/core/grid.c → opm/core/grid/grid.c
View File

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

@@ -27,10 +27,10 @@
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <algorithm>
#include <cmath>

328
opm/core/pressure/HybridPressureSolver.hpp
View File

@@ -1,328 +0,0 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED
#define OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED
#include <opm/core/pressure/fsh.h>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/pressure/mimetic/mimetic.h>
#include <opm/core/GridAdapter.hpp>
#include <stdexcept>
/// @brief
/// Encapsulates the ifsh (= incompressible flow solver hybrid) solver modules.
class HybridPressureSolver
{
public:
/// @brief
/// Default constructor, does nothing.
HybridPressureSolver()
: state_(Uninitialized), data_(0)
{
}
/// @brief
/// Destructor.
~HybridPressureSolver()
{
fsh_destroy(data_);
}
/// @brief
/// Initialize the solver's structures for a given grid (at some point also well pattern).
/// @tparam Grid This must conform to the SimpleGrid concept.
/// @param grid The grid object.
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
template <class Grid>
void init(const Grid& grid, const double* perm, const double* gravity)
{
// Build C grid structure.
grid_.init(grid);
// Checking if grids are properly initialized. Move this code to a test somewhere.
// GridAdapter grid2;
// grid2.init(grid_);
// if (grid2 == grid_) {
// std::cout << "Grids are equal." << std::endl;
// } else {
// std::cout << "Grids are NOT equal." << std::endl;
// }
// Build (empty for now) C well structure.
well_t* w = 0;
// Initialize ifsh data.
data_ = ifsh_construct(grid_.c_grid(), w);
if (!data_) {
throw std::runtime_error("Failed to initialize ifsh solver.");
}
// Compute inner products, gravity contributions.
int num_cells = grid.numCells();
int ngconn = grid_.c_grid()->cell_facepos[num_cells];
gpress_.clear();
gpress_.resize(ngconn, 0.0);
int ngconn2 = data_->sum_ngconn2;
Binv_.resize(ngconn2);
ncf_.resize(num_cells);
typename Grid::Vector grav;
std::copy(gravity, gravity + Grid::dimension, &grav[0]);
int count = 0;
for (int cell = 0; cell < num_cells; ++cell) {
int num_local_faces = grid.numCellFaces(cell);
ncf_[cell] = num_local_faces;
typename Grid::Vector cc = grid.cellCentroid(cell);
for (int local_ix = 0; local_ix < num_local_faces; ++local_ix) {
int face = grid.cellFace(cell, local_ix);
typename Grid::Vector fc = grid.faceCentroid(face);
gpress_[count++] = grav*(fc - cc);
}
}
assert(count == ngconn);
UnstructuredGrid* g = grid_.c_grid();
mim_ip_simple_all(g->number_of_cells, g->dimensions,
data_->max_ngconn,
g->cell_facepos, g->cell_faces,
g->face_cells, g->face_centroids,
g->face_normals, g->face_areas,
g->cell_centroids, g->cell_volumes,
const_cast<double*>(perm), &Binv_[0]);
state_ = Initialized;
}
enum FlowBCTypes { FBC_UNSET, FBC_PRESSURE, FBC_FLUX };
/// @brief
/// Assemble the sparse system.
/// You must call init() prior to calling assemble().
/// @param sources Source terms, one per cell. Positive numbers
/// are sources, negative are sinks.
/// @param total_mobilities Scalar total mobilities, one per cell.
/// @param omegas Gravity term, one per cell. In a multi-phase
/// flow setting this is equal to
/// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f]
/// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a
/// phase density and \f$\lambda_t\f$ is the total mobility.
void assemble(const std::vector<double>& sources,
const std::vector<double>& total_mobilities,
const std::vector<double>& omegas,
const std::vector<FlowBCTypes>& bctypes,
const std::vector<double>& bcvalues)
{
if (state_ == Uninitialized) {
throw std::runtime_error("Error in HybridPressureSolver::assemble(): You must call init() prior to calling assemble().");
}
// Boundary conditions.
assert (bctypes.size() ==
static_cast<std::vector<FlowBCTypes>::size_type>(grid_.numFaces()));
FlowBoundaryConditions *bc = gather_boundary_conditions(bctypes, bcvalues);
// Source terms from user.
double* src = const_cast<double*>(&sources[0]); // Ugly? Yes. Safe? I think so.
// All well related things are zero.
well_control_t* wctrl = 0;
double* WI = 0;
double* wdp = 0;
// Scale inner products and gravity terms by saturation-dependent factors.
UnstructuredGrid* g = grid_.c_grid();
Binv_mobilityweighted_.resize(Binv_.size());
mim_ip_mobility_update(g->number_of_cells, g->cell_facepos, &total_mobilities[0],
&Binv_[0], &Binv_mobilityweighted_[0]);
gpress_omegaweighted_.resize(gpress_.size());
mim_ip_density_update(g->number_of_cells, g->cell_facepos, &omegas[0],
&gpress_[0], &gpress_omegaweighted_[0]);
// Zero the linalg structures.
csrmatrix_zero(data_->A);
for (std::size_t i = 0; i < data_->A->m; i++) {
data_->b[i] = 0.0;
}
// Assemble the embedded linear system.
ifsh_assemble(bc, src, &Binv_mobilityweighted_[0], &gpress_omegaweighted_[0],
wctrl, WI, wdp, data_);
state_ = Assembled;
flow_conditions_destroy(bc);
}
/// Encapsulate a sparse linear system in CSR format.
struct LinearSystem
{
int n;
int nnz;
int* ia;
int* ja;
double* sa;
double* b;
double* x;
};
/// @brief
/// Access the linear system assembled.
/// You must call assemble() prior to calling linearSystem().
/// @param[out] s The linear system encapsulation to modify.
/// After this call, s will point to linear system structures
/// that are owned and allocated internally.
void linearSystem(LinearSystem& s)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in HybridPressureSolver::linearSystem(): "
"You must call assemble() prior to calling linearSystem().");
}
s.n = data_->A->m;
s.nnz = data_->A->nnz;
s.ia = data_->A->ia;
s.ja = data_->A->ja;
s.sa = data_->A->sa;
s.b = data_->b;
s.x = data_->x;
}
/// @brief
/// Compute cell pressures and face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// computePressuresAndFluxes().
/// @param[out] cell_pressures Cell pressure values.
/// @param[out] face_areas Face flux values.
void computePressuresAndFluxes(std::vector<double>& cell_pressures,
std::vector<double>& face_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in HybridPressureSolver::computePressuresAndFluxes(): "
"You must call assemble() (and solve the linear system) "
"prior to calling computePressuresAndFluxes().");
}
int num_cells = grid_.c_grid()->number_of_cells;
int num_faces = grid_.c_grid()->number_of_faces;
cell_pressures.clear();
cell_pressures.resize(num_cells, 0.0);
face_fluxes.clear();
face_fluxes.resize(num_faces, 0.0);
fsh_press_flux(grid_.c_grid(), &Binv_mobilityweighted_[0], &gpress_omegaweighted_[0],
data_, &cell_pressures[0], &face_fluxes[0], 0, 0);
}
/// @brief
/// Compute cell fluxes from face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// faceFluxToCellFlux().
/// @param face_fluxes
/// @param face_areas Face flux values (usually output from computePressuresAndFluxes()).
/// @param[out] cell_fluxes Cell-wise flux values.
/// They are given in cell order, and for each cell there is
/// one value for each adjacent face (in the same order as the
/// cell-face topology of the grid). Positive values represent
/// fluxes out of the cell.
void faceFluxToCellFlux(const std::vector<double>& face_fluxes,
std::vector<double>& cell_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in HybridPressureSolver::faceFluxToCellFlux(): "
"You must call assemble() (and solve the linear system) "
"prior to calling faceFluxToCellFlux().");
}
const UnstructuredGrid& g = *(grid_.c_grid());
int num_cells = g.number_of_cells;
cell_fluxes.resize(g.cell_facepos[num_cells]);
for (int cell = 0; cell < num_cells; ++cell) {
for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
int face = g.cell_faces[hface];
bool pos = (g.face_cells[2*face] == cell);
cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
}
}
}
/// @brief
/// Access the number of connections (faces) per cell. Deprecated, will be removed.
const std::vector<int>& numCellFaces()
{
return ncf_;
}
private:
// Disabling copy and assigment for now.
HybridPressureSolver(const HybridPressureSolver&);
HybridPressureSolver& operator=(const HybridPressureSolver&);
enum State { Uninitialized, Initialized, Assembled };
State state_;
// Solver data.
fsh_data* data_;
// Grid.
GridAdapter grid_;
// Number of faces per cell.
std::vector<int> ncf_;
// B^{-1} storage.
std::vector<double> Binv_;
std::vector<double> Binv_mobilityweighted_;
// Gravity contributions.
std::vector<double> gpress_;
std::vector<double> gpress_omegaweighted_;
FlowBoundaryConditions*
gather_boundary_conditions(const std::vector<FlowBCTypes>& bctypes ,
const std::vector<double>& bcvalues)
{
FlowBoundaryConditions* fbc = flow_conditions_construct(0);
int ok = fbc != 0;
std::vector<FlowBCTypes>::size_type i;
for (i = 0; ok && (i < bctypes.size()); ++i) {
if (bctypes[ i ] == FBC_PRESSURE) {
ok = flow_conditions_append(BC_PRESSURE,
static_cast<int>(i),
bcvalues[ i ],
fbc);
}
else if (bctypes[ i ] == FBC_FLUX) {
ok = flow_conditions_append(BC_FLUX_TOTVOL,
static_cast<int>(i),
bcvalues[ i ],
fbc);
}
}
return fbc;
}
}; // class HybridPressureSolver
#endif // OPM_HYBRIDPRESSURESOLVER_HEADER_INCLUDED

6
opm/core/pressure/IncompTpfa.cpp
View File

@@ -19,8 +19,8 @@
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/mimetic/mimetic.h>
@@ -31,7 +31,7 @@
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <iomanip>
#include <cmath>
#include <algorithm>

527
opm/core/pressure/TPFACompressiblePressureSolver.hpp
View File

@@ -1,527 +0,0 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED
#define OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED
#include <opm/core/pressure/tpfa/cfs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/well.h>
#include <opm/core/pressure/tpfa/compr_quant.h>
#include <opm/core/GridAdapter.hpp>
#include <stdexcept>
/// @brief
/// Encapsulates the cfs_tpfa (= compressible flow solver
/// two-point flux approximation) solver modules.
class TPFACompressiblePressureSolver
{
public:
/// @brief
/// Default constructor, does nothing.
TPFACompressiblePressureSolver()
: state_(Uninitialized), data_(0), bc_(0)
{
wells_.number_of_wells = 0;
}
/// @brief
/// Destructor.
~TPFACompressiblePressureSolver()
{
flow_conditions_destroy(bc_);
cfs_tpfa_destroy(data_);
}
/// @brief
/// Initialize the solver's structures for a given grid, for well setup also call initWells().
/// @tparam Grid This must conform to the SimpleGrid concept.
/// @tparam Wells This must conform to the SimpleWells concept.
/// @param grid The grid object.
/// @param wells Well specifications.
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
/// @param perm Porosity by cell.
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
template <class Grid, class Wells>
void init(const Grid& grid, const Wells& wells, const double* perm, const double* porosity,
const typename Grid::Vector& gravity)
{
initWells(wells);
init(grid, perm, porosity, gravity);
}
/// @brief
/// Initialize the solver's structures for a given grid, for well setup also call initWells().
/// @tparam Grid This must conform to the SimpleGrid concept.
/// @param grid The grid object.
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
template <class Grid>
void init(const Grid& grid, const double* perm, const double* porosity, const typename Grid::Vector& gravity)
{
// Build C grid structure.
grid_.init(grid);
// Initialize data.
int num_phases = 3;
well_t* w = 0;
if (wells_.number_of_wells != 0) {
w = &wells_;
}
data_ = cfs_tpfa_construct(grid_.c_grid(), w, num_phases);
if (!data_) {
throw std::runtime_error("Failed to initialize cfs_tpfa solver.");
}
// Compute half-transmissibilities
int num_cells = grid.numCells();
int ngconn = grid_.c_grid()->cell_facepos[num_cells];
ncf_.resize(num_cells);
for (int cell = 0; cell < num_cells; ++cell) {
int num_local_faces = grid.numCellFaces(cell);
ncf_[cell] = num_local_faces;
}
htrans_.resize(ngconn);
tpfa_htrans_compute(grid_.c_grid(), perm, &htrans_[0]);
// Compute transmissibilities.
trans_.resize(grid_.numFaces());
tpfa_trans_compute(grid_.c_grid(), &htrans_[0], &trans_[0]);
// Compute pore volumes.
porevol_.resize(num_cells);
for (int i = 0; i < num_cells; ++i) {
porevol_[i] = porosity[i]*grid.cellVolume(i);
}
// Set gravity.
if (Grid::dimension != 3) {
throw std::logic_error("Only 3 dimensions supported currently.");
}
std::copy(gravity.begin(), gravity.end(), gravity_);
state_ = Initialized;
}
/// Boundary condition types.
enum FlowBCTypes { FBC_UNSET = BC_NOFLOW ,
FBC_PRESSURE = BC_PRESSURE ,
FBC_FLUX = BC_FLUX_TOTVOL };
/// @brief
/// Assemble the sparse system.
/// You must call init() prior to calling assemble().
/// @param sources Source terms, one per cell. Positive numbers
/// are sources, negative are sinks.
/// @param total_mobilities Scalar total mobilities, one per cell.
/// @param omegas Gravity term, one per cell. In a multi-phase
/// flow setting this is equal to
/// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f]
/// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a
/// phase density and \f$\lambda_t\f$ is the total mobility.
void assemble(const double* sources,
const FlowBCTypes* bctypes,
const double* bcvalues,
const double dt,
const double* totcompr,
const double* voldiscr,
const double* cellA, // num phases^2 * num cells, fortran ordering!
const double* faceA, // num phases^2 * num faces, fortran ordering!
const double* wellperfA,
const double* phasemobf,
const double* phasemobwellperf,
const double* cell_pressure,
const double* gravcapf,
const double* wellperf_gpot,
const double* surf_dens)
{
if (state_ == Uninitialized) {
throw std::runtime_error("Error in TPFACompressiblePressureSolver::assemble(): You must call init() prior to calling assemble().");
}
UnstructuredGrid* g = grid_.c_grid();
// Boundary conditions.
gather_boundary_conditions(bctypes, bcvalues);
// Source terms from user.
double* src = const_cast<double*>(sources); // Ugly? Yes. Safe? I think so.
// Wells.
well_t* wells = NULL;
well_control_t* wctrl = NULL;
struct completion_data* wcompl = NULL;
if (wells_.number_of_wells != 0) {
wells = &wells_;
wctrl = &wctrl_;
wcompl = &wcompl_;
// The next objects already have the correct sizes.
std::copy(wellperf_gpot, wellperf_gpot + well_gpot_storage_.size(), well_gpot_storage_.begin());
std::copy(wellperfA, wellperfA + well_A_storage_.size(), well_A_storage_.begin());
std::copy(phasemobwellperf, phasemobwellperf + well_phasemob_storage_.size(), well_phasemob_storage_.begin());
}
// Assemble the embedded linear system.
compr_quantities cq = { 3 ,
const_cast<double *>(totcompr ) ,
const_cast<double *>(voldiscr ) ,
const_cast<double *>(cellA ) ,
const_cast<double *>(faceA ) ,
const_cast<double *>(phasemobf) };
// Call the assembly routine. After this, linearSystem() may be called.
cfs_tpfa_assemble(g, dt, wells, bc_, src,
&cq, &trans_[0], gravcapf,
wctrl, wcompl,
cell_pressure, &porevol_[0],
data_);
phasemobf_.assign(phasemobf, phasemobf + grid_.numFaces()*3);
gravcapf_.assign(gravcapf, gravcapf + grid_.numFaces()*3);
state_ = Assembled;
}
/// Encapsulate a sparse linear system in CSR format.
struct LinearSystem
{
int n;
int nnz;
int* ia;
int* ja;
double* sa;
double* b;
double* x;
};
/// @brief
/// Access the linear system assembled.
/// You must call assemble() prior to calling linearSystem().
/// @param[out] s The linear system encapsulation to modify.
/// After this call, s will point to linear system structures
/// that are owned and allocated internally.
void linearSystem(LinearSystem& s)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TPFACompressiblePressureSolver::linearSystem(): "
"You must call assemble() prior to calling linearSystem().");
}
s.n = data_->A->m;
s.nnz = data_->A->nnz;
s.ia = data_->A->ia;
s.ja = data_->A->ja;
s.sa = data_->A->sa;
s.b = data_->b;
s.x = data_->x;
}
/// @brief
/// Compute cell pressures and face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// computePressuresAndFluxes().
/// @param[out] cell_pressures Cell pressure values.
/// @param[out] face_areas Face flux values.
void computePressuresAndFluxes(std::vector<double>& cell_pressures,
std::vector<double>& face_pressures,
std::vector<double>& face_fluxes,
std::vector<double>& well_pressures,
std::vector<double>& well_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TPFACompressiblePressureSolver::computePressuresAndFluxes(): "
"You must call assemble() (and solve the linear system) "
"prior to calling computePressuresAndFluxes().");
}
int num_cells = grid_.c_grid()->number_of_cells;
int num_faces = grid_.c_grid()->number_of_faces;
cell_pressures.clear();
cell_pressures.resize(num_cells, 0.0);
face_pressures.clear();
face_pressures.resize(num_faces, 0.0);
face_fluxes.clear();
face_fluxes.resize(num_faces, 0.0);
int np = 3; // Number of phases.
// Wells.
well_t* wells = NULL;
struct completion_data* wcompl = NULL;
double* wpress = 0;
double* wflux = 0;
if (wells_.number_of_wells != 0) {
wells = &wells_;
wcompl = &wcompl_;
well_pressures.resize(wells_.number_of_wells);
well_fluxes.resize(well_cells_storage_.size());
wpress = &well_pressures[0];
wflux = &well_fluxes[0];
}
cfs_tpfa_press_flux(grid_.c_grid(),
bc_, wells,
np, &trans_[0], &phasemobf_[0], &gravcapf_[0],
wcompl,
data_, &cell_pressures[0], &face_fluxes[0],
wpress, wflux);
cfs_tpfa_fpress(grid_.c_grid(), bc_, np, &htrans_[0],
&phasemobf_[0], &gravcapf_[0], data_, &cell_pressures[0],
&face_fluxes[0], &face_pressures[0]);
}
/// @brief
/// Explicit IMPES time step limit.
double explicitTimestepLimit(const double* faceA, // num phases^2 * num faces, fortran ordering!
const double* phasemobf,
const double* phasemobf_deriv,
const double* surf_dens)
{
compr_quantities cq = { 3, // nphases
0, // totcompr
0, // voldiscr
0, // Ac
const_cast<double *>(faceA) ,
const_cast<double *>(phasemobf) };
return cfs_tpfa_impes_maxtime(grid_.c_grid(), &cq, &trans_[0], &porevol_[0], data_,
phasemobf_deriv, surf_dens, gravity_);
}
/// @brief
/// Explicit IMPES transport.
void explicitTransport(const double dt,
double* cell_surfvols)
{
int np = 3; // Number of phases.
well_t* wells = NULL;
if (wells_.number_of_wells != 0) {
wells = &wells_;
}
cfs_tpfa_expl_mass_transport(grid_.c_grid(), wells, &wcompl_,
np, dt, &porevol_[0],
data_, cell_surfvols);
}
/// @brief
/// Compute cell fluxes from face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// faceFluxToCellFlux().
/// @param face_fluxes
/// @param face_areas Face flux values (usually output from computePressuresAndFluxes()).
/// @param[out] cell_fluxes Cell-wise flux values.
/// They are given in cell order, and for each cell there is
/// one value for each adjacent face (in the same order as the
/// cell-face topology of the grid). Positive values represent
/// fluxes out of the cell.
void faceFluxToCellFlux(const std::vector<double>& face_fluxes,
std::vector<double>& cell_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TPFACompressiblePressureSolver::faceFluxToCellFlux(): "
"You must call assemble() (and solve the linear system) "
"prior to calling faceFluxToCellFlux().");
}
const UnstructuredGrid& g = *(grid_.c_grid());
int num_cells = g.number_of_cells;
cell_fluxes.resize(g.cell_facepos[num_cells]);
for (int cell = 0; cell < num_cells; ++cell) {
for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
int face = g.cell_faces[hface];
bool pos = (g.face_cells[2*face] == cell);
cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
}
}
}
/// @brief
/// Access the number of connections (faces) per cell. Deprecated, will be removed.
const std::vector<int>& numCellFaces() const
{
return ncf_;
}
const std::vector<double>& faceTransmissibilities() const
{
return trans_;
}
private:
// Disabling copy and assigment for now.
TPFACompressiblePressureSolver(const TPFACompressiblePressureSolver&);
TPFACompressiblePressureSolver& operator=(const TPFACompressiblePressureSolver&);
enum State { Uninitialized, Initialized, Assembled };
State state_;
// Solver data.
cfs_tpfa_data* data_;
// Grid.
GridAdapter grid_;
// Number of faces per cell.
std::vector<int> ncf_;
// Transmissibility storage.
std::vector<double> htrans_;
std::vector<double> trans_;
// Pore volumes.
std::vector<double> porevol_;
// Phase mobilities per face.
std::vector<double> phasemobf_;
// Gravity and capillary contributions (per face).
std::vector<double> gravcapf_;
// Gravity
double gravity_[3];
// Boundary conditions.
FlowBoundaryConditions *bc_;
// Well data
well_t wells_;
std::vector<int> well_connpos_storage_;
std::vector<int> well_cells_storage_;
well_control_t wctrl_;
std::vector<well_type> wctrl_type_storage_;
std::vector<well_control> wctrl_ctrl_storage_;
std::vector<double> wctrl_target_storage_;
struct completion_data wcompl_;
std::vector<double> well_prodind_storage_;
std::vector<double> well_gpot_storage_;
std::vector<double> well_A_storage_;
std::vector<double> well_phasemob_storage_;
/// @brief
/// Initialize wells in solver structure.
/// @tparam Wells
/// This must conform to the SimpleWells concept.
/// @param w
/// The well object.
template <class Wells>
void initWells(const Wells& w)
{
int num_wells = w.numWells();
if (num_wells == 0) {
wells_.number_of_wells = 0;
return;
}
wctrl_type_storage_.resize(num_wells);
wctrl_ctrl_storage_.resize(num_wells);
wctrl_target_storage_.resize(num_wells);
for (int i = 0; i < num_wells; ++i) {
wctrl_type_storage_[i] = (w.type(i) == Wells::Injector) ? INJECTOR : PRODUCER;
wctrl_ctrl_storage_[i] = (w.control(i) == Wells::Rate) ? RATE : BHP;
wctrl_target_storage_[i] = w.target(i);
int num_perf = w.numPerforations(i);
well_connpos_storage_.push_back(well_cells_storage_.size());
for (int j = 0; j < num_perf; ++j) {
well_cells_storage_.push_back(w.wellCell(i, j));
well_prodind_storage_.push_back(w.wellIndex(i, j));
}
}
well_connpos_storage_.push_back(well_cells_storage_.size());
int tot_num_perf = well_prodind_storage_.size();
well_gpot_storage_.resize(tot_num_perf*3);
well_A_storage_.resize(3*3*tot_num_perf);
well_phasemob_storage_.resize(3*tot_num_perf);
// Setup 'wells_'
wells_.number_of_wells = num_wells;
wells_.well_connpos = &well_connpos_storage_[0];
wells_.well_cells = &well_cells_storage_[0];
// Setup 'wctrl_'
wctrl_.type = &wctrl_type_storage_[0];
wctrl_.ctrl = &wctrl_ctrl_storage_[0];
wctrl_.target = &wctrl_target_storage_[0];
// Setup 'wcompl_'
wcompl_.WI = &well_prodind_storage_[0];
wcompl_.gpot = &well_gpot_storage_[0];
wcompl_.A = &well_A_storage_[0];
wcompl_.phasemob = &well_phasemob_storage_[0];
}
void
gather_boundary_conditions(const FlowBCTypes* bctypes ,
const double* bcvalues)
{
if (bc_ == 0) {
bc_ = flow_conditions_construct(0);
}
else {
flow_conditions_clear(bc_);
}
int ok = bc_ != 0;
for (std::size_t i = 0, nf = grid_.numFaces(); ok && (i < nf); ++i) {
if (bctypes[ i ] == FBC_PRESSURE) {
ok = flow_conditions_append(BC_PRESSURE,
static_cast<int>(i),
bcvalues[ i ],
bc_);
}
else if (bctypes[ i ] == FBC_FLUX) {
ok = flow_conditions_append(BC_FLUX_TOTVOL,
static_cast<int>(i),
bcvalues[ i ],
bc_);
}
}
if (! ok) {
flow_conditions_destroy(bc_);
bc_ = 0;
}
}
}; // class TPFACompressiblePressureSolver
#endif // OPM_TPFACOMPRESSIBLEPRESSURESOLVER_HEADER_INCLUDED

294
opm/core/pressure/TPFAPressureSolver.hpp
View File

@@ -1,294 +0,0 @@
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED
#define OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/pressure/mimetic/mimetic.h> // for updating gpress
#include <opm/core/GridAdapter.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <stdexcept>
#include <cassert>
/// @brief
/// Encapsulates the ifs_tpfa (= incompressible flow solver
/// two-point flux approximation) solver modules.
class TPFAPressureSolver
{
public:
/// @brief
/// Default constructor, does nothing.
TPFAPressureSolver()
: state_(Uninitialized), data_(0)
{
}
/// @brief
/// Destructor.
~TPFAPressureSolver()
{
ifs_tpfa_destroy(data_);
}
/// @brief
/// Initialize the solver's structures for a given grid (at some point also well pattern).
/// @tparam Grid This must conform to the SimpleGrid concept.
/// @param grid The grid object.
/// @param perm Permeability. It should contain dim*dim entries (a full tensor) for each cell.
/// @param gravity Array containing gravity acceleration vector. It should contain dim entries.
template <class Grid>
void init(const Grid& grid, const double* perm, const double* gravity)
{
// Build C grid structure.
grid_.init(grid);
// Build (empty for now) C well structure.
// well_t* w = 0;
// Initialize data.
data_ = ifs_tpfa_construct(grid_.c_grid(), 0);
if (!data_) {
throw std::runtime_error("Failed to initialize ifs_tpfa solver.");
}
// Compute half-transmissibilities, gravity contributions.
int num_cells = grid.numCells();
int ngconn = grid_.c_grid()->cell_facepos[num_cells];
gpress_.clear();
gpress_.resize(ngconn, 0.0);
ncf_.resize(num_cells);
typename Grid::Vector grav;
std::copy(gravity, gravity + Grid::dimension, &grav[0]);
int count = 0;
for (int cell = 0; cell < num_cells; ++cell) {
int num_local_faces = grid.numCellFaces(cell);
ncf_[cell] = num_local_faces;
typename Grid::Vector cc = grid.cellCentroid(cell);
for (int local_ix = 0; local_ix < num_local_faces; ++local_ix) {
int face = grid.cellFace(cell, local_ix);
typename Grid::Vector fc = grid.faceCentroid(face);
gpress_[count++] = grav*(fc - cc);
}
}
assert(count == ngconn);
htrans_.resize(ngconn);
tpfa_htrans_compute(grid_.c_grid(), perm, &htrans_[0]);
state_ = Initialized;
}
enum FlowBCTypes { FBC_UNSET, FBC_PRESSURE, FBC_FLUX };
/// @brief
/// Assemble the sparse system.
/// You must call init() prior to calling assemble().
/// @param sources Source terms, one per cell. Positive numbers
/// are sources, negative are sinks.
/// @param total_mobilities Scalar total mobilities, one per cell.
/// @param omegas Gravity term, one per cell. In a multi-phase
/// flow setting this is equal to
/// \f[ \omega = \sum_{p} \frac{\lambda_p}{\lambda_t} \rho_p \f]
/// where \f$\lambda_p\f$ is a phase mobility, \f$\rho_p\f$ is a
/// phase density and \f$\lambda_t\f$ is the total mobility.
void assemble(const std::vector<double>& sources,
const std::vector<double>& total_mobilities,
const std::vector<double>& omegas,
const std::vector<FlowBCTypes>& bctypes,
const std::vector<double>& bcvalues)
{
if (state_ == Uninitialized) {
throw std::runtime_error("Error in TPFAPressureSolver::assemble(): You must call init() prior to calling assemble().");
}
UnstructuredGrid* g = grid_.c_grid();
// Source terms from user.
double* src = const_cast<double*>(&sources[0]); // Ugly? Yes. Safe? I think so.
// All well related things are zero.
// well_control_t* wctrl = 0;
// double* WI = 0;
// double* wdp = 0;
// Compute effective transmissibilities.
eff_trans_.resize(grid_.numFaces());
tpfa_eff_trans_compute(g, &total_mobilities[0], &htrans_[0], &eff_trans_[0]);
// Update gravity term.
gpress_omegaweighted_.resize(gpress_.size());
mim_ip_density_update(g->number_of_cells, g->cell_facepos, &omegas[0],
&gpress_[0], &gpress_omegaweighted_[0]);
// Zero the linalg structures.
csrmatrix_zero(data_->A);
for (std::size_t i = 0; i < data_->A->m; i++) {
data_->b[i] = 0.0;
}
forces_.src = &src[0];
forces_.bc = 0;
forces_.W = 0;
// Assemble the embedded linear system.
int ok = ifs_tpfa_assemble(g, &forces_, &eff_trans_[0], &gpress_[0], data_);
if (!ok) {
THROW("Assembly of pressure system failed.");
}
state_ = Assembled;
}
/// Encapsulate a sparse linear system in CSR format.
struct LinearSystem
{
int n;
int nnz;
int* ia;
int* ja;
double* sa;
double* b;
double* x;
};
/// @brief
/// Access the linear system assembled.
/// You must call assemble() prior to calling linearSystem().
/// @param[out] s The linear system encapsulation to modify.
/// After this call, s will point to linear system structures
/// that are owned and allocated internally.
void linearSystem(LinearSystem& s)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TPFAPressureSolver::linearSystem(): "
"You must call assemble() prior to calling linearSystem().");
}
s.n = data_->A->m;
s.nnz = data_->A->nnz;
s.ia = data_->A->ia;
s.ja = data_->A->ja;
s.sa = data_->A->sa;
s.b = data_->b;
s.x = data_->x;
}
/// @brief
/// Compute cell pressures and face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// computePressuresAndFluxes().
/// @param[out] cell_pressures Cell pressure values.
/// @param[out] face_areas Face flux values.
void computePressuresAndFluxes(std::vector<double>& cell_pressures,
std::vector<double>& face_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TPFAPressureSolver::computePressuresAndFluxes(): "
"You must call assemble() (and solve the linear system) "
"prior to calling computePressuresAndFluxes().");
}
int num_cells = grid_.c_grid()->number_of_cells;
int num_faces = grid_.c_grid()->number_of_faces;
cell_pressures.clear();
cell_pressures.resize(num_cells, 0.0);
face_fluxes.clear();
face_fluxes.resize(num_faces, 0.0);
ifs_tpfa_solution soln = { 0 };
soln.cell_press = &cell_pressures[0];
soln.face_flux = &face_fluxes [0];
ifs_tpfa_press_flux(grid_.c_grid(), &forces_, &eff_trans_[0],
data_, &soln);
}
/// @brief
/// Compute cell fluxes from face fluxes.
/// You must call assemble() (and solve the linear system accessed
/// by calling linearSystem()) prior to calling
/// faceFluxToCellFlux().
/// @param face_fluxes
/// @param face_areas Face flux values (usually output from computePressuresAndFluxes()).
/// @param[out] cell_fluxes Cell-wise flux values.
/// They are given in cell order, and for each cell there is
/// one value for each adjacent face (in the same order as the
/// cell-face topology of the grid). Positive values represent
/// fluxes out of the cell.
void faceFluxToCellFlux(const std::vector<double>& face_fluxes,
std::vector<double>& cell_fluxes)
{
if (state_ != Assembled) {
throw std::runtime_error("Error in TPFAPressureSolver::faceFluxToCellFlux(): "
"You must call assemble() (and solve the linear system) "
"prior to calling faceFluxToCellFlux().");
}
const UnstructuredGrid& g = *(grid_.c_grid());
int num_cells = g.number_of_cells;
cell_fluxes.resize(g.cell_facepos[num_cells]);
for (int cell = 0; cell < num_cells; ++cell) {
for (int hface = g.cell_facepos[cell]; hface < g.cell_facepos[cell + 1]; ++hface) {
int face = g.cell_faces[hface];
bool pos = (g.face_cells[2*face] == cell);
cell_fluxes[hface] = pos ? face_fluxes[face] : -face_fluxes[face];
}
}
}
/// @brief
/// Access the number of connections (faces) per cell. Deprecated, will be removed.
const std::vector<int>& numCellFaces()
{
return ncf_;
}
private:
// Disabling copy and assigment for now.
TPFAPressureSolver(const TPFAPressureSolver&);
TPFAPressureSolver& operator=(const TPFAPressureSolver&);
enum State { Uninitialized, Initialized, Assembled };
State state_;
// Solver data.
ifs_tpfa_data* data_;
ifs_tpfa_forces forces_;
// Grid.
GridAdapter grid_;
// Number of faces per cell.
std::vector<int> ncf_;
// Transmissibility storage.
std::vector<double> htrans_;
std::vector<double> eff_trans_;
// Gravity contributions.
std::vector<double> gpress_;
std::vector<double> gpress_omegaweighted_;
// Total mobilities.
std::vector<double> totmob_;
// Gravity coefficients (\omega = sum_{i = 1}^{num phases}f_i \rho_i[TODO: check this]).
std::vector<double> omega_;
};
#endif // OPM_TPFAPRESSURESOLVER_HEADER_INCLUDED

2
opm/core/pressure/fsh.h
View File

@@ -52,7 +52,7 @@
*/
#include <opm/core/grid.h>
#include <opm/core/well.h>
#include <opm/core/pressure/legacy_well.h>
#include <opm/core/pressure/flow_bc.h>
#ifdef __cplusplus

2
opm/core/pressure/fsh_common_impl.c
View File

@@ -25,7 +25,7 @@
#include <string.h>
#include <opm/core/grid.h>
#include <opm/core/well.h>
#include <opm/core/pressure/legacy_well.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/pressure/fsh.h>

2
opm/core/pressure/well.c → opm/core/pressure/legacy_well.c
View File

@@ -19,7 +19,7 @@
#include <stdlib.h>
#include <opm/core/well.h>
#include <opm/core/pressure/legacy_well.h>
/* Release memory resources for cell->well mapping. */

6
opm/core/well.h → opm/core/pressure/legacy_well.h
View File

@@ -17,8 +17,8 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELL_HEADER_INCLUDED
#define OPM_WELL_HEADER_INCLUDED
#ifndef OPM_LEGACY_WELL_HEADER_INCLUDED
#define OPM_LEGACY_WELL_HEADER_INCLUDED
/**
* \file
@@ -133,4 +133,4 @@ derive_cell_wells(int nc, well_t *W, int *cwpos, int *cwells);
}
#endif
#endif /* OPM_WELL_HEADER_INCLUDED */
#endif /* OPM_LEGACY_WELL_HEADER_INCLUDED */

2
opm/core/pressure/mimetic/hybsys_global.h
View File

@@ -93,7 +93,7 @@ extern "C" {
#endif
#include <opm/core/grid.h>
#include <opm/core/well.h>
#include <opm/core/pressure/legacy_well.h>
#include <opm/core/linalg/sparse_sys.h>

2
opm/core/pressure/tpfa/cfs_tpfa.c
View File

@@ -4,7 +4,7 @@
#include <stdlib.h>
#include <string.h>
#include <opm/core/well.h>
#include <opm/core/pressure/legacy_well.h>
#include <opm/core/linalg/blas_lapack.h>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/pressure/flow_bc.h>

2
opm/core/pressure/tpfa/cfs_tpfa.h
View File

@@ -21,7 +21,7 @@
#define OPM_CFS_TPFA_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/well.h>
#include <opm/core/pressure/legacy_well.h>
#include <opm/core/pressure/flow_bc.h>
#ifdef __cplusplus

2
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@@ -4,7 +4,7 @@
#include <stdlib.h>
#include <string.h>
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <opm/core/linalg/blas_lapack.h>
#include <opm/core/linalg/sparse_sys.h>

2
opm/core/pressure/tpfa/cfs_tpfa_residual.h
View File

@@ -21,7 +21,7 @@
#define OPM_CFS_TPFA_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/tpfa/compr_source.h>

2
opm/core/pressure/tpfa/ifs_tpfa.c
View File

@@ -6,7 +6,7 @@
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/pressure/tpfa/ifs_tpfa.h>

0
opm/core/fluid/blackoil/BlackoilPhases.hpp → opm/core/props/BlackoilPhases.hpp
View File

2
opm/core/fluid/BlackoilPropertiesBasic.cpp → opm/core/props/BlackoilPropertiesBasic.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>

8
opm/core/fluid/BlackoilPropertiesBasic.hpp → opm/core/props/BlackoilPropertiesBasic.hpp
View File

@@ -21,10 +21,10 @@
#define OPM_BLACKOILPROPERTIESBASIC_HEADER_INCLUDED
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/RockBasic.hpp>
#include <opm/core/fluid/PvtPropertiesBasic.hpp>
#include <opm/core/fluid/SaturationPropsBasic.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/rock/RockBasic.hpp>
#include <opm/core/props/pvt/PvtPropertiesBasic.hpp>
#include <opm/core/props/satfunc/SaturationPropsBasic.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
namespace Opm

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

@@ -17,16 +17,19 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
namespace Opm
{
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid)
const UnstructuredGrid& grid,
bool init_rock)
{
rock_.init(deck, grid);
if (init_rock){
rock_.init(deck, grid);
}
pvt_.init(deck, 200);
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
@@ -41,9 +44,13 @@ namespace Opm
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param)
const parameter::ParameterGroup& param,
bool init_rock)
{
rock_.init(deck, grid);
if(init_rock){
rock_.init(deck, grid);
}
const int pvt_samples = param.getDefault("pvt_tab_size", 200);
pvt_.init(deck, pvt_samples);

13
opm/core/fluid/BlackoilPropertiesFromDeck.hpp → opm/core/props/BlackoilPropertiesFromDeck.hpp
View File

@@ -21,10 +21,10 @@
#define OPM_BLACKOILPROPERTIESFROMDECK_HEADER_INCLUDED
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/RockFromDeck.hpp>
#include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <boost/scoped_ptr.hpp>
@@ -45,7 +45,7 @@ namespace Opm
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid);
const UnstructuredGrid& grid, bool init_rock=true );
/// Initialize from deck, grid and parameters.
/// \param[in] deck Deck input parser
@@ -60,7 +60,8 @@ namespace Opm
/// be done, and the input fluid data used directly for linear interpolation.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param);
const parameter::ParameterGroup& param,
bool init_rock=true);
/// Destructor.
virtual ~BlackoilPropertiesFromDeck();

0
opm/core/fluid/BlackoilPropertiesInterface.hpp → opm/core/props/BlackoilPropertiesInterface.hpp
View File

2
opm/core/fluid/IncompPropertiesBasic.cpp → opm/core/props/IncompPropertiesBasic.cpp
View File

@@ -19,7 +19,7 @@
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/props/IncompPropertiesBasic.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>

8
opm/core/fluid/IncompPropertiesBasic.hpp → opm/core/props/IncompPropertiesBasic.hpp
View File

@@ -20,10 +20,10 @@
#ifndef OPM_INCOMPPROPERTIESBASIC_HEADER_INCLUDED
#define OPM_INCOMPPROPERTIESBASIC_HEADER_INCLUDED
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/RockBasic.hpp>
#include <opm/core/fluid/PvtPropertiesBasic.hpp>
#include <opm/core/fluid/SaturationPropsBasic.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/rock/RockBasic.hpp>
#include <opm/core/props/pvt/PvtPropertiesBasic.hpp>
#include <opm/core/props/satfunc/SaturationPropsBasic.hpp>
namespace Opm
{

2
opm/core/fluid/IncompPropertiesFromDeck.cpp → opm/core/props/IncompPropertiesFromDeck.cpp
View File

@@ -18,7 +18,7 @@
*/
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <opm/core/props/IncompPropertiesFromDeck.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>

8
opm/core/fluid/IncompPropertiesFromDeck.hpp → opm/core/props/IncompPropertiesFromDeck.hpp
View File

@@ -20,10 +20,10 @@
#ifndef OPM_INCOMPPROPERTIESFROMDECK_HEADER_INCLUDED
#define OPM_INCOMPPROPERTIESFROMDECK_HEADER_INCLUDED
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/RockFromDeck.hpp>
#include <opm/core/fluid/PvtPropertiesIncompFromDeck.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
struct UnstructuredGrid;

0
opm/core/fluid/IncompPropertiesInterface.hpp → opm/core/props/IncompPropertiesInterface.hpp
View File

2
opm/core/fluid/blackoil/phaseUsageFromDeck.hpp → opm/core/props/phaseUsageFromDeck.hpp
View File

@@ -22,7 +22,7 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm

14
opm/core/fluid/blackoil/BlackoilPvtProperties.cpp → opm/core/props/pvt/BlackoilPvtProperties.cpp
View File

@@ -19,13 +19,13 @@
#include <opm/core/fluid/blackoil/BlackoilPvtProperties.hpp>
#include <opm/core/fluid/blackoil/SinglePvtDead.hpp>
#include <opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp>
#include <opm/core/fluid/blackoil/SinglePvtLiveOil.hpp>
#include <opm/core/fluid/blackoil/SinglePvtLiveGas.hpp>
#include <opm/core/fluid/blackoil/SinglePvtConstCompr.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/props/pvt/SinglePvtDead.hpp>
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/props/pvt/SinglePvtConstCompr.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>

4
opm/core/fluid/blackoil/BlackoilPvtProperties.hpp → opm/core/props/pvt/BlackoilPvtProperties.hpp
View File

@@ -22,8 +22,8 @@
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <string>
#include <tr1/memory>

2
opm/core/fluid/PvtPropertiesBasic.cpp → opm/core/props/pvt/PvtPropertiesBasic.cpp
View File

@@ -19,7 +19,7 @@
#include <opm/core/fluid/PvtPropertiesBasic.hpp>
#include <opm/core/props/pvt/PvtPropertiesBasic.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>

0
opm/core/fluid/PvtPropertiesBasic.hpp → opm/core/props/pvt/PvtPropertiesBasic.hpp
View File

6
opm/core/fluid/PvtPropertiesIncompFromDeck.cpp → opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
View File

@@ -18,12 +18,12 @@
*/
#include <opm/core/fluid/PvtPropertiesIncompFromDeck.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm

0
opm/core/fluid/PvtPropertiesIncompFromDeck.hpp → opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
View File

2
opm/core/fluid/blackoil/SinglePvtConstCompr.hpp → opm/core/props/pvt/SinglePvtConstCompr.hpp
View File

@@ -21,7 +21,7 @@
#define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <vector>
#include <algorithm>

2
opm/core/fluid/blackoil/SinglePvtDead.cpp → opm/core/props/pvt/SinglePvtDead.cpp
View File

@@ -18,7 +18,7 @@
*/
#include <opm/core/fluid/blackoil/SinglePvtDead.hpp>
#include <opm/core/props/pvt/SinglePvtDead.hpp>
#include <algorithm>
// Extra includes for debug dumping of tables.

2
opm/core/fluid/blackoil/SinglePvtDead.hpp → opm/core/props/pvt/SinglePvtDead.hpp
View File

@@ -21,7 +21,7 @@
#define OPM_SINGLEPVTDEAD_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <vector>

2
opm/core/fluid/blackoil/SinglePvtDeadSpline.cpp → opm/core/props/pvt/SinglePvtDeadSpline.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp>
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <algorithm>

2
opm/core/fluid/blackoil/SinglePvtDeadSpline.hpp → opm/core/props/pvt/SinglePvtDeadSpline.hpp
View File

@@ -21,7 +21,7 @@
#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <vector>

2
opm/core/fluid/blackoil/SinglePvtInterface.cpp → opm/core/props/pvt/SinglePvtInterface.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
namespace Opm

2
opm/core/fluid/blackoil/SinglePvtInterface.hpp → opm/core/props/pvt/SinglePvtInterface.hpp
View File

@@ -21,7 +21,7 @@
#define OPM_SINGLEPVTINTERFACE_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm

2
opm/core/fluid/blackoil/SinglePvtLiveGas.cpp → opm/core/props/pvt/SinglePvtLiveGas.cpp
View File

@@ -28,7 +28,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/blackoil/SinglePvtLiveGas.hpp>
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <algorithm>

2
opm/core/fluid/blackoil/SinglePvtLiveGas.hpp → opm/core/props/pvt/SinglePvtLiveGas.hpp
View File

@@ -20,7 +20,7 @@
#ifndef OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
#define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <vector>
namespace Opm

2
opm/core/fluid/blackoil/SinglePvtLiveOil.cpp → opm/core/props/pvt/SinglePvtLiveOil.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/blackoil/SinglePvtLiveOil.hpp>
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linInt.hpp>
#include <algorithm>

2
opm/core/fluid/blackoil/SinglePvtLiveOil.hpp → opm/core/props/pvt/SinglePvtLiveOil.hpp
View File

@@ -21,7 +21,7 @@
#define OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/SinglePvtInterface.hpp>
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <vector>
namespace Opm

2
opm/core/fluid/RockBasic.cpp → opm/core/props/rock/RockBasic.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/RockBasic.hpp>
#include <opm/core/props/rock/RockBasic.hpp>
namespace Opm
{

0
opm/core/fluid/RockBasic.hpp → opm/core/props/rock/RockBasic.hpp
View File

2
opm/core/fluid/RockCompressibility.cpp → opm/core/props/rock/RockCompressibility.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>

0
opm/core/fluid/RockCompressibility.hpp → opm/core/props/rock/RockCompressibility.hpp
View File

2
opm/core/fluid/RockFromDeck.cpp → opm/core/props/rock/RockFromDeck.cpp
View File

@@ -18,7 +18,7 @@
*/
#include <opm/core/fluid/RockFromDeck.hpp>
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/grid.h>
#include <tr1/array>

0
opm/core/fluid/RockFromDeck.hpp → opm/core/props/rock/RockFromDeck.hpp
View File

8
opm/core/fluid/SatFuncGwseg.cpp → opm/core/props/satfunc/SatFuncGwseg.cpp
View File

@@ -17,11 +17,11 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/SatFuncGwseg.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/props/satfunc/SatFuncGwseg.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>

2
opm/core/fluid/SatFuncGwseg.hpp → opm/core/props/satfunc/SatFuncGwseg.hpp
View File

@@ -22,7 +22,7 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
namespace Opm

8
opm/core/fluid/SatFuncSimple.cpp → opm/core/props/satfunc/SatFuncSimple.cpp
View File

@@ -17,11 +17,11 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/SatFuncSimple.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>

2
opm/core/fluid/SatFuncSimple.hpp → opm/core/props/satfunc/SatFuncSimple.hpp
View File

@@ -22,7 +22,7 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
namespace Opm

8
opm/core/fluid/SatFuncStone2.cpp → opm/core/props/satfunc/SatFuncStone2.cpp
View File

@@ -17,11 +17,11 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/SatFuncStone2.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>

2
opm/core/fluid/SatFuncStone2.hpp → opm/core/props/satfunc/SatFuncStone2.hpp
View File

@@ -22,7 +22,7 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
namespace Opm

2
opm/core/fluid/SaturationPropsBasic.cpp → opm/core/props/satfunc/SaturationPropsBasic.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/SaturationPropsBasic.hpp>
#include <opm/core/props/satfunc/SaturationPropsBasic.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <iostream>

0
opm/core/fluid/SaturationPropsBasic.hpp → opm/core/props/satfunc/SaturationPropsBasic.hpp
View File

2
opm/core/fluid/SaturationPropsFromDeck.cpp → opm/core/props/satfunc/SaturationPropsFromDeck.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/fluid/SaturationPropsFromDeck.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/ErrorMacros.hpp>

12
opm/core/fluid/SaturationPropsFromDeck.hpp → opm/core/props/satfunc/SaturationPropsFromDeck.hpp
View File

@@ -20,13 +20,13 @@
#ifndef OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
#define OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED
#include <opm/core/fluid/SaturationPropsInterface.hpp>
#include <opm/core/props/satfunc/SaturationPropsInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/fluid/SatFuncStone2.hpp>
#include <opm/core/fluid/SatFuncSimple.hpp>
#include <opm/core/fluid/SatFuncGwseg.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
#include <opm/core/props/satfunc/SatFuncGwseg.hpp>
#include <vector>
struct UnstructuredGrid;
@@ -136,7 +136,7 @@ namespace Opm
} // namespace Opm
#include <opm/core/fluid/SaturationPropsFromDeck_impl.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp>
#endif // OPM_SATURATIONPROPSFROMDECK_HEADER_INCLUDED

2
opm/core/fluid/SaturationPropsFromDeck_impl.hpp → opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
View File

@@ -23,7 +23,7 @@
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/grid.h>
namespace Opm

2
opm/core/fluid/SaturationPropsInterface.hpp → opm/core/props/satfunc/SaturationPropsInterface.hpp
View File

@@ -20,7 +20,7 @@
#ifndef OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
#define OPM_SATURATIONPROPSINTERFACE_HEADER_INCLUDED
#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm

2
opm/core/simulator/BlackoilState.hpp
View File

@@ -21,7 +21,7 @@
#define OPM_BLACKOILSTATE_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <vector>
namespace Opm

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

@@ -28,7 +28,7 @@
#include <opm/core/pressure/CompressibleTpfa.hpp>
#include <opm/core/grid.h>
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
@@ -40,13 +40,13 @@
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/grid/ColumnExtract.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportModelCompressibleTwophase.hpp>
#include <opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
@@ -101,7 +101,7 @@ namespace Opm
const double* gravity_;
// Solvers
CompressibleTpfa psolver_;
TransportModelCompressibleTwophase tsolver_;
TransportSolverCompressibleTwophaseReorder tsolver_;
// Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_;
// Misc. data

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

@@ -28,7 +28,7 @@
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/grid.h>
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <opm/core/pressure/flow_bc.h>
#include <opm/core/simulator/SimulatorReport.hpp>
@@ -39,14 +39,13 @@
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/RockCompressibility.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/core/utility/ColumnExtract.hpp>
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
#include <opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp>
#include <boost/filesystem.hpp>
#include <boost/scoped_ptr.hpp>
#include <boost/lexical_cast.hpp>
@@ -87,6 +86,7 @@ namespace Opm
int max_well_control_iterations_;
// Parameters for transport solver.
int num_transport_substeps_;
bool use_reorder_;
bool use_segregation_split_;
// Observed objects.
const UnstructuredGrid& grid_;
@@ -98,9 +98,7 @@ namespace Opm
const FlowBoundaryConditions* bcs_;
// Solvers
IncompTpfa psolver_;
TransportModelTwophase tsolver_;
// Needed by column-based gravity segregation solver.
std::vector< std::vector<int> > columns_;
boost::scoped_ptr<TransportSolverTwophaseInterface> tsolver_;
// Misc. data
std::vector<int> allcells_;
};
@@ -318,7 +316,9 @@ namespace Opm
const FlowBoundaryConditions* bcs,
LinearSolverInterface& linsolver,
const double* gravity)
: grid_(grid),
: use_reorder_(param.getDefault("use_reorder", true)),
use_segregation_split_(param.getDefault("use_segregation_split", false)),
grid_(grid),
props_(props),
rock_comp_props_(rock_comp_props),
wells_manager_(wells_manager),
@@ -329,11 +329,33 @@ namespace Opm
param.getDefault("nl_pressure_residual_tolerance", 0.0),
param.getDefault("nl_pressure_change_tolerance", 1.0),
param.getDefault("nl_pressure_maxiter", 10),
gravity, wells_manager.c_wells(), src, bcs),
tsolver_(grid, props,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30))
gravity, wells_manager.c_wells(), src, bcs)
{
// Initialize transport solver.
if (use_reorder_) {
tsolver_.reset(new Opm::TransportSolverTwophaseReorder(grid,
props,
use_segregation_split_ ? gravity : NULL,
param.getDefault("nl_tolerance", 1e-9),
param.getDefault("nl_maxiter", 30)));
} else {
if (rock_comp_props && rock_comp_props->isActive()) {
THROW("The implicit pressure solver cannot handle rock compressibility.");
}
if (use_segregation_split_) {
THROW("The implicit pressure solver is not set up to use segregation splitting.");
}
std::vector<double> porevol;
computePorevolume(grid, props.porosity(), porevol);
tsolver_.reset(new Opm::TransportSolverTwophaseImplicit(grid,
props,
porevol,
gravity,
psolver_.getHalfTrans(),
param));
}
// For output.
output_ = param.getDefault("output", true);
if (output_) {
@@ -356,11 +378,6 @@ namespace Opm
// Transport related init.
num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
use_segregation_split_ = param.getDefault("use_segregation_split", false);
if (gravity != 0 && use_segregation_split_){
tsolver_.initGravity(gravity);
extractColumn(grid_, columns_);
}
// Misc init.
const int num_cells = grid.number_of_cells;
@@ -427,9 +444,12 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputVectorMatlab(std::string("reorder_it"),
tsolver_.getReorderIterations(),
timer.currentStepNum(), output_dir_);
if (use_reorder_) {
// This use of dynamic_cast is not ideal, but should be safe.
outputVectorMatlab(std::string("reorder_it"),
dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
timer.currentStepNum(), output_dir_);
}
}
SimulatorReport sreport;
@@ -523,8 +543,8 @@ namespace Opm
double injected[2] = { 0.0 };
double produced[2] = { 0.0 };
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0],
stepsize, state.saturation());
tsolver_->solve(&initial_porevol[0], &transport_src[0], stepsize, state);
double substep_injected[2] = { 0.0 };
double substep_produced[2] = { 0.0 };
Opm::computeInjectedProduced(props_, state.saturation(), transport_src, stepsize,
@@ -533,8 +553,11 @@ namespace Opm
injected[1] += substep_injected[1];
produced[0] += substep_produced[0];
produced[1] += substep_produced[1];
if (use_segregation_split_) {
tsolver_.solveGravity(columns_, &initial_porevol[0], stepsize, state.saturation());
if (use_reorder_ && use_segregation_split_) {
// Again, unfortunate but safe use of dynamic_cast.
// Possible solution: refactor gravity solver to its own class.
dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
.solveGravity(&initial_porevol[0], stepsize, state);
}
watercut.push(timer.currentTime() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
@@ -571,9 +594,12 @@ namespace Opm
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
}
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
outputVectorMatlab(std::string("reorder_it"),
tsolver_.getReorderIterations(),
timer.currentStepNum(), output_dir_);
if (use_reorder_) {
// This use of dynamic_cast is not ideal, but should be safe.
outputVectorMatlab(std::string("reorder_it"),
dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
timer.currentStepNum(), output_dir_);
}
outputWaterCut(watercut, output_dir_);
if (wells_) {
outputWellReport(wellreport, output_dir_);

2
opm/core/simulator/TwophaseState.hpp
View File

@@ -21,7 +21,7 @@
#define OPM_TWOPHASESTATE_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <vector>
namespace Opm

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

@@ -20,7 +20,7 @@
#ifndef OPM_WELLSTATE_HEADER_INCLUDED
#define OPM_WELLSTATE_HEADER_INCLUDED
#include <opm/core/newwells.h>
#include <opm/core/wells.h>
#include <vector>
namespace Opm

2
opm/core/utility/initState.hpp → opm/core/simulator/initState.hpp
View File

@@ -114,6 +114,6 @@ namespace Opm
} // namespace Opm
#include <opm/core/utility/initState_impl.hpp>
#include <opm/core/simulator/initState_impl.hpp>
#endif // OPM_INITSTATE_HEADER_INCLUDED

6
opm/core/utility/initState_impl.hpp → opm/core/simulator/initState_impl.hpp
View File

@@ -27,9 +27,9 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/MonotCubicInterpolator.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/fluid/IncompPropertiesInterface.hpp>
#include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
#include <opm/core/props/IncompPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <cmath>
namespace Opm

2
opm/core/transport/reorder/DGBasis.cpp → opm/core/tof/DGBasis.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/DGBasis.hpp>
#include <opm/core/tof/DGBasis.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <numeric>

0
opm/core/transport/reorder/DGBasis.hpp → opm/core/tof/DGBasis.hpp
View File

32
opm/core/transport/reorder/TransportModelTracerTofDiscGal.cpp → opm/core/tof/TofDiscGalReorder.cpp
View File

@@ -19,8 +19,8 @@
#include <opm/core/grid/CellQuadrature.hpp>
#include <opm/core/grid/FaceQuadrature.hpp>
#include <opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp>
#include <opm/core/transport/reorder/DGBasis.hpp>
#include <opm/core/tof/TofDiscGalReorder.hpp>
#include <opm/core/tof/DGBasis.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/VelocityInterpolation.hpp>
@@ -34,10 +34,6 @@ namespace Opm
{
// --------------- Methods of TransportModelTracerTofDiscGal ---------------
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] param Parameters for the solver.
@@ -55,8 +51,8 @@ namespace Opm
/// computing (unlimited) solution.
/// AsSimultaneousPostProcess Apply to each cell independently, using un-
/// limited solution in neighbouring cells.
TransportModelTracerTofDiscGal::TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const parameter::ParameterGroup& param)
TofDiscGalReorder::TofDiscGalReorder(const UnstructuredGrid& grid,
const parameter::ParameterGroup& param)
: grid_(grid),
use_cvi_(false),
use_limiter_(false),
@@ -127,7 +123,7 @@ namespace Opm
/// cell comes before the K coefficients corresponding
/// to the second cell etc.
/// K depends on degree and grid dimension.
void TransportModelTracerTofDiscGal::solveTof(const double* darcyflux,
void TofDiscGalReorder::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof_coeff)
@@ -173,7 +169,7 @@ namespace Opm
void TransportModelTracerTofDiscGal::solveSingleCell(const int cell)
void TofDiscGalReorder::solveSingleCell(const int cell)
{
// Residual:
// For each cell K, basis function b_j (spanning V_h),
@@ -396,7 +392,7 @@ namespace Opm
void TransportModelTracerTofDiscGal::solveMultiCell(const int num_cells, const int* cells)
void TofDiscGalReorder::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
@@ -407,7 +403,7 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyLimiter(const int cell, double* tof)
void TofDiscGalReorder::applyLimiter(const int cell, double* tof)
{
switch (limiter_method_) {
case MinUpwindFace:
@@ -424,7 +420,7 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyMinUpwindLimiter(const int cell, const bool face_min, double* tof)
void TofDiscGalReorder::applyMinUpwindLimiter(const int cell, const bool face_min, double* tof)
{
if (basis_func_->degree() != 1) {
THROW("This limiter only makes sense for our DG1 implementation.");
@@ -513,12 +509,12 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyLimiterAsPostProcess()
void TofDiscGalReorder::applyLimiterAsPostProcess()
{
// Apply the limiter sequentially to all cells.
// This means that a cell's limiting behaviour may be affected by
// any limiting applied to its upstream cells.
const std::vector<int>& seq = TransportModelInterface::sequence();
const std::vector<int>& seq = ReorderSolverInterface::sequence();
const int nc = seq.size();
ASSERT(nc == grid_.number_of_cells);
for (int i = 0; i < nc; ++i) {
@@ -530,7 +526,7 @@ namespace Opm
void TransportModelTracerTofDiscGal::applyLimiterAsSimultaneousPostProcess()
void TofDiscGalReorder::applyLimiterAsSimultaneousPostProcess()
{
// Apply the limiter simultaneously to all cells.
// This means that each cell is limited independently from all other cells,
@@ -547,7 +543,7 @@ namespace Opm
double TransportModelTracerTofDiscGal::totalFlux(const int cell) const
double TofDiscGalReorder::totalFlux(const int cell) const
{
// Find total upstream/downstream fluxes.
double upstream_flux = 0.0;
@@ -575,7 +571,7 @@ namespace Opm
double TransportModelTracerTofDiscGal::minCornerVal(const int cell, const int face) const
double TofDiscGalReorder::minCornerVal(const int cell, const int face) const
{
// Evaluate the solution in all corners.
const int dim = grid_.dimensions;

16
opm/core/transport/reorder/TransportModelTracerTofDiscGal.hpp → opm/core/tof/TofDiscGalReorder.hpp
View File

@@ -17,10 +17,10 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOFDISCGAL_HEADER_INCLUDED
#ifndef OPM_TOFDISCGALREORDER_HEADER_INCLUDED
#define OPM_TOFDISCGALREORDER_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/transport/reorder/ReorderSolverInterface.hpp>
#include <boost/shared_ptr.hpp>
#include <vector>
#include <map>
@@ -45,7 +45,7 @@ namespace Opm
/// \tau is specified to be zero on all inflow boundaries.
/// The user may specify the polynomial degree of the basis function space
/// used, but only degrees 0 and 1 are supported so far.
class TransportModelTracerTofDiscGal : public TransportModelInterface
class TofDiscGalReorder : public ReorderSolverInterface
{
public:
/// Construct solver.
@@ -68,8 +68,8 @@ namespace Opm
/// computing (unlimited) solution.
/// AsSimultaneousPostProcess Apply to each cell independently, using un-
/// limited solution in neighbouring cells.
TransportModelTracerTofDiscGal(const UnstructuredGrid& grid,
const parameter::ParameterGroup& param);
TofDiscGalReorder(const UnstructuredGrid& grid,
const parameter::ParameterGroup& param);
/// Solve for time-of-flight.
@@ -95,8 +95,8 @@ namespace Opm
private:
// Disable copying and assignment.
TransportModelTracerTofDiscGal(const TransportModelTracerTofDiscGal&);
TransportModelTracerTofDiscGal& operator=(const TransportModelTracerTofDiscGal&);
TofDiscGalReorder(const TofDiscGalReorder&);
TofDiscGalReorder& operator=(const TofDiscGalReorder&);
// Data members
const UnstructuredGrid& grid_;

32
opm/core/transport/reorder/TransportModelTracerTof.cpp → opm/core/tof/TofReorder.cpp
View File

@@ -17,7 +17,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/TransportModelTracerTof.hpp>
#include <opm/core/tof/TofReorder.hpp>
#include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
@@ -30,8 +30,8 @@ namespace Opm
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
TransportModelTracerTof::TransportModelTracerTof(const UnstructuredGrid& grid,
const bool use_multidim_upwind)
TofReorder::TofReorder(const UnstructuredGrid& grid,
const bool use_multidim_upwind)
: grid_(grid),
darcyflux_(0),
porevolume_(0),
@@ -53,7 +53,7 @@ namespace Opm
/// (+) inflow flux,
/// (-) outflow flux.
/// \param[out] tof Array of time-of-flight values.
void TransportModelTracerTof::solveTof(const double* darcyflux,
void TofReorder::solveTof(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof)
@@ -93,11 +93,11 @@ namespace Opm
/// \param[out] tof Array of time-of-flight values (1 per cell).
/// \param[out] tracer Array of tracer values (N per cell, where N is
/// the number of cells c for which source[c] > 0.0).
void TransportModelTracerTof::solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof,
std::vector<double>& tracer)
void TofReorder::solveTofTracer(const double* darcyflux,
const double* porevolume,
const double* source,
std::vector<double>& tof,
std::vector<double>& tracer)
{
darcyflux_ = darcyflux;
porevolume_ = porevolume;
@@ -137,7 +137,7 @@ namespace Opm
void TransportModelTracerTof::solveSingleCell(const int cell)
void TofReorder::solveSingleCell(const int cell)
{
if (use_multidim_upwind_) {
solveSingleCellMultidimUpwind(cell);
@@ -194,7 +194,7 @@ namespace Opm
void TransportModelTracerTof::solveSingleCellMultidimUpwind(const int cell)
void TofReorder::solveSingleCellMultidimUpwind(const int cell)
{
// Compute flux terms.
// Sources have zero tof, and therefore do not contribute
@@ -243,7 +243,7 @@ namespace Opm
void TransportModelTracerTof::solveMultiCell(const int num_cells, const int* cells)
void TofReorder::solveMultiCell(const int num_cells, const int* cells)
{
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
@@ -260,10 +260,10 @@ namespace Opm
// tof(face) = face_term + cell_term_factor*tof(upwind_cell).
// It is not computed here, since these factors are needed to
// compute the tof(upwind_cell) itself.
void TransportModelTracerTof::multidimUpwindTerms(const int face,
const int upwind_cell,
double& face_term,
double& cell_term_factor) const
void TofReorder::multidimUpwindTerms(const int face,
const int upwind_cell,
double& face_term,
double& cell_term_factor) const
{
// Implements multidim upwind according to
// "Multidimensional upstream weighting for multiphase transport on general grids"

12
opm/core/transport/reorder/TransportModelTracerTof.hpp → opm/core/tof/TofReorder.hpp
View File

@@ -17,10 +17,10 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#define OPM_TRANSPORTMODELTRACERTOF_HEADER_INCLUDED
#ifndef OPM_TOFREORDER_HEADER_INCLUDED
#define OPM_TOFREORDER_HEADER_INCLUDED
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/transport/reorder/ReorderSolverInterface.hpp>
#include <vector>
#include <map>
#include <ostream>
@@ -38,14 +38,14 @@ namespace Opm
/// where v is the fluid velocity, \tau is time-of-flight and
/// \phi is the porosity. This is a boundary value problem, where
/// \tau is specified to be zero on all inflow boundaries.
class TransportModelTracerTof : public TransportModelInterface
class TofReorder : public ReorderSolverInterface
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding.
TransportModelTracerTof(const UnstructuredGrid& grid,
const bool use_multidim_upwind = false);
TofReorder(const UnstructuredGrid& grid,
const bool use_multidim_upwind = false);
/// Solve for time-of-flight.
/// \param[in] darcyflux Array of signed face fluxes.

72
opm/core/transport/SimpleFluid2pWrapper.hpp
View File

@@ -1,72 +0,0 @@
/*===========================================================================
//
// File: SimpleFluid2pWrapper.hpp
//
// Created: 2011-09-30 11:38:28+0200
//
// Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
// Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_SimpleFluid2pWrapper_HPP_HEADER
#define OPM_SimpleFluid2pWrapper_HPP_HEADER
#include <cmath>
namespace Opm {
template <class ReservoirProperties>
class SimpleFluid2pWrapper {
public:
SimpleFluid2pWrapper (const ReservoirProperties& resprop)
{
resprop_ = resprop;
}
double density(int p) const { return 0; }
template <class Sat ,
class Mob ,
class DMob>
void
mobility(int c, const Sat& s, Mob& mob, DMob& dmob) const {
const double s1 = s[0];
const double s2 = 1 - s1;
mob[0] = resprop_.mobilityFirstPhase(c, s[0]) ;
mob[1] = resprop_.mobilitySecondPhase(c, s[0]) ;
dmob[0] = resprop_.dmobilityFirstPhase(c, s[0]) ;
dmob[3] = -resprop_.dmobilitySecondPhase(c, s[0]) ;
dmob[1] = dmob[2] = 0.0;
}
private:
ReservoirProperties& resprop_;
};
}
#endif /* OPM_SimpleFluid2pWrapper_HPP_HEADER */

34
opm/core/transport/TransportSolverTwophaseInterface.cpp Normal file
View File

@@ -0,0 +1,34 @@
/*===========================================================================
//
// File: TwoPhaseTransportSolver.cpp
//
// Author: hnil <hnil@sintef.no>
//
// Created: 9 Nov 2012
//==========================================================================*/
/*
Copyright 2011 SINTEF ICT, Applied Mathematics.
Copyright 2011 Statoil ASA.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "TransportSolverTwophaseInterface.hpp"
Opm::TransportSolverTwophaseInterface::~TransportSolverTwophaseInterface()
{
}

50
opm/core/transport/TransportSolverTwophaseInterface.hpp Normal file
View File

@@ -0,0 +1,50 @@
/*
Copyright 2012, 2013 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media Project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_TRANSPORTSOLVERTWOPHASEINTERFACE_HEADER_INCLUDED
#define OPM_TRANSPORTSOLVERTWOPHASEINTERFACE_HEADER_INCLUDED
#include <opm/core/simulator/TwophaseState.hpp>
#include <opm/core/simulator/WellState.hpp>
namespace Opm
{
/// Base class for two-phase incompressible transport solvers.
class TransportSolverTwophaseInterface
{
public:
/// Virtual destructor to enable inheritance.
virtual ~TransportSolverTwophaseInterface();
/// Solve for saturation at next timestep.
/// \param[in] porevolume Array of pore volumes.
/// \param[in] source Transport source term. For interpretation see Opm::computeTransportSource().
/// \param[in] dt Time step.
/// \param[in, out] state Reservoir state. Calling solve() will read state.faceflux() and
/// read and write state.saturation().
virtual void solve(const double* porevolume,
const double* source,
const double dt,
TwophaseState& state) = 0;
};
}
#endif // OPM_TRANSPORTSOLVERTWOPHASEINTERFACE_HEADER_INCLUDED

0
opm/core/transport/CSRMatrixBlockAssembler.hpp → opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp
View File

0
opm/core/transport/CSRMatrixUmfpackSolver.hpp → opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp
View File

0
opm/core/transport/GravityColumnSolver.hpp → opm/core/transport/implicit/GravityColumnSolver.hpp
View File

Some files were not shown because too many files have changed in this diff Show More