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

replace the blackoil PVT classes by the ones of opm-material

Browse Source 
the opm-material classes are the ones which are now used by
opm-autodiff and this patch makes it much easier to keep the opm-core
and opm-autodiff results consistent. Also, the opm-material classes
seem to be a bit faster than the opm-core ones (see
https://github.com/OPM/opm-autodiff/pull/576)

I ran the usual array of tests with `flow`: SPE1, SPE3, SPE9 and Norne
all produce the same results at the identical runtime (modulo noise)
and also "Model 2" seems to work.
This commit is contained in:
Andreas Lauser
2016-02-15 16:03:30 +01:00
parent 0c05fc4dc0
commit 23088f987f
23 changed files with 882 additions and 4365 deletions
Download Patch File Download Diff File Expand all files Collapse all files

110
CMakeLists_files.cmake
View File

@@ -29,30 +29,35 @@
# originally generated with the command:
# find opm -name '*.c*' -printf '\t%p\n' | sort
list (APPEND MAIN_SOURCE_FILES
opm/core/flowdiagnostics/AnisotropicEikonal.cpp
opm/core/flowdiagnostics/DGBasis.cpp
opm/core/flowdiagnostics/FlowDiagnostics.cpp
opm/core/flowdiagnostics/TofDiscGalReorder.cpp
opm/core/flowdiagnostics/TofReorder.cpp
opm/core/grid/GridHelpers.cpp
opm/core/grid/GridManager.cpp
opm/core/grid/GridUtilities.cpp
opm/core/grid/grid.c
opm/core/grid/grid_equal.cpp
opm/core/grid/cart_grid.c
opm/core/grid/cornerpoint_grid.c
opm/core/grid/cpgpreprocess/facetopology.c
opm/core/grid/cpgpreprocess/geometry.c
opm/core/grid/cpgpreprocess/preprocess.c
opm/core/grid/cpgpreprocess/uniquepoints.c
opm/core/grid/grid.c
opm/core/grid/grid_equal.cpp
opm/core/io/OutputWriter.cpp
opm/core/io/eclipse/EclipseGridInspector.cpp
opm/core/io/eclipse/EclipseWriter.cpp
opm/core/io/eclipse/EclipseReader.cpp
opm/core/io/eclipse/EclipseWriteRFTHandler.cpp
opm/core/io/eclipse/EclipseWriter.cpp
opm/core/io/eclipse/writeECLData.cpp
opm/core/io/OutputWriter.cpp
opm/core/io/vag/vag.cpp
opm/core/io/vtk/writeVtkData.cpp
opm/core/linalg/LinearSolverFactory.cpp
opm/core/linalg/LinearSolverInterface.cpp
opm/core/linalg/LinearSolverIstl.cpp
opm/core/linalg/LinearSolverUmfpack.cpp
opm/core/linalg/LinearSolverPetsc.cpp
opm/core/linalg/LinearSolverUmfpack.cpp
opm/core/linalg/call_umfpack.c
opm/core/linalg/sparse_sys.c
opm/core/pressure/CompressibleTpfa.cpp
@@ -64,6 +69,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/pressure/fsh.c
opm/core/pressure/fsh_common_impl.c
opm/core/pressure/ifsh.c
opm/core/pressure/legacy_well.c
opm/core/pressure/mimetic/hybsys.c
opm/core/pressure/mimetic/hybsys_global.c
opm/core/pressure/mimetic/mimetic.c
@@ -73,6 +79,7 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/pressure/msmfem/hash_set.c
opm/core/pressure/msmfem/ifsh_ms.c
opm/core/pressure/msmfem/partition.c
opm/core/pressure/tpfa/TransTpfa.cpp
opm/core/pressure/tpfa/cfs_tpfa.c
opm/core/pressure/tpfa/cfs_tpfa_residual.c
opm/core/pressure/tpfa/compr_bc.c
@@ -80,61 +87,49 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/pressure/tpfa/compr_quant_general.c
opm/core/pressure/tpfa/compr_source.c
opm/core/pressure/tpfa/ifs_tpfa.c
opm/core/pressure/tpfa/TransTpfa.cpp
opm/core/pressure/tpfa/trans_tpfa.c
opm/core/pressure/legacy_well.c
opm/core/props/BlackoilPropertiesBasic.cpp
opm/core/props/BlackoilPropertiesFromDeck.cpp
opm/core/props/IncompPropertiesBasic.cpp
opm/core/props/IncompPropertiesFromDeck.cpp
opm/core/props/IncompPropertiesSinglePhase.cpp
opm/core/props/pvt/BlackoilPvtProperties.cpp
opm/core/props/pvt/PvtPropertiesBasic.cpp
opm/core/props/pvt/PvtPropertiesIncompFromDeck.cpp
opm/core/props/pvt/PvtDead.cpp
opm/core/props/pvt/PvtDeadSpline.cpp
opm/core/props/pvt/PvtInterface.cpp
opm/core/props/pvt/PvtLiveGas.cpp
opm/core/props/pvt/PvtLiveOil.cpp
opm/core/props/rock/RockBasic.cpp
opm/core/props/rock/RockCompressibility.cpp
opm/core/props/rock/RockFromDeck.cpp
opm/core/props/satfunc/RelpermDiagnostics.cpp
opm/core/props/satfunc/SaturationPropsBasic.cpp
opm/core/props/satfunc/SaturationPropsFromDeck.cpp
opm/core/props/satfunc/RelpermDiagnostics.cpp
opm/core/simulator/AdaptiveSimulatorTimer.cpp
opm/core/simulator/BlackoilState.cpp
opm/core/simulator/TimeStepControl.cpp
opm/core/simulator/SimulatorCompressibleTwophase.cpp
opm/core/simulator/SimulatorIncompTwophase.cpp
opm/core/simulator/SimulatorOutput.cpp
opm/core/simulator/SimulatorReport.cpp
opm/core/simulator/SimulatorState.cpp
opm/core/simulator/SimulatorTimer.cpp
opm/core/flowdiagnostics/AnisotropicEikonal.cpp
opm/core/flowdiagnostics/DGBasis.cpp
opm/core/flowdiagnostics/FlowDiagnostics.cpp
opm/core/flowdiagnostics/TofReorder.cpp
opm/core/flowdiagnostics/TofDiscGalReorder.cpp
opm/core/simulator/TimeStepControl.cpp
opm/core/transport/TransportSolverTwophaseInterface.cpp
opm/core/transport/implicit/TransportSolverTwophaseImplicit.cpp
opm/core/transport/implicit/transport_source.c
opm/core/transport/minimal/spu_explicit.c
opm/core/transport/minimal/spu_implicit.c
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
opm/core/transport/reorder/ReorderSolverInterface.cpp
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.cpp
opm/core/transport/reorder/TransportSolverTwophaseReorder.cpp
opm/core/transport/reorder/reordersequence.cpp
opm/core/transport/reorder/tarjan.c
opm/core/utility/compressedToCartesian.cpp
opm/core/utility/Event.cpp
opm/core/utility/MonotCubicInterpolator.cpp
opm/core/utility/NullStream.cpp
opm/core/utility/StopWatch.cpp
opm/core/utility/VelocityInterpolation.cpp
opm/core/utility/WachspressCoord.cpp
opm/core/utility/compressedToCartesian.cpp
opm/core/utility/extractPvtTableIndex.cpp
opm/core/utility/miscUtilities.cpp
opm/core/utility/miscUtilitiesBlackoil.cpp
opm/core/utility/NullStream.cpp
opm/core/utility/parameters/Parameter.cpp
opm/core/utility/parameters/ParameterGroup.cpp
opm/core/utility/parameters/ParameterTools.cpp
@@ -149,8 +144,8 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/wells/WellCollection.cpp
opm/core/wells/WellsGroup.cpp
opm/core/wells/WellsManager.cpp
opm/core/wells/wells.c
opm/core/wells/well_controls.c
opm/core/wells/wells.c
)
# originally generated with the command:
@@ -182,7 +177,6 @@ list (APPEND TEST_SOURCE_FILES
tests/test_linearsolver.cpp
tests/test_parallel_linearsolver.cpp
tests/test_param.cpp
tests/test_blackoilfluid.cpp
tests/test_satfunc.cpp
tests/test_shadow.cpp
tests/test_equil.cpp
@@ -280,6 +274,11 @@ list (APPEND PROGRAM_SOURCE_FILES
# find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort
list (APPEND PUBLIC_HEADER_FILES
opm/core/doxygen_main.hpp
opm/core/flowdiagnostics/AnisotropicEikonal.hpp
opm/core/flowdiagnostics/DGBasis.hpp
opm/core/flowdiagnostics/FlowDiagnostics.hpp
opm/core/flowdiagnostics/TofDiscGalReorder.hpp
opm/core/flowdiagnostics/TofReorder.hpp
opm/core/grid.h
opm/core/grid/CellQuadrature.hpp
opm/core/grid/ColumnExtract.hpp
@@ -295,28 +294,26 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/grid/cpgpreprocess/geometry.h
opm/core/grid/cpgpreprocess/preprocess.h
opm/core/grid/cpgpreprocess/uniquepoints.h
opm/core/io/eclipse/CornerpointChopper.hpp
opm/core/io/eclipse/EclipseIOUtil.hpp
opm/core/io/eclipse/EclipseGridInspector.hpp
opm/core/io/eclipse/EclipseUnits.hpp
opm/core/io/eclipse/EclipseWriter.hpp
opm/core/io/eclipse/EclipseReader.hpp
opm/core/io/eclipse/EclipseWriteRFTHandler.hpp
opm/core/io/eclipse/writeECLData.hpp
opm/core/io/OutputWriter.hpp
opm/core/io/eclipse/CornerpointChopper.hpp
opm/core/io/eclipse/EclipseGridInspector.hpp
opm/core/io/eclipse/EclipseIOUtil.hpp
opm/core/io/eclipse/EclipseReader.hpp
opm/core/io/eclipse/EclipseUnits.hpp
opm/core/io/eclipse/EclipseWriteRFTHandler.hpp
opm/core/io/eclipse/EclipseWriter.hpp
opm/core/io/eclipse/writeECLData.hpp
opm/core/io/vag/vag.hpp
opm/core/io/vtk/writeVtkData.hpp
opm/core/linalg/LinearSolverFactory.hpp
opm/core/linalg/LinearSolverInterface.hpp
opm/core/linalg/LinearSolverIstl.hpp
opm/core/linalg/LinearSolverUmfpack.hpp
opm/core/linalg/LinearSolverPetsc.hpp
opm/core/linalg/LinearSolverUmfpack.hpp
opm/core/linalg/ParallelIstlInformation.hpp
opm/core/linalg/blas_lapack.h
opm/core/linalg/call_umfpack.h
opm/core/linalg/sparse_sys.h
opm/core/wells.h
opm/core/well_controls.h
opm/core/pressure/CompressibleTpfa.hpp
opm/core/pressure/FlowBCManager.hpp
opm/core/pressure/IncompTpfa.hpp
@@ -334,6 +331,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/pressure/msmfem/hash_set.h
opm/core/pressure/msmfem/ifsh_ms.h
opm/core/pressure/msmfem/partition.h
opm/core/pressure/tpfa/TransTpfa.hpp
opm/core/pressure/tpfa/TransTpfa_impl.hpp
opm/core/pressure/tpfa/cfs_tpfa.h
opm/core/pressure/tpfa/cfs_tpfa_residual.h
opm/core/pressure/tpfa/compr_bc.h
@@ -341,8 +340,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/pressure/tpfa/compr_quant_general.h
opm/core/pressure/tpfa/compr_source.h
opm/core/pressure/tpfa/ifs_tpfa.h
opm/core/pressure/tpfa/TransTpfa.hpp
opm/core/pressure/tpfa/TransTpfa_impl.hpp
opm/core/pressure/tpfa/trans_tpfa.h
opm/core/props/BlackoilPhases.hpp
opm/core/props/BlackoilPropertiesBasic.hpp
@@ -355,25 +352,18 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/IncompPropertiesShadow_impl.hpp
opm/core/props/IncompPropertiesSinglePhase.hpp
opm/core/props/phaseUsageFromDeck.hpp
opm/core/props/pvt/BlackoilPvtProperties.hpp
opm/core/props/pvt/PvtPropertiesBasic.hpp
opm/core/props/pvt/PvtPropertiesIncompFromDeck.hpp
opm/core/props/pvt/PvtConstCompr.hpp
opm/core/props/pvt/PvtDead.hpp
opm/core/props/pvt/PvtDeadSpline.hpp
opm/core/props/pvt/PvtInterface.hpp
opm/core/props/pvt/PvtLiveGas.hpp
opm/core/props/pvt/PvtLiveOil.hpp
opm/core/props/pvt/ThermalWaterPvtWrapper.hpp
opm/core/props/pvt/ThermalOilPvtWrapper.hpp
opm/core/props/pvt/ThermalGasPvtWrapper.hpp
opm/core/props/pvt/ThermalOilPvtWrapper.hpp
opm/core/props/pvt/ThermalWaterPvtWrapper.hpp
opm/core/props/rock/RockBasic.hpp
opm/core/props/rock/RockCompressibility.hpp
opm/core/props/rock/RockFromDeck.hpp
opm/core/props/satfunc/RelpermDiagnostics.hpp
opm/core/props/satfunc/SaturationPropsBasic.hpp
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
opm/core/props/satfunc/SaturationPropsInterface.hpp
opm/core/props/satfunc/RelpermDiagnostics.hpp
opm/core/props/satfunc/RelpermDiagnostics_impl.hpp
opm/core/simulator/AdaptiveSimulatorTimer.hpp
opm/core/simulator/AdaptiveTimeStepping.hpp
@@ -383,27 +373,22 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/EquilibrationHelpers.hpp
opm/core/simulator/ExplicitArraysFluidState.hpp
opm/core/simulator/ExplicitArraysSatDerivativesFluidState.hpp
opm/core/simulator/TimeStepControl.hpp
opm/core/simulator/SimulatorCompressibleTwophase.hpp
opm/core/simulator/SimulatorIncompTwophase.hpp
opm/core/simulator/SimulatorOutput.hpp
opm/core/simulator/SimulatorReport.hpp
opm/core/simulator/SimulatorState.hpp
opm/core/simulator/SimulatorTimerInterface.hpp
opm/core/simulator/SimulatorTimer.hpp
opm/core/simulator/SimulatorTimerInterface.hpp
opm/core/simulator/TimeStepControl.hpp
opm/core/simulator/TimeStepControlInterface.hpp
opm/core/simulator/TwophaseState.hpp
opm/core/simulator/TwophaseState_impl.hpp
opm/core/simulator/WellState.hpp
opm/core/simulator/initState.hpp
opm/core/simulator/initState_impl.hpp
opm/core/simulator/initStateEquil.hpp
opm/core/simulator/initStateEquil_impl.hpp
opm/core/flowdiagnostics/AnisotropicEikonal.hpp
opm/core/flowdiagnostics/DGBasis.hpp
opm/core/flowdiagnostics/FlowDiagnostics.hpp
opm/core/flowdiagnostics/TofReorder.hpp
opm/core/flowdiagnostics/TofDiscGalReorder.hpp
opm/core/simulator/initState_impl.hpp
opm/core/transport/TransportSolverTwophaseInterface.hpp
opm/core/transport/implicit/CSRMatrixBlockAssembler.hpp
opm/core/transport/implicit/CSRMatrixUmfpackSolver.hpp
@@ -411,21 +396,20 @@ list (APPEND PUBLIC_HEADER_FILES
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/SinglePointUpwindTwoPhase.hpp
opm/core/transport/implicit/TransportSolverTwophaseImplicit.hpp
opm/core/transport/implicit/transport_source.h
opm/core/transport/minimal/spu_explicit.h
opm/core/transport/minimal/spu_implicit.h
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp
opm/core/transport/reorder/ReorderSolverInterface.hpp
opm/core/transport/reorder/TransportSolverCompressibleTwophaseReorder.hpp
opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp
opm/core/transport/reorder/reordersequence.h
opm/core/transport/reorder/tarjan.h
opm/core/utility/Average.hpp
opm/core/utility/CompressedPropertyAccess.hpp
opm/core/utility/compressedToCartesian.hpp
opm/core/utility/DataMap.hpp
opm/core/utility/Event.hpp
opm/core/utility/Event_impl.hpp
@@ -443,11 +427,13 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/utility/VelocityInterpolation.hpp
opm/core/utility/WachspressCoord.hpp
opm/core/utility/buildUniformMonotoneTable.hpp
opm/core/utility/compressedToCartesian.hpp
opm/core/utility/extractPvtTableIndex.hpp
opm/core/utility/have_boost_redef.hpp
opm/core/utility/linearInterpolation.hpp
opm/core/utility/miscUtilities.hpp
opm/core/utility/miscUtilities_impl.hpp
opm/core/utility/miscUtilitiesBlackoil.hpp
opm/core/utility/miscUtilities_impl.hpp
opm/core/utility/parameters/Parameter.hpp
opm/core/utility/parameters/ParameterGroup.hpp
opm/core/utility/parameters/ParameterGroup_impl.hpp
@@ -460,6 +446,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/utility/parameters/tinyxml/tinyxml.h
opm/core/utility/share_obj.hpp
opm/core/utility/thresholdPressures.hpp
opm/core/well_controls.h
opm/core/wells.h
opm/core/wells/InjectionSpecification.hpp
opm/core/wells/ProductionSpecification.hpp
opm/core/wells/WellCollection.hpp

406
opm/core/props/BlackoilPropertiesFromDeck.cpp
View File

@@ -20,14 +20,15 @@
#include "config.h"
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/compressedToCartesian.hpp>
#include <opm/core/utility/extractPvtTableIndex.hpp>
#include <vector>
#include <numeric>
namespace Opm
{
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
const UnstructuredGrid& grid,
@@ -133,16 +134,15 @@ namespace Opm
if (init_rock){
rock_.init(eclState, number_of_cells, global_cell, cart_dims);
}
pvt_.init(deck, eclState, /*numSamples=*/0);
phaseUsage_ = phaseUsageFromDeck(deck);
initSurfaceDensities_(deck);
oilPvt_.initFromDeck(deck, eclState);
gasPvt_.initFromDeck(deck, eclState);
waterPvt_.initFromDeck(deck, eclState);
SaturationPropsFromDeck* ptr
= new SaturationPropsFromDeck();
ptr->init(phaseUsageFromDeck(deck), materialLawManager);
satprops_.reset(ptr);
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
inline void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr deck,
@@ -162,8 +162,11 @@ namespace Opm
rock_.init(eclState, number_of_cells, global_cell, cart_dims);
}
const int pvt_samples = param.getDefault("pvt_tab_size", -1);
pvt_.init(deck, eclState, pvt_samples);
phaseUsage_ = phaseUsageFromDeck(deck);
initSurfaceDensities_(deck);
oilPvt_.initFromDeck(deck, eclState);
gasPvt_.initFromDeck(deck, eclState);
waterPvt_.initFromDeck(deck, eclState);
// Unfortunate lack of pointer smartness here...
std::string threephase_model = param.getDefault<std::string>("threephase_model", "gwseg");
@@ -175,18 +178,12 @@ namespace Opm
= new SaturationPropsFromDeck();
ptr->init(phaseUsageFromDeck(deck), materialLawManager);
satprops_.reset(ptr);
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()
{
}
/// \return D, the number of spatial dimensions.
int BlackoilPropertiesFromDeck::numDimensions() const
{
@@ -219,13 +216,13 @@ namespace Opm
/// \return P, the number of phases (also the number of components).
int BlackoilPropertiesFromDeck::numPhases() const
{
return pvt_.numPhases();
return phaseUsage_.num_phases;
}
/// \return Object describing the active phases.
PhaseUsage BlackoilPropertiesFromDeck::phaseUsage() const
{
return pvt_.phaseUsage();
return phaseUsage_;
}
/// \param[in] n Number of data points.
@@ -244,16 +241,45 @@ namespace Opm
double* mu,
double* dmudp) const
{
if (dmudp) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::viscosity() -- derivatives of viscosity not yet implemented.");
} else {
const int *cellPvtTableIdx = cellPvtRegionIndex();
assert(cellPvtTableIdx != 0);
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
const auto& pu = phaseUsage();
pvt_.mu(n, &pvtTableIdx[0], p, T, z, mu);
enum PressureEvalTag {};
typedef Opm::LocalAd::Evaluation<double, PressureEvalTag, /*size=*/1> LadEval;
LadEval pLad = 0.0;
LadEval TLad = 0.0;
LadEval RsLad = 0.0;
LadEval RvLad = 0.0;
LadEval muLad = 0.0;
pLad.derivatives[0] = 1.0;
for (int i = 0; i < n; ++ i) {
int cellIdx = cells[i];
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
pLad.value = p[i];
TLad.value = T[i];
if (pu.phase_used[BlackoilPhases::Aqua]) {
muLad = waterPvt_.viscosity(pvtRegionIdx, TLad, pLad);
int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua];
mu[offset] = muLad.value;
dmudp[offset] = muLad.derivatives[0];
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
muLad = oilPvt_.viscosity(pvtRegionIdx, TLad, pLad, RsLad);
int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid];
mu[offset] = muLad.value;
dmudp[offset] = muLad.derivatives[0];
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
muLad = gasPvt_.viscosity(pvtRegionIdx, TLad, pLad, RvLad);
int offset = pu.num_phases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour];
mu[offset] = muLad.value;
dmudp[offset] = muLad.derivatives[0];
}
}
}
@@ -278,27 +304,23 @@ namespace Opm
{
const int np = numPhases();
const int *cellPvtTableIdx = cellPvtRegionIndex();
std::vector<int> pvtTableIdx(n);
for (int i = 0; i < n; ++ i)
pvtTableIdx[i] = cellPvtTableIdx[cells[i]];
B_.resize(n*np);
R_.resize(n*np);
if (dAdp) {
dB_.resize(n*np);
dR_.resize(n*np);
pvt_.dBdp(n, &pvtTableIdx[0], p, T, z, &B_[0], &dB_[0]);
pvt_.dRdp(n, &pvtTableIdx[0], p, z, &R_[0], &dR_[0]);
this->compute_dBdp_(n, p, T, z, cells, &B_[0], &dB_[0]);
this->compute_dRdp_(n, p, T, z, cells, &R_[0], &dR_[0]);
} else {
pvt_.B(n, &pvtTableIdx[0], p, T, z, &B_[0]);
pvt_.R(n, &pvtTableIdx[0], p, z, &R_[0]);
this->compute_B_(n, p, T, z, cells, &B_[0]);
this->compute_R_(n, p, T, z, cells, &R_[0]);
}
const int* phase_pos = pvt_.phasePosition();
bool oil_and_gas = pvt_.phaseUsed()[BlackoilPhases::Liquid] &&
pvt_.phaseUsed()[BlackoilPhases::Vapour];
const int o = phase_pos[BlackoilPhases::Liquid];
const int g = phase_pos[BlackoilPhases::Vapour];
const auto& pu = phaseUsage();
bool oil_and_gas = pu.phase_pos[BlackoilPhases::Liquid] &&
pu.phase_pos[BlackoilPhases::Vapour];
const int o = pu.phase_pos[BlackoilPhases::Liquid];
const int g = pu.phase_pos[BlackoilPhases::Vapour];
// Compute A matrix
// #pragma omp parallel for
@@ -360,6 +382,276 @@ namespace Opm
}
}
void BlackoilPropertiesFromDeck::compute_B_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* B) const
{
const auto& pu = phaseUsage();
typedef double LadEval;
LadEval pLad = 0.0;
LadEval TLad = 0.0;
LadEval RsLad = 0.0;
LadEval RvLad = 0.0;
for (int i = 0; i < n; ++ i) {
int cellIdx = cells[i];
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
pLad = p[i];
TLad = T[i];
int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua];
if (pu.phase_used[BlackoilPhases::Aqua]) {
LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
B[waterOffset] = BLad;
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
double currentRs = 0.0;
double maxRs = 0.0;
if (pu.phase_used[BlackoilPhases::Vapour]) {
currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad);
}
LadEval BLad;
if (currentRs >= maxRs) {
BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
}
else {
RsLad = currentRs;
BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
}
B[oilOffset] = BLad;
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
double currentRv = 0.0;
double maxRv = 0.0;
if (pu.phase_used[BlackoilPhases::Liquid]) {
currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad);
}
LadEval BLad;
if (currentRv >= maxRv) {
BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
}
else {
RvLad = currentRv;
BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad);
}
B[gasOffset] = BLad;
}
}
}
void BlackoilPropertiesFromDeck::compute_dBdp_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* B,
double* dBdp) const
{
const auto& pu = phaseUsage();
enum PressureEvalTag {};
typedef Opm::LocalAd::Evaluation<double, PressureEvalTag, /*size=*/1> LadEval;
LadEval pLad = 0.0;
LadEval TLad = 0.0;
LadEval RsLad = 0.0;
LadEval RvLad = 0.0;
pLad.derivatives[0] = 1.0;
for (int i = 0; i < n; ++ i) {
int cellIdx = cells[i];
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
pLad.value = p[i];
TLad.value = T[i];
int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua];
if (pu.phase_used[BlackoilPhases::Aqua]) {
LadEval BLad = 1.0/waterPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
B[waterOffset] = BLad.value;
dBdp[waterOffset] = BLad.derivatives[0];
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
double currentRs = 0.0;
double maxRs = 0.0;
if (pu.phase_used[BlackoilPhases::Vapour]) {
currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
maxRs = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad.value, pLad.value);
}
LadEval BLad;
if (currentRs >= maxRs) {
BLad = 1.0/oilPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
}
else {
RsLad.value = currentRs;
BLad = 1.0/oilPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RsLad);
}
B[oilOffset] = BLad.value;
dBdp[oilOffset] = BLad.derivatives[0];
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
double currentRv = 0.0;
double maxRv = 0.0;
if (pu.phase_used[BlackoilPhases::Liquid]) {
currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
maxRv = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad.value, pLad.value);
}
LadEval BLad;
if (currentRv >= maxRv) {
BLad = 1.0/gasPvt_.saturatedInverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad);
}
else {
RvLad.value = currentRv;
BLad = 1.0/gasPvt_.inverseFormationVolumeFactor(pvtRegionIdx, TLad, pLad, RvLad);
}
B[gasOffset] = BLad.value;
dBdp[gasOffset] = BLad.derivatives[0];
}
}
}
void BlackoilPropertiesFromDeck::compute_R_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* R) const
{
const auto& pu = phaseUsage();
typedef double LadEval;
LadEval pLad = 0.0;
LadEval TLad = 0.0;
for (int i = 0; i < n; ++ i) {
int cellIdx = cells[i];
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
pLad = p[i];
TLad = T[i];
int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua];
if (pu.phase_used[BlackoilPhases::Aqua]) {
R[waterOffset] = 0.0; // water is always immiscible!
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad);
double currentRs = 0.0;
if (pu.phase_used[BlackoilPhases::Vapour]) {
currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
}
RsSatLad = std::min(RsSatLad, currentRs);
R[oilOffset] = RsSatLad;
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad);
double currentRv = 0.0;
if (pu.phase_used[BlackoilPhases::Liquid]) {
currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
}
RvSatLad = std::min(RvSatLad, currentRv);
R[gasOffset] = RvSatLad;
}
}
}
void BlackoilPropertiesFromDeck::compute_dRdp_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* R,
double* dRdp) const
{
const auto& pu = phaseUsage();
enum PressureEvalTag {};
typedef Opm::LocalAd::Evaluation<double, PressureEvalTag, /*size=*/1> LadEval;
typedef Opm::MathToolbox<LadEval> Toolbox;
LadEval pLad = 0.0;
LadEval TLad = 0.0;
pLad.derivatives[0] = 1.0;
for (int i = 0; i < n; ++ i) {
int cellIdx = cells[i];
int pvtRegionIdx = cellPvtRegionIdx_[cellIdx];
pLad.value = p[i];
TLad.value = T[i];
int oilOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Liquid];
int gasOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Vapour];
int waterOffset = pu.num_phases*i + pu.phase_pos[BlackoilPhases::Aqua];
if (pu.phase_used[BlackoilPhases::Aqua]) {
R[waterOffset] = 0.0; // water is always immiscible!
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
LadEval RsSatLad = oilPvt_.saturatedGasDissolutionFactor(pvtRegionIdx, TLad, pLad);
LadEval currentRs = 0.0;
if (pu.phase_used[BlackoilPhases::Vapour]) {
currentRs = (z[oilOffset] == 0.0) ? 0.0 : z[gasOffset]/z[oilOffset];
}
RsSatLad = Toolbox::min(RsSatLad, currentRs);
R[oilOffset] = RsSatLad.value;
dRdp[oilOffset] = RsSatLad.derivatives[0];
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
LadEval RvSatLad = gasPvt_.saturatedOilVaporizationFactor(pvtRegionIdx, TLad, pLad);
LadEval currentRv = 0.0;
if (pu.phase_used[BlackoilPhases::Liquid]) {
currentRv = (z[gasOffset] == 0.0) ? 0.0 : z[oilOffset]/z[gasOffset];
}
RvSatLad = Toolbox::min(RvSatLad, currentRv);
R[gasOffset] = RvSatLad.value;
dRdp[gasOffset] = RvSatLad.derivatives[0];
}
}
}
/// \param[in] n Number of data points.
/// \param[in] A Array of nP^2 values, where the P^2 values for a cell give the
/// matrix A = RB^{-1} which relates z to u by z = Au. The matrices
@@ -376,8 +668,7 @@ namespace Opm
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int cellIdx = cells?cells[i]:i;
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
const double* sdens = pvt_.surfaceDensities(pvtRegionIdx);
const double *sdens = surfaceDensity(cellIdx);
for (int phase = 0; phase < np; ++phase) {
rho[np*i + phase] = 0.0;
for (int comp = 0; comp < np; ++comp) {
@@ -391,8 +682,37 @@ namespace Opm
/// \return Array of P density values.
const double* BlackoilPropertiesFromDeck::surfaceDensity(int cellIdx) const
{
const auto& pu = phaseUsage();
int pvtRegionIdx = getTableIndex_(cellPvtRegionIndex(), cellIdx);
return pvt_.surfaceDensities(pvtRegionIdx);
return &surfaceDensities_[pvtRegionIdx*pu.num_phases];
}
void BlackoilPropertiesFromDeck::initSurfaceDensities_(Opm::DeckConstPtr deck)
{
const auto& pu = phaseUsage();
int np = pu.num_phases;
int numPvtRegions = 1;
if (deck->hasKeyword("TABDIMS")) {
const auto& tabdimsKeyword = deck->getKeyword("TABDIMS");
numPvtRegions = tabdimsKeyword.getRecord(0).getItem("NTPVT").template get<int>(0);
}
const auto& densityKeyword = deck->getKeyword("DENSITY");
surfaceDensities_.resize(np*numPvtRegions);
for (int pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) {
if (pu.phase_used[BlackoilPhases::Aqua])
surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Aqua]] =
densityKeyword.getRecord(pvtRegionIdx).getItem("WATER").getSIDouble(0);
if (pu.phase_used[BlackoilPhases::Liquid])
surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Liquid]] =
densityKeyword.getRecord(pvtRegionIdx).getItem("OIL").getSIDouble(0);
if (pu.phase_used[BlackoilPhases::Vapour])
surfaceDensities_[np*pvtRegionIdx + pu.phase_pos[BlackoilPhases::Vapour]] =
densityKeyword.getRecord(pvtRegionIdx).getItem("GAS").getSIDouble(0);
}
}
/// \param[in] n Number of data points.

57
opm/core/props/BlackoilPropertiesFromDeck.hpp
View File

@@ -23,9 +23,11 @@
#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/utility/parameters/ParameterGroup.hpp>
#include <opm/material/fluidsystems/blackoilpvt/OilPvtMultiplexer.hpp>
#include <opm/material/fluidsystems/blackoilpvt/GasPvtMultiplexer.hpp>
#include <opm/material/fluidsystems/blackoilpvt/WaterPvtMultiplexer.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
@@ -233,10 +235,19 @@ namespace Opm
const double pcow,
double & swat);
// return a reference to the "raw" PVT fluid object for a phase.
const PvtInterface& pvt(int phaseIdx) const
const OilPvtMultiplexer<double>& oilPvt() const
{
return pvt_.pvt(phaseIdx);
return oilPvt_;
}
const GasPvtMultiplexer<double>& gasPvt() const
{
return gasPvt_;
}
const WaterPvtMultiplexer<double>& waterPvt() const
{
return waterPvt_;
}
private:
@@ -247,6 +258,38 @@ namespace Opm
return pvtTableIdx[cellIdx];
}
void initSurfaceDensities_(Opm::DeckConstPtr deck);
void compute_B_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* B) const;
void compute_dBdp_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* B,
double* dBdp) const;
void compute_R_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* R) const;
void compute_dRdp_(const int n,
const double* p,
const double* T,
const double* z,
const int* cells,
double* R,
double* dRdp) const;
void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclState,
std::shared_ptr<MaterialLawManager> materialLawManager,
@@ -265,10 +308,14 @@ namespace Opm
bool init_rock);
RockFromDeck rock_;
PhaseUsage phaseUsage_;
std::vector<int> cellPvtRegionIdx_;
BlackoilPvtProperties pvt_;
OilPvtMultiplexer<double> oilPvt_;
GasPvtMultiplexer<double> gasPvt_;
WaterPvtMultiplexer<double> waterPvt_;
std::shared_ptr<MaterialLawManager> materialLawManager_;
std::shared_ptr<SaturationPropsInterface> satprops_;
std::vector<double> surfaceDensities_;
mutable std::vector<double> B_;
mutable std::vector<double> dB_;
mutable std::vector<double> R_;

2
opm/core/props/IncompPropertiesFromDeck.hpp
View File

@@ -23,8 +23,8 @@
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.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/rock/RockFromDeck.hpp>
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
struct UnstructuredGrid;

255
opm/core/props/pvt/BlackoilPvtProperties.cpp
View File

@@ -1,255 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/props/pvt/PvtDead.hpp>
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
#include <opm/core/props/pvt/PvtLiveOil.hpp>
#include <opm/core/props/pvt/PvtLiveGas.hpp>
#include <opm/core/props/pvt/PvtConstCompr.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
namespace Opm
{
BlackoilPvtProperties::BlackoilPvtProperties()
{
}
void BlackoilPvtProperties::init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int numSamples)
{
phase_usage_ = phaseUsageFromDeck(deck);
// Surface densities. Accounting for different orders in eclipse and our code.
const auto& densityKeyword = deck->getKeyword("DENSITY");
int numRegions = densityKeyword.size();
densities_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
if (phase_usage_.phase_used[Liquid]) {
densities_[regionIdx][phase_usage_.phase_pos[Liquid]]
= densityKeyword.getRecord(regionIdx).getItem("OIL").getSIDouble(0);
}
if (phase_usage_.phase_used[Aqua]) {
densities_[regionIdx][phase_usage_.phase_pos[Aqua]]
= densityKeyword.getRecord(regionIdx).getItem("WATER").getSIDouble(0);
}
if (phase_usage_.phase_used[Vapour]) {
densities_[regionIdx][phase_usage_.phase_pos[Vapour]]
= densityKeyword.getRecord(regionIdx).getItem("GAS").getSIDouble(0);
}
}
// Resize the property objects container
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
// if water is used, we require the presence of the "PVTW"
// keyword for now...
std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
pvtw->initFromWater(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
}
{
auto tables = eclipseState->getTableManager();
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
// for oil, we support the "PVDO", "PVTO" and "PVCDO"
// keywords...
const auto &pvdoTables = tables->getPvdoTables();
const auto &pvtoTables = tables->getPvtoTables();
if (pvdoTables.size() > 0) {
if (numSamples > 0) {
auto splinePvt = std::shared_ptr<PvtDeadSpline>(new PvtDeadSpline);
splinePvt->initFromOil(pvdoTables, numSamples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
auto deadPvt = std::shared_ptr<PvtDead>(new PvtDead);
deadPvt->initFromOil(pvdoTables);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (pvtoTables.size() > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(pvtoTables));
} else if (deck->hasKeyword("PVCDO")) {
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
}
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
// gas can be specified using the "PVDG" or "PVTG" keywords...
const auto &pvdgTables = tables->getPvdgTables();
const auto &pvtgTables = tables->getPvtgTables();
if (pvdgTables.size() > 0) {
if (numSamples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgTables, numSamples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgTables);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (pvtgTables.size() > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(pvtgTables));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
}
}
}
const double* BlackoilPvtProperties::surfaceDensities(int regionIdx) const
{
return &densities_[regionIdx][0];
}
PhaseUsage BlackoilPvtProperties::phaseUsage() const
{
return phase_usage_;
}
int BlackoilPvtProperties::numPhases() const
{
return phase_usage_.num_phases;
}
const int* BlackoilPvtProperties::phaseUsed() const
{
return phase_usage_.phase_used;
}
const int* BlackoilPvtProperties::phasePosition() const
{
return phase_usage_.phase_pos;
}
void BlackoilPvtProperties::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->mu(n, pvtTableIdx, p, T, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_mu[phase_usage_.num_phases*i + phase] = data1_[i];
}
}
}
void BlackoilPvtProperties::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->B(n, pvtTableIdx, p, T, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
}
}
}
void BlackoilPvtProperties::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const
{
data1_.resize(n);
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dBdp(n, pvtTableIdx, p, T, z, &data1_[0], &data2_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[phase_usage_.num_phases*i + phase] = data1_[i];
output_dBdp[phase_usage_.num_phases*i + phase] = data2_[i];
}
}
}
void BlackoilPvtProperties::R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const
{
data1_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->R(n, pvtTableIdx, p, z, &data1_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
}
}
}
void BlackoilPvtProperties::dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
data1_.resize(n);
data2_.resize(n);
for (int phase = 0; phase < phase_usage_.num_phases; ++phase) {
props_[phase]->dRdp(n, pvtTableIdx, p, z, &data1_[0], &data2_[0]);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[phase_usage_.num_phases*i + phase] = data1_[i];
output_dRdp[phase_usage_.num_phases*i + phase] = data2_[i];
}
}
}
} // namespace Opm

146
opm/core/props/pvt/BlackoilPvtProperties.hpp
View File

@@ -1,146 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
#define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <string>
#include <memory>
#include <array>
namespace Opm
{
/// Class collecting the pvt properties for all active phases.
/// For all the methods, the following apply:
/// - p and z are expected to be of size n and n*num_phases, respectively.
/// - pvtTableIdx specifies the PVT table to be used for each data
/// point and is thus expected to be an array of size n
/// - Output arrays shall be of size n*num_phases, and must be valid
/// before calling the method.
/// NOTE: The difference between this interface and the one defined
/// by PvtInterface is that this collects all phases' properties,
/// and therefore the output arrays are of size n*num_phases as opposed
/// to size n in PvtInterface.
class BlackoilPvtProperties : public BlackoilPhases
{
public:
/// Default constructor.
BlackoilPvtProperties();
/// Initialize from deck.
///
/// \param deck An input deck from the opm-parser module.
void init(Opm::DeckConstPtr deck,
Opm::EclipseStateConstPtr eclipseState,
int samples);
/// \return Object describing the active phases.
PhaseUsage phaseUsage() const;
/// Number of active phases.
int numPhases() const;
/// For each canonical phase, indicates if it is
/// active or not (boolean usage of int).
/// \return Array of size MaxNumPhases
const int* phaseUsed() const;
/// Positions of canonical phases in arrays of phase
/// properties (saturations etc.).
/// \return Array of size MaxNumPhases
const int* phasePosition() const;
/// Densities of stock components at surface conditions.
/// \return Array of size numPhases().
const double* surfaceDensities(int regionIdx = 0) const;
/// Viscosity as a function of p, T and z.
void mu(const int n,
const int *pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const;
/// Formation volume factor as a function of p, T and z.
void B(const int n,
const int *pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p, T and z.
void dBdp(const int n,
const int *pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const;
/// Solution factor as a function of p and z.
void R(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
void dRdp(const int n,
const int *pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
// return a reference to the "raw" PVT fluid object for a phase.
const PvtInterface& pvt(int phaseIdx) const
{
return *props_[phaseIdx];
}
private:
// Disabling copying (just to avoid surprises, since we use shared_ptr).
BlackoilPvtProperties(const BlackoilPvtProperties&);
BlackoilPvtProperties& operator=(const BlackoilPvtProperties&);
PhaseUsage phase_usage_;
// The PVT properties. We need to store one object per PVT
// region per active fluid phase.
std::vector<std::shared_ptr<PvtInterface> > props_;
std::vector<std::array<double, MaxNumPhases> > densities_;
mutable std::vector<double> data1_;
mutable std::vector<double> data2_;
};
}
#endif // OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED

314
opm/core/props/pvt/PvtConstCompr.hpp
View File

@@ -1,314 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PVTCONSTCOMPR_HEADER_INCLUDED
#define OPM_PVTCONSTCOMPR_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/Deck/DeckItem.hpp>
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
#include <vector>
#include <algorithm>
namespace Opm
{
/// Class for constant compressible phases (PVTW or PVCDO).
/// The PVT properties can either be given as a function of
/// pressure (p) and surface volume (z) or pressure (p) and gas
/// resolution factor (r). Also, since this class supports
/// multiple PVT regions, the concrete table to be used for each
/// data point needs to be specified via the pvtTableIdx argument
/// of the respective method. For all the virtual methods, the
/// following apply: pvtTableIdx, p, r and z are expected to be of
/// size n, size n, size n and n*num_phases, respectively. Output
/// arrays shall be of size n, and must be valid before calling
/// the method.
class PvtConstCompr : public PvtInterface
{
public:
PvtConstCompr()
{}
void initFromWater(const Opm::DeckKeyword& pvtwKeyword)
{
int numRegions = pvtwKeyword.size();
ref_press_.resize(numRegions);
ref_B_.resize(numRegions);
comp_.resize(numRegions);
viscosity_.resize(numRegions);
visc_comp_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
const auto& pvtwRecord = pvtwKeyword.getRecord(regionIdx);
ref_press_[regionIdx] = pvtwRecord.getItem("P_REF").getSIDouble(0);
ref_B_[regionIdx] = pvtwRecord.getItem("WATER_VOL_FACTOR").getSIDouble(0);
comp_[regionIdx] = pvtwRecord.getItem("WATER_COMPRESSIBILITY").getSIDouble(0);
viscosity_[regionIdx] = pvtwRecord.getItem("WATER_VISCOSITY").getSIDouble(0);
visc_comp_[regionIdx] = pvtwRecord.getItem("WATER_VISCOSIBILITY").getSIDouble(0);
}
}
void initFromOil(const DeckKeyword& pvcdoKeyword)
{
int numRegions = pvcdoKeyword.size();
ref_press_.resize(numRegions);
ref_B_.resize(numRegions);
comp_.resize(numRegions);
viscosity_.resize(numRegions);
visc_comp_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
const auto& pvcdoRecord = pvcdoKeyword.getRecord(regionIdx);
ref_press_[regionIdx] = pvcdoRecord.getItem("P_REF").getSIDouble(0);
ref_B_[regionIdx] = pvcdoRecord.getItem("OIL_VOL_FACTOR").getSIDouble(0);
comp_[regionIdx] = pvcdoRecord.getItem("OIL_COMPRESSIBILITY").getSIDouble(0);
viscosity_[regionIdx] = pvcdoRecord.getItem("OIL_VISCOSITY").getSIDouble(0);
visc_comp_[regionIdx] = pvcdoRecord.getItem("OIL_VISCOSIBILITY").getSIDouble(0);
}
}
/*!
* \brief Create a PVT object with a given viscosity that
* assumes all fluid phases to be incompressible.
*/
explicit PvtConstCompr(double visc)
: ref_press_(1, 0.0),
ref_B_(1, 1.0),
comp_(1, 0.0),
viscosity_(1, visc),
visc_comp_(1, 0.0)
{
}
virtual ~PvtConstCompr()
{
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
}
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
output_mu[i] = viscosity_[tableIdx]/d;
output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
output_mu[i] = viscosity_[tableIdx]/d;
output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
output_B[i] = ref_B_[tableIdx]/(1.0 + x + 0.5*x*x);
}
}
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
output_B[i] = ref_B_[tableIdx]/d;
output_dBdp[i] = (-ref_B_[tableIdx]/(d*d))*(1 + x) * comp_[tableIdx];
}
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
// b = 1/B = d/ref_B_B;
output_b[i] = d/ref_B_[tableIdx];
output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
// Computing a polynomial approximation to the exponential.
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double x = comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
double d = (1.0 + x + 0.5*x*x);
// b = 1/B = d/ref_B_[tableIdx]B;
output_b[i] = d/ref_B_[tableIdx];
output_dbdp[i] = (1 + x) * comp_[tableIdx]/ref_B_[tableIdx];
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
virtual void rsSat(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
virtual void rvSat(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
virtual void R(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
virtual void dRdp(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// The PVT properties. We need to store one value per PVT
// region.
std::vector<double> ref_press_;
std::vector<double> ref_B_;
std::vector<double> comp_;
std::vector<double> viscosity_;
std::vector<double> visc_comp_;
};
}
#endif // OPM_PVTCONSTCOMPR_HEADER_INCLUDED

294
opm/core/props/pvt/PvtDead.cpp
View File

@@ -1,294 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtDead.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdoTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdgTable.hpp>
#include <algorithm>
// Extra includes for debug dumping of tables.
// #include <boost/lexical_cast.hpp>
// #include <string>
// #include <fstream>
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
/// Constructor
void PvtDead::initFromOil(const TableContainer& pvdoTables)
{
int numRegions = pvdoTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
inverseBmu_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdoTable& pvdoTable = pvdoTables.getTable<PvdoTable>(regionIdx);
// Copy data
const auto& press = pvdoTable.getColumn("P");
const auto& b_var = pvdoTable.getColumn("BO");
const auto& visc = pvdoTable.getColumn("MUO");
const int sz = b_var.size();
std::vector<double> inverseB(sz);
for (int i = 0; i < sz; ++i) {
inverseB[i] = 1.0 / b_var[i];
}
std::vector<double> inverseBmu(sz);
for (int i = 0; i < sz; ++i) {
inverseBmu[i] = 1.0 / (b_var[i] * visc[i]);
}
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
}
}
void PvtDead::initFromGas(const TableContainer& pvdgTables)
{
int numRegions = pvdgTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
inverseBmu_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdgTable& pvdgTable = pvdgTables.getTable<PvdgTable>(regionIdx);
// Copy data
const auto& press = pvdgTable.getColumn("P");
const auto& b = pvdgTable.getColumn("BG");
const auto& visc = pvdgTable.getColumn("MUG");
const int sz = b.size();
std::vector<double> inverseB(sz);
for (int i = 0; i < sz; ++i) {
inverseB[i] = 1.0 / b[i];
}
std::vector<double> inverseBmu(sz);
for (int i = 0; i < sz; ++i) {
inverseBmu[i] = 1.0 / (b[i] * visc[i]);
}
b_[regionIdx] = NonuniformTableLinear<double>(press, inverseB);
viscosity_[regionIdx] = NonuniformTableLinear<double>(press, visc);
inverseBmu_[regionIdx] = NonuniformTableLinear<double>(press, inverseBmu);
}
}
// Destructor
PvtDead::~PvtDead()
{
}
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
double tempInvB = b_[regionIdx](p[i]);
double tempInvBmu = inverseBmu_[regionIdx](p[i]);
output_mu[i] = tempInvB / tempInvBmu;
}
}
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
double tempInvB = b_[regionIdx](p[i]);
double tempInvBmu = inverseBmu_[regionIdx](p[i]);
output_mu[i] = tempInvB / tempInvBmu;
output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
- tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void PvtDead::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
double tempInvB = b_[regionIdx](p[i]);
double tempInvBmu = inverseBmu_[regionIdx](p[i]);
output_mu[i] = tempInvB / tempInvBmu;
output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
- tempInvB * inverseBmu_[regionIdx].derivative(p[i]))
/ (tempInvBmu * tempInvBmu);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void PvtDead::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_B) const
{
// #pragma omp parallel for
// B = 1/b
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = 1.0/b_[regionIdx](p[i]);
}
}
void PvtDead::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, pvtTableIdx, p, T, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
}
}
void PvtDead::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void PvtDead::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void PvtDead::rsSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
void PvtDead::rvSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
void PvtDead::R(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
void PvtDead::dRdp(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
}

171
opm/core/props/pvt/PvtDead.hpp
View File

@@ -1,171 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PVTDEAD_HEADER_INCLUDED
#define OPM_PVTDEAD_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableContainer.hpp>
#include <vector>
namespace Opm
{
/// Class for immiscible dead oil and dry gas.
/// The PVT properties can either be given as a function of
/// pressure (p) and surface volume (z) or pressure (p) and gas
/// resolution factor (r). Also, since this class supports
/// multiple PVT regions, the concrete table to be used for each
/// data point needs to be specified via the pvtTableIdx argument
/// of the respective method. For all the virtual methods, the
/// following apply: pvtTableIdx, p, r and z are expected to be of
/// size n, size n, size n and n*num_phases, respectively. Output
/// arrays shall be of size n, and must be valid before calling
/// the method.
class PvtDead : public PvtInterface
{
public:
PvtDead() {};
void initFromOil(const TableContainer& pvdoTables);
void initFromGas(const TableContainer& pvdgTables);
virtual ~PvtDead();
/// Viscosity as a function of p, T and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const;
/// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p, T and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// PVT properties of dry gas or dead oil. We need to store one
// table per PVT region.
std::vector<NonuniformTableLinear<double> > b_;
std::vector<NonuniformTableLinear<double> > viscosity_;
std::vector<NonuniformTableLinear<double> > inverseBmu_;
};
}
#endif // OPM_PVTDEAD_HEADER_INCLUDED

268
opm/core/props/pvt/PvtDeadSpline.cpp
View File

@@ -1,268 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdoTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvdgTable.hpp>
#include <algorithm>
// Extra includes for debug dumping of tables.
// #include <boost/lexical_cast.hpp>
// #include <string>
// #include <fstream>
namespace Opm
{
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtDeadSpline::PvtDeadSpline()
{}
void PvtDeadSpline::initFromOil(const TableContainer& pvdoTables,
int numSamples)
{
int numRegions = pvdoTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdoTable& pvdoTable = pvdoTables.getTable<PvdoTable>(regionIdx);
int numRows = pvdoTable.numRows();
// Copy data
std::vector<double> press = pvdoTable.getColumn("P").vectorCopy();
std::vector<double> B_var = pvdoTable.getColumn("BO").vectorCopy();
std::vector<double> visc = pvdoTable.getColumn("MUO").vectorCopy();
std::vector<double> B_inv(numRows);
for (int i = 0; i < numRows; ++i) {
B_inv[i] = 1.0 / B_var[i];
}
buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
}
}
void PvtDeadSpline::initFromGas(const TableContainer& pvdgTables,
int numSamples)
{
int numRegions = pvdgTables.size();
// resize the attributes of the object
b_.resize(numRegions);
viscosity_.resize(numRegions);
for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
const Opm::PvdgTable& pvdgTable = pvdgTables.getTable<PvdgTable>(regionIdx);
int numRows = pvdgTable.numRows();
// Copy data
std::vector<double> press = pvdgTable.getColumn("P").vectorCopy();
std::vector<double> B = pvdgTable.getColumn("BG").vectorCopy();
std::vector<double> visc = pvdgTable.getColumn("MUG").vectorCopy();
std::vector<double> B_inv(numRows);
for (int i = 0; i < numRows; ++i) {
B_inv[i] = 1.0 / B[i];
}
buildUniformMonotoneTable(press, B_inv, numSamples, b_[regionIdx]);
buildUniformMonotoneTable(press, visc, numSamples, viscosity_[regionIdx]);
}
}
// Destructor
PvtDeadSpline::~PvtDeadSpline()
{
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
}
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void PvtDeadSpline::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_mu[i] = viscosity_[regionIdx](p[i]);
output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
}
std::fill(output_dmudr, output_dmudr + n, 0.0);
}
void PvtDeadSpline::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* /*z*/,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = 1.0/b_[regionIdx](p[i]);
}
}
void PvtDeadSpline::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* /*z*/,
double* output_B,
double* output_dBdp) const
{
B(n, pvtTableIdx, p, T, 0, output_B);
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
double Bg = output_B[i];
output_dBdp[i] = -Bg*Bg*b_[regionIdx].derivative(p[i]);
}
}
void PvtDeadSpline::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /* T */,
const double* /*r*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void PvtDeadSpline::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /* T */,
const double* /*r*/,
const PhasePresence* /*cond*/,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int regionIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = b_[regionIdx](p[i]);
output_dbdp[i] = b_[regionIdx].derivative(p[i]);
}
std::fill(output_dbdr, output_dbdr + n, 0.0);
}
void PvtDeadSpline::rsSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
void PvtDeadSpline::rvSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
void PvtDeadSpline::R(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R) const
{
std::fill(output_R, output_R + n, 0.0);
}
void PvtDeadSpline::dRdp(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
const double* /*z*/,
double* output_R,
double* output_dRdp) const
{
std::fill(output_R, output_R + n, 0.0);
std::fill(output_dRdp, output_dRdp + n, 0.0);
}
}

169
opm/core/props/pvt/PvtDeadSpline.hpp
View File

@@ -1,169 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PVTDEADSPLINE_HEADER_INCLUDED
#define OPM_PVTDEADSPLINE_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableContainer.hpp>
#include <vector>
namespace Opm
{
/// Class for immiscible dead oil and dry gas.
/// The PVT properties can either be given as a function of pressure (p) and surface volume (z)
/// or pressure (p) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class PvtDeadSpline : public PvtInterface
{
public:
PvtDeadSpline();
void initFromOil(const TableContainer& pvdoTables,
int numSamples);
void initFromGas(const TableContainer& pvdgTables,
int numSamples);
virtual ~PvtDeadSpline();
/// Viscosity as a function of p, T and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const;
/// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p, T and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
// PVT properties of dry gas or dead oil. We need to store one
// table per PVT region.
std::vector<UniformTableLinear<double> > b_;
std::vector<UniformTableLinear<double> > viscosity_;
};
}
#endif // OPM_PVTDEADSPLINE_HEADER_INCLUDED

67
opm/core/props/pvt/PvtInterface.cpp
View File

@@ -1,67 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
#include <cassert>
namespace Opm
{
PvtInterface::PvtInterface()
: num_phases_(MaxNumPhases)
{
for (int i = 0; i < MaxNumPhases; ++i) {
phase_pos_[i] = i;
}
}
PvtInterface::~PvtInterface()
{
}
void PvtInterface::setPhaseConfiguration(const int num_phases, const int* phase_pos)
{
num_phases_ = num_phases;
for (int i = 0; i < MaxNumPhases; ++i) {
phase_pos_[i] = phase_pos[i];
}
}
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int *compressedToCartesianCellIdx)
{
//Get the PVTNUM data
const std::vector<int>& pvtnumData = eclState->getIntGridProperty("PVTNUM")->getData();
// Convert this into an array of compressed cells
// Eclipse uses Fortran-style indices which start at 1
// instead of 0, we subtract 1.
pvtTableIdx.resize(numCompressed);
for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) {
size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx;
assert(cartesianCellIdx < pvtnumData.size());
pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1;
}
}
} // namespace Opm

187
opm/core/props/pvt/PvtInterface.hpp
View File

@@ -1,187 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PVTINTERFACE_HEADER_INCLUDED
#define OPM_PVTINTERFACE_HEADER_INCLUDED
#include <vector>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
namespace Opm
{
class PvtInterface : public BlackoilPhases
{
public:
PvtInterface();
virtual ~PvtInterface();
/// \param[in] num_phases The number of active phases.
/// \param[in] phase_pos Array of BlackpoilPhases::MaxNumPhases
/// integers, giving the relative
/// positions of the three canonical
/// phases A, L, V in order to handle
/// arbitrary two-phase and three-phase situations.
void setPhaseConfiguration(const int num_phases, const int* phase_pos);
/// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply:
/// - pvtRegionIdx is an array of size n and represents the
/// index of the PVT table which should be used to calculate
/// the output. NULL can also be passed and is interpreted
/// such that the first table should be used for the output
/// - p, r and z are expected to be of size n, size n and
/// n*num_phases, respectively.
/// - Output arrays shall be of size n, and must be valid before
/// calling the method.
/// Viscosity as a function of p, T and z.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const = 0;
/// Viscosity as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const = 0;
/// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const = 0;
/// Formation volume factor as a function of p, T and z.
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const = 0;
/// Formation volume factor and p-derivative as functions of p and z.
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const = 0;
/// The inverse of the volume factor b = 1 / B as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dpdr) const = 0;
/// The inverse of the volume factor b = 1 / B as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dpdr) const = 0;
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const = 0;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const = 0;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const = 0;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const = 0;
protected:
int num_phases_;
int phase_pos_[MaxNumPhases];
};
/*!
* \brief Helper function to create an array containing the (C-Style)
* PVT table index for each compressed cell from an Eclipse deck.
*
* This function assumes that the degrees of freedom where PVT
* properties need to be calculated are grid cells. The main point
* of this function is to avoid code duplication because the
* Eclipse deck only contains Fortran-style PVT table indices
* which start at 1 instead of 0 and -- more relevantly -- it uses
* logically cartesian cell indices to specify the table index of
* a cell while the classes which use the PvtInterface
* implementations usually use compressed cells.
*/
void extractPvtTableIndex(std::vector<int>& pvtTableIdx,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int* compressedToCartesianIdx);
} // namespace Opm
#endif // OPM_PVTINTERFACE_HEADER_INCLUDED

613
opm/core/props/pvt/PvtLiveGas.cpp
View File

@@ -1,613 +0,0 @@
//===========================================================================
//
// File: MiscibilityLiveGas.cpp
//
// Created: Wed Feb 10 09:21:53 2010
//
// Author: Bjørn Spjelkavik <bsp@sintef.no>
//
// Revision: $Id$
//
//===========================================================================
/*
Copyright 2010 SINTEF ICT, Applied Mathematics.
Copyright 2015 IRIS AS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtLiveGas.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvtgTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
#include <algorithm>
#include <iostream>
namespace Opm
{
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveGas::PvtLiveGas(const std::vector<Opm::PvtgTable>& pvtgTables)
{
int numTables = pvtgTables.size();
saturated_gas_table_.resize(numTables);
undersat_gas_tables_.resize(numTables);
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
const Opm::PvtgTable& pvtgTable = pvtgTables[pvtTableIdx];
const auto& saturatedPvtg = pvtgTable.getSaturatedTable( );
// GAS, P
// saturated_gas_table_[pvtTableIdx].resize(4);
// Adding one extra line to PVTG to store 1./(Bg*mu_g)
saturated_gas_table_[pvtTableIdx].resize(5);
for (int k=0; k<5; ++k) {
saturated_gas_table_[pvtTableIdx][k].resize(saturatedPvtg.numRows());
}
for (size_t row=0; row < saturatedPvtg.numRows(); row++) {
saturated_gas_table_[pvtTableIdx][0][row] = saturatedPvtg.get("PG" , row); // Pressure
saturated_gas_table_[pvtTableIdx][1][row] = saturatedPvtg.get("BG" , row); // Bg
saturated_gas_table_[pvtTableIdx][2][row] = saturatedPvtg.get("MUG" , row); // mu_g
// 1/Bg will go in [3]
saturated_gas_table_[pvtTableIdx][4][row] = saturatedPvtg.get("RV" , row); // Rv
}
int sz = saturatedPvtg.numRows();
undersat_gas_tables_[pvtTableIdx].resize(sz);
for (int i=0; i<sz; ++i) {
const auto &undersatTable = pvtgTable.getUnderSaturatedTable(i);
undersat_gas_tables_[pvtTableIdx][i].resize(4);
for (size_t j=0; j < 4; j++) {
undersat_gas_tables_[pvtTableIdx][i][j].resize( undersatTable.numRows() );
}
for (size_t row=0; row < undersatTable.numRows(); row++) {
undersat_gas_tables_[pvtTableIdx][i][0][row] = undersatTable.get("RV" , row); // Rv
undersat_gas_tables_[pvtTableIdx][i][1][row] = undersatTable.get("BG" , row); // Bg
undersat_gas_tables_[pvtTableIdx][i][2][row] = undersatTable.get("MUG" , row); // mu_g
}
}
// Bg -> 1/Bg
for (int i=0; i<sz; ++i) {
saturated_gas_table_[pvtTableIdx][3][i] = 1.0 / (saturated_gas_table_[pvtTableIdx][1][i]
* saturated_gas_table_[pvtTableIdx][2][i]); // 1/(Bg*mu_g)
saturated_gas_table_[pvtTableIdx][1][i] = 1.0 / saturated_gas_table_[pvtTableIdx][1][i];
for (size_t j=0; j<undersat_gas_tables_[pvtTableIdx][i][1].size(); ++j) {
undersat_gas_tables_[pvtTableIdx][i][3][j] = 1.0 / (undersat_gas_tables_[pvtTableIdx][i][1][j]
* undersat_gas_tables_[pvtTableIdx][i][2][j]); // 1/(Bg*mu_g)
undersat_gas_tables_[pvtTableIdx][i][1][j] = 1.0 / undersat_gas_tables_[pvtTableIdx][i][1][j];
}
}
// Complete undersaturated tables by extrapolating from existing data
int iNext = -1;
for (int i = 0; i < sz; ++i) {
// Skip records already containing undersaturated data
if (undersat_gas_tables_[pvtTableIdx][i][0].size() > 1) {
continue;
}
// Look ahead for next record containing undersaturated data
if (iNext < i) {
iNext = i+1;
while (iNext < sz && undersat_gas_tables_[pvtTableIdx][iNext][0].size() < 2) {
++iNext;
}
if (iNext == sz) {
OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
}
}
// Undersaturated data is added to the current record in the same way as for liveoil.
// It is unclear whether this expantion maintains the compressibility and viscosibility,
// but it seems like it reproduces eclipse results for spe3.
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
auto& from_table = undersat_gas_tables_[pvtTableIdx][iNext];
auto& to_table = undersat_gas_tables_[pvtTableIdx][i];
enum {RV = 0, BG = 1, MUG = 2, BGMUG = 3};
for (sz_t j = 1; j < from_table[0].size(); ++j) {
double diffSolubility = from_table[RV][j] - from_table[RV][j-1];
double solubility = to_table[RV].back() + diffSolubility;
to_table[RV].push_back(solubility);
double compr = (1.0/from_table[BG][j] - 1.0/from_table[BG][j-1])
/ (0.5*(1.0/from_table[BG][j] + 1.0/from_table[BG][j-1]));
double B_var = (1.0/to_table[BG].back()) * (1.0+0.5*compr) / (1.0-0.5*compr);
to_table[BG].push_back(1.0/B_var);
double visc = (from_table[MUG][j] - from_table[MUG][j-1])
/ (0.5*(from_table[MUG][j] + from_table[MUG][j-1]));
double mu_var = (to_table[MUG].back()) * (1.0+0.5*visc) / (1.0-0.5*visc);
to_table[MUG].push_back(mu_var);
// A try to expolate the 1/BMu with the expolated mu and B
double inverseBMu = 1.0 / (B_var*mu_var);
to_table[BGMUG].push_back(inverseBMu);
}
}
}
}
// Destructor
PvtLiveGas::~PvtLiveGas()
{
}
void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
double inverseB = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 1, false);
double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false);
output_mu[i] = inverseB / inverseBMu;
}
}
/// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::mu(const int /*n*/,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*T*/,
const double* /*r*/,
double* /*output_mu*/,
double* /*output_dmudp*/,
double* /*output_dmudr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
int tableIdx = getTableIndex_(pvtRegionIdx, i);
double inverseBMu = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 0);
double inverseB = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0);
output_mu[i] = inverseB / inverseBMu;
double dinverseBmudp = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 1);
double dinverseBdp = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1);
output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
/ (inverseBMu * inverseBMu);
double dinverseBmudr = miscible_gas(p[i], r[i], cnd, tableIdx, 3, 2);
double dinverseBdr = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2);
output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
/ (inverseBMu * inverseBMu);
}
}
/// Formation volume factor as a function of p, T and z.
void PvtLiveGas::B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_B[i] = evalB(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
}
}
/// Formation volume factor and p-derivative as functions of p and z.
void PvtLiveGas::dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalBDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_B[i], output_dBdp[i]);
}
}
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::b(const int /*n*/,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
const double* /*T*/,
const double* /*r*/,
double* /*output_b*/,
double* /*output_dbdp*/,
double* /*output_dbdr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
void PvtLiveGas::b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
int tableIdx = getTableIndex_(pvtRegionIdx, i);
output_b[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 0);
output_dbdp[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 1);
output_dbdr[i] = miscible_gas(p[i], r[i], cnd, tableIdx, 1, 2);
}
}
/// Gas resolution and its derivatives at bublepoint as a function of p.
void PvtLiveGas::rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const
{
for (int i = 0; i < n; ++i) {
int pvtTableIdx = getTableIndex_(pvtRegionIdx, i);
output_rvSat[i] = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][4],p[i]);
output_drvSatdp[i] = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][4],p[i]);
}
}
void PvtLiveGas::rsSat(const int n,
const int* /*pvtRegionIdx*/,
const double* /*p*/,
double* output_rsSat,
double* output_drsSatdp) const
{
std::fill(output_rsSat, output_rsSat + n, 0.0);
std::fill(output_drsSatdp, output_drsSatdp + n, 0.0);
}
/// Solution factor as a function of p and z.
void PvtLiveGas::R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
output_R[i] = evalR(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i));
}
}
/// Solution factor and p-derivative as functions of p and z.
void PvtLiveGas::dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
evalRDeriv(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), output_R[i], output_dRdp[i]);
}
}
// ---- Private methods ----
double PvtLiveGas::evalB(const double press, const double* surfvol, int pvtTableIdx) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
// To handle no-gas case.
return 1.0;
}
return 1.0/miscible_gas(press, surfvol, pvtTableIdx, 1, false);
}
void PvtLiveGas::evalBDeriv(const double press, const double* surfvol, int pvtTableIdx,
double& Bval, double& dBdpval) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
// To handle no-gas case.
Bval = 1.0;
dBdpval = 0.0;
return;
}
Bval = evalB(press, surfvol, pvtTableIdx);
dBdpval = -Bval*Bval*miscible_gas(press, surfvol, pvtTableIdx, 1, true);
}
double PvtLiveGas::evalR(const double press, const double* surfvol, int pvtTableIdx) const
{
if (surfvol[phase_pos_[Liquid]] == 0.0) {
// To handle no-gas case.
return 0.0;
}
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][4], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
return satR;
} else {
// Undersaturated case
return maxR;
}
}
void PvtLiveGas::evalRDeriv(const double press, const double* surfvol, int pvtTableIdx,
double& Rval, double& dRdpval) const
{
if (surfvol[phase_pos_[Liquid]] == 0.0) {
// To handle no-gas case.
Rval = 0.0;
dRdpval = 0.0;
return;
}
double satR = linearInterpolation(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][4], press);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (satR < maxR ) {
// Saturated case
Rval = satR;
dRdpval = linearInterpolationDerivative(saturated_gas_table_[pvtTableIdx][0],
saturated_gas_table_[pvtTableIdx][4],
press);
} else {
// Undersaturated case
Rval = maxR;
dRdpval = 0.0;
}
}
double PvtLiveGas::miscible_gas(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv) const
{
const std::vector<std::vector<double> > &saturatedGasTable =
saturated_gas_table_[pvtTableIdx];
const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
undersat_gas_tables_[pvtTableIdx];
int section;
double Rval = linearInterpolation(saturatedGasTable[0],
saturatedGasTable[4], press,
section);
double maxR = surfvol[phase_pos_[Liquid]]/surfvol[phase_pos_[Vapour]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolationDerivative(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = section;
if (undersatGasTables[is][0].size() < 2) {
double val = (saturatedGasTable[item][is+1]
- saturatedGasTable[item][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
maxR);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
maxR);
double val = (val2 - val1)/
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
return val;
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolation(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = section;
// Extrapolate from first table section
if (is == 0 && press < saturatedGasTable[0][0]) {
return linearInterpolation(undersatGasTables[0][0],
undersatGasTables[0][item],
maxR);
}
// Extrapolate from last table section
int ltp = saturatedGasTable[0].size() - 1;
if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
return linearInterpolation(undersatGasTables[ltp][0],
undersatGasTables[ltp][item],
maxR);
}
// Interpolate between table sections
double w = (press - saturatedGasTable[0][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
if (undersatGasTables[is][0].size() < 2) {
double val = saturatedGasTable[item][is] +
w*(saturatedGasTable[item][is+1] -
saturatedGasTable[item][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
maxR);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
maxR);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double PvtLiveGas::miscible_gas(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv) const
{
const std::vector<std::vector<double> > &saturatedGasTable =
saturated_gas_table_[pvtTableIdx];
const std::vector<std::vector<std::vector<double> > > &undersatGasTables =
undersat_gas_tables_[pvtTableIdx];
const bool isSat = cond.hasFreeOil();
// Derivative w.r.t p
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturatedGasTable[0], press);
if (undersatGasTables[is][0].size() < 2) {
double val = (saturatedGasTable[item][is+1]
- saturatedGasTable[item][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
r);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
r);
double val = (val2 - val1)/
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
return val;
}
} else if (deriv == 2){
if (isSat) {
return 0;
} else {
int is = tableIndex(saturatedGasTable[0], press);
double w = (press - saturatedGasTable[0][is]) /
(saturatedGasTable[0][is+1] - saturatedGasTable[0][is]);
assert(undersatGasTables[is][0].size() >= 2);
assert(undersatGasTables[is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersatGasTables[is][0],
undersatGasTables[is][item],
r);
double val2 =
linearInterpolationDerivative(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
r);
double val = val1 + w * (val2 - val1);
return val;
}
} else {
if (isSat) { // Saturated case
return linearInterpolation(saturatedGasTable[0],
saturatedGasTable[item],
press);
} else { // Undersaturated case
int is = tableIndex(saturatedGasTable[0], press);
// Extrapolate from last table section
//int ltp = saturatedGasTable[0].size() - 1;
//if (is+1 == ltp && press > saturatedGasTable[0][ltp]) {
// return linearInterpolation(undersatGasTables[ltp][0],
// undersatGasTables[ltp][item],
// r);
//}
// Interpolate between table sections
double w = (press - saturatedGasTable[0][is]) /
(saturatedGasTable[0][is+1] -
saturatedGasTable[0][is]);
if (undersatGasTables[is][0].size() < 2) {
double val = saturatedGasTable[item][is] +
w*(saturatedGasTable[item][is+1] -
saturatedGasTable[item][is]);
return val;
}
double val1 =
linearInterpolation(undersatGasTables[is][0],
undersatGasTables[is][item],
r);
double val2 =
linearInterpolation(undersatGasTables[is+1][0],
undersatGasTables[is+1][item],
r);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

182
opm/core/props/pvt/PvtLiveGas.hpp
View File

@@ -1,182 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PVTLIVEGAS_HEADER_INCLUDED
#define OPM_PVTLIVEGAS_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm
{
class PvtgTable;
/// Class for miscible wet gas (with vaporized oil in vapour phase).
/// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class PvtLiveGas : public PvtInterface
{
public:
PvtLiveGas(const std::vector<Opm::PvtgTable>& pvtgTables);
virtual ~PvtLiveGas();
/// Viscosity as a function of p, T and z.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const;
/// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p, T and z.
virtual void B(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its p and r derivatives as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtRegionIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p, T.
virtual void rsSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtRegionIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtRegionIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
protected:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
double evalB(double press, const double* surfvol, int pvtTableIdx) const;
void evalBDeriv(double press, const double* surfvol, int pvtTableIdx, double& B, double& dBdp) const;
double evalR(double press, const double* surfvol, int pvtTableIdx) const;
void evalRDeriv(double press, const double* surfvol, int pvtTableIdx, double& R, double& dRdp) const;
// item: 1=>B 2=>mu;
double miscible_gas(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv = false) const;
double miscible_gas(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv = 0) const;
// PVT properties of wet gas (with vaporised oil). We need to
// store one table per PVT region.
std::vector< std::vector<std::vector<double> > > saturated_gas_table_;
std::vector< std::vector<std::vector<std::vector<double> > > > undersat_gas_tables_;
};
}
#endif // OPM_PVTLIVEGAS_HEADER_INCLUDED

653
opm/core/props/pvt/PvtLiveOil.cpp
View File

@@ -1,653 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/props/pvt/PvtLiveOil.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/PvtoTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
#include <algorithm>
namespace Opm
{
using Opm::linearInterpolation;
using Opm::linearInterpolationDerivative;
using Opm::tableIndex;
//------------------------------------------------------------------------
// Member functions
//-------------------------------------------------------------------------
PvtLiveOil::PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables)
{
int numTables = pvtoTables.size();
saturated_oil_table_.resize(numTables);
undersat_oil_tables_.resize(numTables);
for (int pvtTableIdx = 0; pvtTableIdx < numTables; ++pvtTableIdx) {
const Opm::PvtoTable& pvtoTable = pvtoTables[pvtTableIdx];
const auto& saturatedPvto = pvtoTable.getSaturatedTable();
// OIL, PVTO
// saturated_oil_table_[pvtTableIdx].resize(4);
// adding a extra colummn to the PVTO to store 1/(B*mu)
saturated_oil_table_[pvtTableIdx].resize(5);
const int sz = saturatedPvto.numRows();
// for (int k=0; k<4; ++k) {
for (int k=0; k<5; ++k) {
saturated_oil_table_[pvtTableIdx][k].resize(sz);
}
for (int i=0; i<sz; ++i) {
saturated_oil_table_[pvtTableIdx][0][i] = saturatedPvto.get("P" , i); // P
saturated_oil_table_[pvtTableIdx][1][i] = 1.0/saturatedPvto.get("BO" , i); // BO
saturated_oil_table_[pvtTableIdx][2][i] = saturatedPvto.get("MU" , i); // MU
saturated_oil_table_[pvtTableIdx][3][i] = 1.0 / (saturatedPvto.get("BO" , i)
* saturatedPvto.get("MU" , i));; // 1/(Bo*mu_o)
saturated_oil_table_[pvtTableIdx][4][i] = saturatedPvto.get("RS" , i); // RS
}
undersat_oil_tables_[pvtTableIdx].resize(sz);
for (int i=0; i<sz; ++i) {
const auto& undersaturatedPvto = pvtoTable.getUnderSaturatedTable( i );
// undersat_oil_tables_[pvtTableIdx][i].resize(3);
// adding a extra colummn to the PVTO to store 1/(B*mu)
undersat_oil_tables_[pvtTableIdx][i].resize(4);
int tsize = undersaturatedPvto.numRows();
undersat_oil_tables_[pvtTableIdx][i][0].resize(tsize);
undersat_oil_tables_[pvtTableIdx][i][1].resize(tsize);
undersat_oil_tables_[pvtTableIdx][i][2].resize(tsize);
undersat_oil_tables_[pvtTableIdx][i][3].resize(tsize);
for (int j=0; j<tsize; ++j) {
undersat_oil_tables_[pvtTableIdx][i][0][j] = undersaturatedPvto.get("P" , j); // P
undersat_oil_tables_[pvtTableIdx][i][1][j] = 1.0/undersaturatedPvto.get("BO" , j); // BO
undersat_oil_tables_[pvtTableIdx][i][2][j] = undersaturatedPvto.get("MU" , j); // MU
undersat_oil_tables_[pvtTableIdx][i][3][j] = 1.0 / (undersaturatedPvto.get("BO", j)*
undersaturatedPvto.get("MU", j)); // 1/(Bo*mu_o)
}
}
// Complete undersaturated tables by extrapolating from existing data
// as is done in Eclipse and Mrst
// TODO: check if the following formulations applying to 1/(Bo*mu_o)
int iNext = -1;
for (int i=0; i<sz; ++i) {
// Skip records already containing undersaturated data
if (undersat_oil_tables_[pvtTableIdx][i][0].size() > 1) {
continue;
}
// Look ahead for next record containing undersaturated data
if (iNext < i) {
iNext = i+1;
while (iNext<sz && undersat_oil_tables_[pvtTableIdx][iNext][0].size() < 2) {
++iNext;
}
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
}
// Add undersaturated data to current record while maintaining compressibility and viscosibility
// TODO: How to add 1/(B*mu) in this way?
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
for (sz_t j=1; j<undersat_oil_tables_[pvtTableIdx][iNext][0].size(); ++j) {
double diffPressure = undersat_oil_tables_[pvtTableIdx][iNext][0][j]-undersat_oil_tables_[pvtTableIdx][iNext][0][j-1];
double pressure = undersat_oil_tables_[pvtTableIdx][i][0].back()+diffPressure;
undersat_oil_tables_[pvtTableIdx][i][0].push_back(pressure);
double compr = (1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]-1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1])
/ (0.5*(1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j]+1.0/undersat_oil_tables_[pvtTableIdx][iNext][1][j-1]));
double B_var = (1.0/undersat_oil_tables_[pvtTableIdx][i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
undersat_oil_tables_[pvtTableIdx][i][1].push_back(1.0/B_var);
double visc = (undersat_oil_tables_[pvtTableIdx][iNext][2][j]-undersat_oil_tables_[pvtTableIdx][iNext][2][j-1])
/ (0.5*(undersat_oil_tables_[pvtTableIdx][iNext][2][j]+undersat_oil_tables_[pvtTableIdx][iNext][2][j-1]));
double mu_var = (undersat_oil_tables_[pvtTableIdx][i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
undersat_oil_tables_[pvtTableIdx][i][2].push_back(mu_var);
// A try to expolate the 1/BMu with the expolated mu and B
double inverseBMu = 1.0 / (B_var*mu_var);
undersat_oil_tables_[pvtTableIdx][i][3].push_back(inverseBMu);
}
}
}
}
/// Destructor.
PvtLiveOil::~PvtLiveOil()
{
}
/// Viscosity as a function of p, T and z.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_mu) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
double inverseB = miscible_oil(p[i], z + num_phases_*i, tableIdx, 1, false);
double inverseBMu = miscible_oil(p[i], z + num_phases_*i, tableIdx, 3, false);
output_mu[i] = inverseB / inverseBMu;
}
}
/// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
double inverseBMu = miscible_oil(p[i], r[i], tableIdx, 3, 0);
double inverseB = miscible_oil(p[i], r[i], tableIdx, 1, 0);
output_mu[i] = inverseB / inverseBMu;
double dinverseBmudp = miscible_oil(p[i], r[i], tableIdx, 3, 1);
double dinverseBdp = miscible_oil(p[i], r[i], tableIdx, 1, 1);
output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
/ (inverseBMu * inverseBMu);
double dinverseBmudr = miscible_oil(p[i], r[i], tableIdx, 3, 2);
double dinverseBdr = miscible_oil(p[i], r[i], tableIdx, 1, 2);
output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
/ (inverseBMu * inverseBMu);
}
}
/// Viscosity and its p and r derivatives as a function of p, T and r.
void PvtLiveOil::mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
const PhasePresence& cnd = cond[i];
double inverseBMu = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 0);
double inverseB = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0);
output_mu[i] = inverseB / inverseBMu;
double dinverseBmudp = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 1);
double dinverseBdp = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1);
output_dmudp[i] = (inverseBMu * dinverseBdp - inverseB * dinverseBmudp)
/ (inverseBMu * inverseBMu);
double dinverseBmudr = miscible_oil(p[i], r[i], cnd, tableIdx, 3, 2);
double dinverseBdr = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2);
output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
/ (inverseBMu * inverseBMu);
}
}
/// Formation volume factor as a function of p, T and z.
void PvtLiveOil::B(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_B) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_B[i] = evalB(tableIdx, p[i], z + num_phases_*i);
}
}
/// Formation volume factor and p-derivative as functions of p, T and z.
void PvtLiveOil::dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* z,
double* output_B,
double* output_dBdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
evalBDeriv(tableIdx, p[i], z + num_phases_*i, output_B[i], output_dBdp[i]);
}
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = miscible_oil(p[i], r[i], tableIdx, 1, 0);
output_dbdp[i] = miscible_oil(p[i], r[i], tableIdx, 1, 1);
output_dbdr[i] = miscible_oil(p[i], r[i], tableIdx, 1, 2);
}
}
void PvtLiveOil::b(const int n,
const int* pvtTableIdx,
const double* p,
const double* /*T*/,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
const PhasePresence& cnd = cond[i];
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_b[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 0);
output_dbdp[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 1);
output_dbdr[i] = miscible_oil(p[i], r[i], cnd, tableIdx, 1, 2);
}
}
void PvtLiveOil::rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const
{
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_rsSat[i] = linearInterpolation(saturated_oil_table_[tableIdx][0],
saturated_oil_table_[tableIdx][4],p[i]);
output_drsSatdp[i] = linearInterpolationDerivative(saturated_oil_table_[tableIdx][0],
saturated_oil_table_[tableIdx][4],p[i]);
}
}
void PvtLiveOil::rvSat(const int n,
const int* /*pvtTableIdx*/,
const double* /*p*/,
double* output_rvSat,
double* output_drvSatdp) const
{
std::fill(output_rvSat, output_rvSat + n, 0.0);
std::fill(output_drvSatdp, output_drvSatdp + n, 0.0);
}
/// Solution factor as a function of p and z.
void PvtLiveOil::R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
output_R[i] = evalR(tableIdx, p[i], z + num_phases_*i);
}
}
/// Solution factor and p-derivative as functions of p and z.
void PvtLiveOil::dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const
{
// #pragma omp parallel for
for (int i = 0; i < n; ++i) {
int tableIdx = getTableIndex_(pvtTableIdx, i);
evalRDeriv(tableIdx, p[i], z + num_phases_*i, output_R[i], output_dRdp[i]);
}
}
// ---- Private methods ----
double PvtLiveOil::evalB(size_t pvtTableIdx,
double press, const double* surfvol) const
{
// if (surfvol[phase_pos_[Liquid]] == 0.0) return 1.0; // To handle no-oil case.
return 1.0/miscible_oil(press, surfvol, pvtTableIdx, 1, false);
}
void PvtLiveOil::evalBDeriv(size_t pvtTableIdx,
const double press, const double* surfvol,
double& Bval, double& dBdpval) const
{
Bval = evalB(pvtTableIdx, press, surfvol);
dBdpval = -Bval*Bval*miscible_oil(press, surfvol, pvtTableIdx, 1, true);
}
double PvtLiveOil::evalR(size_t pvtTableIdx,
double press, const double* surfvol) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
return 0.0;
}
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][4], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) { // Saturated case
return Rval;
} else {
return maxR; // Undersaturated case
}
}
void PvtLiveOil::evalRDeriv(size_t pvtTableIdx,
const double press, const double* surfvol,
double& Rval, double& dRdpval) const
{
if (surfvol[phase_pos_[Vapour]] == 0.0) {
Rval = 0.0;
dRdpval = 0.0;
return;
}
Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][4], press);
double maxR = surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (Rval < maxR ) {
// Saturated case
dRdpval = linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][4],
press);
} else {
// Undersaturated case
Rval = maxR;
dRdpval = 0.0;
}
}
// TODO: Check if this function need to be adapted to the 1/(B*mu) interpolation.
double PvtLiveOil::miscible_oil(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv) const
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][4],
press, section);
double maxR = (surfvol[phase_pos_[Liquid]] == 0.0) ? 0.0 : surfvol[phase_pos_[Vapour]]/surfvol[phase_pos_[Liquid]];
if (deriv) {
if (Rval < maxR ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR);
double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
} else {
if (Rval < maxR ) { // Saturated case
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], maxR);
double w = (maxR - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double PvtLiveOil::miscible_oil(const double press,
const double r,
const int pvtTableIdx,
const int item,
const int deriv) const
{
int section;
double Rval = linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][4],
press, section);
// derivative with respect to frist component (pressure)
if (deriv == 1) {
if (Rval <= r ) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
// derivative with respect to second component (r)
} else if (deriv == 2) {
if (Rval <= r ) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]);
return val;
}
} else {
if (Rval <= r ) { // Saturated case
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
double PvtLiveOil::miscible_oil(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv) const
{
const bool isSat = cond.hasFreeGas();
// derivative with respect to frist component (pressure)
if (deriv == 1) {
if (isSat) { // Saturated case
return linearInterpolationDerivative(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolationDerivative(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
// derivative with respect to second component (r)
} else if (deriv == 2) {
if (isSat) { // Saturated case
return 0;
} else { // Undersaturated case
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = (val2 - val1)/(saturated_oil_table_[pvtTableIdx][4][is+1]-saturated_oil_table_[pvtTableIdx][4][is]);
return val;
}
}
// no derivative
else {
if (isSat) { // Saturated case
return linearInterpolation(saturated_oil_table_[pvtTableIdx][0],
saturated_oil_table_[pvtTableIdx][item],
press);
} else { // Undersaturated case
// Interpolate between table sections
int is = tableIndex(saturated_oil_table_[pvtTableIdx][4], r);
double w = (r - saturated_oil_table_[pvtTableIdx][4][is]) /
(saturated_oil_table_[pvtTableIdx][4][is+1] - saturated_oil_table_[pvtTableIdx][4][is]);
assert(undersat_oil_tables_[pvtTableIdx][is][0].size() >= 2);
assert(undersat_oil_tables_[pvtTableIdx][is+1][0].size() >= 2);
double val1 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is][0],
undersat_oil_tables_[pvtTableIdx][is][item],
press);
double val2 =
linearInterpolation(undersat_oil_tables_[pvtTableIdx][is+1][0],
undersat_oil_tables_[pvtTableIdx][is+1][item],
press);
double val = val1 + w*(val2 - val1);
return val;
}
}
}
} // namespace Opm

190
opm/core/props/pvt/PvtLiveOil.hpp
View File

@@ -1,190 +0,0 @@
/*
Copyright 2010, 2011, 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PVTLIVEOIL_HEADER_INCLUDED
#define OPM_PVTLIVEOIL_HEADER_INCLUDED
#include <opm/core/props/pvt/PvtInterface.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm
{
class PvtoTable;
/// Class for miscible live oil (with dissolved gas in liquid phase).
/// The PVT properties can either be given as a function of pressure (p), temperature (T) and surface volume (z)
/// or pressure (p), temperature (T) and gas resolution factor (r).
/// For all the virtual methods, the following apply: p, r and z
/// are expected to be of size n, size n and n*num_phases, respectively.
/// Output arrays shall be of size n, and must be valid before
/// calling the method.
class PvtLiveOil : public PvtInterface
{
public:
PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables);
virtual ~PvtLiveOil();
/// Viscosity as a function of p, T and z.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_mu) const;
/// Viscosity and its p and r derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Viscosity as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void mu(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const;
/// Formation volume factor as a function of p, T and z.
virtual void B(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B) const;
/// Formation volume factor and p-derivative as functions of p, T and z.
virtual void dBdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* z,
double* output_B,
double* output_dBdp) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p, T and r.
/// The fluid is considered saturated if r >= rsSat(p).
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p, T and r.
/// State condition determined by 'cond'.
virtual void b(const int n,
const int* pvtTableIdx,
const double* p,
const double* T,
const double* r,
const PhasePresence* cond,
double* output_b,
double* output_dbdp,
double* output_dbdr) const;
/// Solution gas/oil ratio and its derivatives at saturated conditions as a function of p.
virtual void rsSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rsSat,
double* output_drsSatdp) const;
/// Vapor oil/gas ratio and its derivatives at saturated conditions as a function of p.
virtual void rvSat(const int n,
const int* pvtTableIdx,
const double* p,
double* output_rvSat,
double* output_drvSatdp) const;
/// Solution factor as a function of p and z.
virtual void R(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R) const;
/// Solution factor and p-derivative as functions of p and z.
virtual void dRdp(const int n,
const int* pvtTableIdx,
const double* p,
const double* z,
double* output_R,
double* output_dRdp) const;
private:
int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
{
if (!pvtTableIdx)
return 0;
return pvtTableIdx[cellIdx];
}
double evalB(size_t pvtTableIdx, double press, const double* surfvol) const;
void evalBDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& B, double& dBdp) const;
double evalR(size_t pvtTableIdx, double press, const double* surfvol) const;
void evalRDeriv(size_t pvtTableIdx, double press, const double* surfvol, double& R, double& dRdp) const;
// item: 1=>1/B 2=>mu;
double miscible_oil(const double press,
const double* surfvol,
const int pvtTableIdx,
const int item,
const bool deriv = false) const;
double miscible_oil(const double press,
const double r,
const int pvtTableIdx,
const int item,
const int deriv = 0) const;
double miscible_oil(const double press,
const double r,
const PhasePresence& cond,
const int pvtTableIdx,
const int item,
const int deriv = 0) const;
// PVT properties of live oil (with dissolved gas). We need to
// store one table per PVT region.
std::vector<std::vector<std::vector<double> > > saturated_oil_table_;
std::vector<std::vector<std::vector<std::vector<double> > > > undersat_oil_tables_;
};
}
#endif // OPM_PVTLIVEOIL_HEADER_INCLUDED

46
opm/core/utility/extractPvtTableIndex.cpp Normal file
View File

@@ -0,0 +1,46 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include "extractPvtTableIndex.hpp"
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/GridProperty.hpp>
#include <vector>
namespace Opm {
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int *compressedToCartesianCellIdx)
{
//Get the PVTNUM data
const std::vector<int>& pvtnumData = eclState->getIntGridProperty("PVTNUM")->getData();
// Convert this into an array of compressed cells
// Eclipse uses Fortran-style indices which start at 1
// instead of 0, we subtract 1.
pvtTableIdx.resize(numCompressed);
for (size_t cellIdx = 0; cellIdx < numCompressed; ++ cellIdx) {
size_t cartesianCellIdx = compressedToCartesianCellIdx ? compressedToCartesianCellIdx[cellIdx]:cellIdx;
assert(cartesianCellIdx < pvtnumData.size());
pvtTableIdx[cellIdx] = pvtnumData[cartesianCellIdx] - 1;
}
}
}

33
opm/core/utility/extractPvtTableIndex.hpp Normal file
View File

@@ -0,0 +1,33 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_EXTRACT_PVT_TABLE_INDEX_HPP
#define OPM_EXTRACT_PVT_TABLE_INDEX_HPP
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <vector>
namespace Opm {
void extractPvtTableIndex(std::vector<int> &pvtTableIdx,
Opm::EclipseStateConstPtr eclState,
size_t numCompressed,
const int *compressedToCartesianCellIdx);
}
#endif // OPM_EXTRACT_PVT_TABLE_INDEX_HPP

74
opm/core/utility/thresholdPressures.hpp
View File

@@ -128,8 +128,6 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
}
// calculate the initial fluid densities for the gravity correction.
std::vector<double> b(numCells);
std::vector<std::vector<double>> rho(numPhases);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
rho[phaseIdx].resize(numCells);
@@ -176,69 +174,55 @@ void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
}
}
// calculate the inverse formation volume factors for the active phases and each cell
// calculate the densities of the active phases for each cell
if (pu.phase_used[BlackoilPhases::Aqua]) {
std::vector<double> dummy(numCells*BlackoilPhases::MaxNumPhases);
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const PvtInterface& pvtw = props.pvt(wpos);
pvtw.b(numCells,
pvtRegion.data(),
phasePressure[wpos].data(),
initialState.temperature().data(),
initialState.gasoilratio().data(),
cond.data(),
b.data(),
dummy.data(),
dummy.data());
const auto& pvtw = props.waterPvt();
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
rho[wpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][wpos]*b[cellIdx];
int pvtRegionIdx = pvtRegion[cellIdx];
double T = initialState.temperature()[cellIdx];
double p = initialState.pressure()[cellIdx];
double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b;
}
}
if (pu.phase_used[BlackoilPhases::Liquid]) {
std::vector<double> dummy(numCells*BlackoilPhases::MaxNumPhases);
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const PvtInterface& pvto = props.pvt(opos);
pvto.b(numCells,
pvtRegion.data(),
phasePressure[opos].data(),
initialState.temperature().data(),
initialState.gasoilratio().data(),
cond.data(),
b.data(),
dummy.data(),
dummy.data());
const auto& pvto = props.oilPvt();
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
rho[opos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][opos]*b[cellIdx];
int pvtRegionIdx = pvtRegion[cellIdx];
double T = initialState.temperature()[cellIdx];
double p = initialState.pressure()[cellIdx];
double Rs = initialState.gasoilratio()[cellIdx];
double b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
if (pu.phase_used[BlackoilPhases::Vapour]) {
int gpos = pu.phase_pos[BlackoilPhases::Vapour];
rho[opos][cellIdx] +=
surfaceDensity[pvtRegion[cellIdx]][gpos]*initialState.gasoilratio()[cellIdx]*b[cellIdx];
rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b;
}
}
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
std::vector<double> dummy(numCells*BlackoilPhases::MaxNumPhases);
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const PvtInterface& pvtg = props.pvt(gpos);
pvtg.b(numCells,
pvtRegion.data(),
phasePressure[gpos].data(),
initialState.temperature().data(),
initialState.rv().data(),
cond.data(),
b.data(),
dummy.data(),
dummy.data());
const int gpos = pu.phase_pos[BlackoilPhases::Liquid];
const auto& pvtg = props.gasPvt();
for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
rho[gpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][gpos]*b[cellIdx];
int pvtRegionIdx = pvtRegion[cellIdx];
double T = initialState.temperature()[cellIdx];
double p = initialState.pressure()[cellIdx];
double Rv = initialState.rv()[cellIdx];
double b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
if (pu.phase_used[BlackoilPhases::Liquid]) {
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
rho[gpos][cellIdx] +=
surfaceDensity[pvtRegion[cellIdx]][opos]*initialState.rv()[cellIdx]*b[cellIdx];
int opos = pu.phase_pos[BlackoilPhases::Liquid];
rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b;
}
}
}

11
tests/norne_pvt.data
View File

@@ -5,11 +5,22 @@
-- Copyright (C) 2015 Statoil
RUNSPEC
TABDIMS
--ntsfun ntpvt nssfun nppvt ntfip nrpvt ntendp
2 2 33 60 16 60 /
DIMENS
1 1 1 /
GRID
PROPS
DENSITY
859.5 1033.0 0.854 / Justert 22/7
860.04 1033.0 0.853 / Justert 22/7
PVTG

388
tests/test_blackoilfluid.cpp
View File

@@ -1,388 +0,0 @@
#include <config.h>
#include <opm/core/props/pvt/PvtConstCompr.hpp>
#include <opm/core/props/pvt/PvtDead.hpp>
#include <opm/core/props/pvt/PvtDeadSpline.hpp>
#include <opm/core/props/pvt/PvtLiveOil.hpp>
#include <opm/core/props/pvt/PvtLiveGas.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/props/pvt/BlackoilPvtProperties.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseMode.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE BlackoilFluidTest
#define BOOST_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <memory>
#include <iostream>
#include <iterator>
#include <vector>
#include <string>
using namespace Opm;
using namespace std;
std::vector<std::shared_ptr<PvtInterface> > getProps(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclipseState, PhaseUsage phase_usage_){
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
int samples = 0;
std::shared_ptr<const TableManager> tables = eclipseState->getTableManager();
std::vector<std::shared_ptr<PvtInterface> > props_;
// Set the properties.
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (deck->hasKeyword("PVTW")) {
std::shared_ptr<PvtConstCompr> pvtw(new PvtConstCompr);
pvtw->initFromWater(deck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]] = pvtw;
} else {
// Eclipse 100 default.
props_[phase_usage_.phase_pos[Aqua]].reset(new PvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
}
}
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
const auto& pvdoTables = tables->getPvdoTables();
const auto& pvtoTables = tables->getPvtoTables();
if (!pvdoTables.empty()) {
if (samples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromOil(pvdoTables, samples);
props_[phase_usage_.phase_pos[Liquid]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromOil(pvdoTables);
props_[phase_usage_.phase_pos[Liquid]] = deadPvt;
}
} else if (!pvtoTables.empty()) {
props_[phase_usage_.phase_pos[Liquid]].reset(new PvtLiveOil(pvtoTables));
} else if (deck->hasKeyword("PVCDO")) {
std::shared_ptr<PvtConstCompr> pvcdo(new PvtConstCompr);
pvcdo->initFromOil(deck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]] = pvcdo;
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDO, PVCDO or PVTO\n");
}
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
const auto& pvdgTables = tables->getPvdgTables();
const auto& pvtgTables = tables->getPvtgTables();
if (!pvdgTables.empty()) {
if (samples > 0) {
std::shared_ptr<PvtDeadSpline> splinePvt(new PvtDeadSpline);
splinePvt->initFromGas(pvdgTables, samples);
props_[phase_usage_.phase_pos[Vapour]] = splinePvt;
} else {
std::shared_ptr<PvtDead> deadPvt(new PvtDead);
deadPvt->initFromGas(pvdgTables);
props_[phase_usage_.phase_pos[Vapour]] = deadPvt;
}
} else if (!pvtgTables.empty()) {
props_[phase_usage_.phase_pos[Vapour]].reset(new PvtLiveGas(pvtgTables));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
}
return props_;
}
void testmu(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> T, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
std::vector<double> mu(n);
std::vector<double> dmudp(n);
std::vector<double> dmudr(n);
std::vector<double> mu_new(n);
double dmudp_diff;
double dmudr_diff;
double dmudp_diff_u;
double dmudr_diff_u;
// test mu
for (int phase = 0; phase < np; ++phase) {
props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &T[0], &r[0], &condition[0], &mu_new[0], &dmudp[0], &dmudr[0]);
props_[phase]->mu(n, &pvtTableIdx[0], &p[0], &T[0], &z[0], &mu[0]);
dmudp_diff = (mu_new[1]-mu_new[0])/(p[1]-p[0]);
dmudr_diff = (mu_new[2]-mu_new[0])/(r[2]-r[0]);
dmudp_diff_u = (mu_new[4]-mu_new[3])/(p[4]-p[3]);
dmudr_diff_u = (mu_new[5]-mu_new[3])/(r[5]-r[3]);
for (int i = 0; i < n; ++i){
BOOST_CHECK_CLOSE(mu_new[i],mu[i],reltol);
}
// saturated case
BOOST_CHECK_CLOSE(dmudp_diff,dmudp[0],reltol);
BOOST_CHECK_CLOSE(dmudr_diff,dmudr[0],reltol);
// unsaturated case
BOOST_CHECK_CLOSE(dmudp_diff_u,dmudp[3],reltol);
BOOST_CHECK_CLOSE(dmudr_diff_u,dmudr[3],reltol);
}
}
void testb(const double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<double> T, std::vector<double> r,std::vector<double> z,
std::vector<std::shared_ptr<PvtInterface> > props_, std::vector<Opm::PhasePresence> condition)
{
// test b
std::vector<double> b(n);
std::vector<double> B(n);
std::vector<double> invB(n);
std::vector<double> dinvBdp(n);
std::vector<double> dBdp(n);
std::vector<double> dbdr(n);
std::vector<double> dbdp(n);
double dbdp_diff;
double dbdr_diff;
double dbdp_diff_u;
double dbdr_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase]->b(n, &pvtTableIdx[0], &p[0], &T[0], &r[0], &condition[0], &b[0], &dbdp[0], &dbdr[0]);
props_[phase]->dBdp(n, &pvtTableIdx[0], &p[0], &T[0], &z[0], &B[0], &dBdp[0]);
dbdp_diff = (b[1]-b[0])/(p[1]-p[0]);
dbdr_diff = (b[2]-b[0])/(r[2]-r[0]);
dbdp_diff_u = (b[4]-b[3])/(p[4]-p[3]);
dbdr_diff_u = (b[5]-b[3])/(r[5]-r[3]);
for (int i = 0; i < n; ++i){
invB[i] = 1/B[i];
dinvBdp[i] = -1/pow(B[i],2) * dBdp[i];
}
for (int i = 0; i < n; ++i){
BOOST_CHECK_CLOSE(invB[i],b[i] , reltol);
BOOST_CHECK_CLOSE(dinvBdp[i],dbdp[i] , reltol);
}
// saturated case
BOOST_CHECK_CLOSE(dbdp_diff,dbdp[0], reltol);
BOOST_CHECK_CLOSE(dbdr_diff,dbdr[0], reltol);
// unsaturated case
BOOST_CHECK_CLOSE(dbdp_diff_u,dbdp[3], reltol);
BOOST_CHECK_CLOSE(dbdr_diff_u,dbdr[3], reltol);
}
}
void testrsSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
// test bublepoint pressure
std::vector<double> rs(n);
std::vector<double> drsdp(n);
double drsdp_diff;
double drsdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rsSat(n, &pvtTableIdx[0], &p[0], &rs[0], &drsdp[0]);
drsdp_diff = (rs[1]-rs[0])/(p[1]-p[0]);
drsdp_diff_u = (rs[4]-rs[3])/(p[4]-p[3]);
// saturated case
BOOST_CHECK_CLOSE(drsdp_diff,drsdp[0], reltol);
// unsaturad case
BOOST_CHECK_CLOSE(drsdp_diff_u,drsdp[3], reltol);
}
}
void testrvSat(double reltol, int n, int np, const std::vector<int> &pvtTableIdx, std::vector<double> p, std::vector<std::shared_ptr<PvtInterface> > props_){
// test rv saturated
std::vector<double> rv(n);
std::vector<double> drvdp(n);
double drvdp_diff;
double drvdp_diff_u;
for (int phase = 0; phase < np; ++phase) {
props_[phase] ->rvSat(n, &pvtTableIdx[0], &p[0], &rv[0], &drvdp[0]);
drvdp_diff = (rv[1]-rv[0])/(p[1]-p[0]);
drvdp_diff_u = (rv[4]-rv[3])/(p[4]-p[3]);
// saturated case
BOOST_CHECK_CLOSE(drvdp_diff,drvdp[0], reltol);
// unsaturad case
BOOST_CHECK_CLOSE(drvdp_diff_u,drvdp[3], reltol);
}
}
BOOST_AUTO_TEST_CASE(test_liveoil)
{
// read eclipse deck
const std::string filename = "liveoil.DATA";
cout << "Reading deck: " << filename << endl;
Opm::ParserPtr parser(new Opm::Parser());
Opm::ParseMode parseMode;
Opm::DeckConstPtr deck(parser->parseFile(filename , parseMode));
Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck, parseMode));
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck, eclipseState, phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
// approximation of the derivatives.
const int n = 6;
const int np = phase_usage_.num_phases;
std::vector<int> pvtRegionIdx(n, 0);
// the relative tolerance in percentage for acceptable difference in values
const double reltolper = 1e-9;
// the relative tolerance in percentage for acceptable difference in values for viscosity
const double reltolpermu = 1e-1;
std::vector<double> p(n);
std::vector<double> T(n, 273.15 + 20);
std::vector<double> r(n);
std::vector<PhasePresence> condition(n);
std::vector<double> z(n * np);
// Used for forward difference calculations
const double h_p = 1e4;
const double h_r = 1;
// saturated
p[0] = 1e7;
p[1] = p[0] + h_p;
p[2] = p[0];
r[0] = 200;
r[1] = 200;
r[2] = 200 + h_r;
condition[0].setFreeGas();
condition[1].setFreeGas();
condition[2].setFreeGas();
// undersaturated
p[3] = p[0];
p[4] = p[1];
p[5] = p[2];
r[3] = 50;
r[4] = 50;
r[5] = 50 +h_r;
// Corresponing z factors, used to compare with the [p,z] interface
for (int i = 0; i < n; ++i) {
z[0+i*np] = 0; z[1+i*np] = 1;
z[2+i*np] = r[i];
}
testmu(reltolpermu, n, np, pvtRegionIdx, p, T, r,z, props_, condition);
testb(reltolper,n,np,pvtRegionIdx,p,T,r,z,props_,condition);
testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);
testrvSat(reltolper,n,np,pvtRegionIdx,p,props_);
}
BOOST_AUTO_TEST_CASE(test_wetgas)
{
// read eclipse deck
const std::string filename = "wetgas.DATA";
cout << "Reading deck: " << filename << endl;
Opm::ParserPtr parser(new Opm::Parser());
Opm::ParseMode parseMode;
Opm::DeckConstPtr deck(parser->parseFile(filename , parseMode));
Opm::EclipseStateConstPtr eclipseState(new EclipseState(deck, parseMode));
// setup pvt interface
PhaseUsage phase_usage_ = phaseUsageFromDeck(deck);
std::vector<std::shared_ptr<PvtInterface> > props_ = getProps(deck,eclipseState,phase_usage_);
// setup a test case. We will check 6 [p,r] pairs and compare them to both the [p,z] interface and a finite difference
// approximation of the derivatives.
const int n = 6;
const int np = phase_usage_.num_phases;
std::vector<int> pvtRegionIdx(n, 0);
// the relative tolerance in percentage for acceptable difference in values
const double reltolper = 1e-9;
// the relative tolerance in percentage for acceptable difference in values for viscosity
const double reltolpermu = 1e-1;
std::vector<double> p(n);
std::vector<double> T(n, 273.15+20);
std::vector<double> r(n);
std::vector<PhasePresence> condition(n);
std::vector<double> z(n * np);
// Used for forward difference calculations
const double h_p = 1e4;
const double h_r = 1e-7;
// saturated
p[0] = 1e7;
p[1] = p[0] + h_p;
p[2] = p[0];
r[0] = 5e-5;
r[1] = 5e-5;
r[2] = 5e-5 + h_r;
condition[0].setFreeOil();
condition[1].setFreeOil();
condition[2].setFreeOil();
// undersaturated
p[3] = p[0];
p[4] = p[1];
p[5] = p[2];
r[3] = 1e-5;
r[4] = 1e-5;
r[5] = 1e-5 +h_r;
// Corresponing z factors, used to compare with the [p,z] interface
for (int i = 0; i < n; ++i) {
z[0+i*np] = 0; z[1+i*np] = r[i];
z[2+i*np] = 1;
}
testmu(reltolpermu, n, np, pvtRegionIdx, p,T, r,z, props_, condition);
testb(reltolper,n,np,pvtRegionIdx,p,T,r,z,props_,condition);
testrsSat(reltolper,n,np,pvtRegionIdx,p,props_);
testrvSat(reltolper,n,np,pvtRegionIdx,p,props_);
}

109
tests/test_norne_pvt.cpp
View File

@@ -39,7 +39,7 @@
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/core/props/pvt/PvtLiveOil.hpp>
#include <opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp>
using namespace Opm;
@@ -57,19 +57,9 @@ using namespace Opm;
further semantic meaning.
*/
void check_vectors( const std::vector<double>& v1 , const std::vector<double>& v2) {
double tol = 1e-5;
for (decltype(v1.size()) i = 0; i < v1.size(); i++) {
BOOST_CHECK_CLOSE( v1[i] , v2[i] , tol );
}
}
void verify_norne_oil_pvt_region2(const TableManager& tableManager) {
auto pvtoTables = tableManager.getPvtoTables();
PvtLiveOil oilPvt(pvtoTables);
void verify_norne_oil_pvt_region1(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState) {
Opm::LiveOilPvt<double> oilPvt;
oilPvt.initFromDeck(deck, eclState);
std::vector<double> rs = {33, 33,
43, 43,
@@ -113,33 +103,37 @@ void verify_norne_oil_pvt_region2(const TableManager& tableManager) {
{
std::vector<int> tableIndex(P.size() , 0);
std::vector<double> mu(P.size());
std::vector<double> dmudp(P.size());
std::vector<double> dmudr(P.size());
std::vector<double> b(P.size());
std::vector<double> dbdp(P.size());
std::vector<double> dbdr(P.size());
// convert the pressures to SI units (bar to Pascal)
for (auto& value : P)
value *= Metric::Pressure;
// convert the gas dissolution factors to SI units
for (auto& value : rs)
value *= Metric::GasDissolutionFactor;
oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data());
oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data());
for (unsigned i = 0; i < P.size(); ++i) {
double mu;
double b;
double RsSat = oilPvt.saturatedGasDissolutionFactor(/*tableIndex=*/0, /*T=*/273.15, P[i]);
if (rs[i] >= RsSat) {
mu = oilPvt.saturatedViscosity(/*tableIndex=*/0, /*T=*/273.15, P[i]);
b = oilPvt.saturatedInverseFormationVolumeFactor(/*tableIndex=*/0, /*T=*/273.15, P[i]);
}
else {
mu = oilPvt.viscosity(/*tableIndex=*/0, /*T=*/273.15, P[i], rs[i]);
b = oilPvt.inverseFormationVolumeFactor(/*tableIndex=*/0, /*T=*/273.15, P[i], rs[i]);
}
check_vectors( mu , mu_expected );
check_vectors( b , b_expected );
BOOST_CHECK_CLOSE( mu , mu_expected[i], 1e-5 );
BOOST_CHECK_CLOSE( b , b_expected[i], 1e-5 );
}
}
}
void verify_norne_oil_pvt_region1(const TableManager& tableManager) {
auto pvtoTables = tableManager.getPvtoTables();
PvtLiveOil oilPvt(pvtoTables);
void verify_norne_oil_pvt_region2(Opm::DeckConstPtr deck, Opm::EclipseStateConstPtr eclState) {
Opm::LiveOilPvt<double> oilPvt;
oilPvt.initFromDeck(deck, eclState);
std::vector<double> rs = {21 , 21,
30 , 30,
@@ -254,50 +248,41 @@ void verify_norne_oil_pvt_region1(const TableManager& tableManager) {
0.57289091037, 0.56019050084,
0.55474601877, 0.55809201119, 0.54526832277};
{
std::vector<int> tableIndex(P.size() , 1);
std::vector<double> mu(P.size());
std::vector<double> dmudp(P.size());
std::vector<double> dmudr(P.size());
std::vector<double> b(P.size());
std::vector<double> dbdp(P.size());
std::vector<double> dbdr(P.size());
// convert the pressures to SI units (bar to Pascal)
for (auto& value : P)
value *= Metric::Pressure;
// convert the gas dissolution factors to SI units
for (auto& value : rs)
value *= Metric::GasDissolutionFactor;
oilPvt.mu( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , mu.data() , dmudp.data() , dmudr.data());
oilPvt.b( P.size() , tableIndex.data() , P.data() , NULL , rs.data() , b.data() , dbdp.data() , dbdr.data());
for (unsigned i = 0; i < P.size(); ++i) {
double mu;
double b;
double RsSat = oilPvt.saturatedGasDissolutionFactor(/*tableIndex=*/1, /*T=*/273.15, P[i]);
if (rs[i] >= RsSat) {
mu = oilPvt.saturatedViscosity(/*tableIndex=*/1, /*T=*/273.15, P[i]);
b = oilPvt.saturatedInverseFormationVolumeFactor(/*tableIndex=*/1, /*T=*/273.15, P[i]);
}
else {
mu = oilPvt.viscosity(/*tableIndex=*/1, /*T=*/273.15, P[i], rs[i]);
b = oilPvt.inverseFormationVolumeFactor(/*tableIndex=*/1, /*T=*/273.15, P[i], rs[i]);
}
check_vectors( mu , mu_expected );
check_vectors( b , b_expected );
BOOST_CHECK_CLOSE( mu , mu_expected[i], 1e-5 );
BOOST_CHECK_CLOSE( b , b_expected[i], 1e-5 );
}
}
TableManager loadTables( const std::string& deck_file) {
BOOST_AUTO_TEST_CASE( Test_Norne_PVT) {
Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }});
Opm::ParserPtr parser(new Parser());
std::shared_ptr<const Deck> deck;
deck = parser->parseFile("norne_pvt.data", parseMode);
deck = parser->parseFile(deck_file, parseMode);
return TableManager(*deck);
}
BOOST_AUTO_TEST_CASE( Test_Norne_PVT) {
TableManager tableManager = loadTables( "norne_pvt.data" );
verify_norne_oil_pvt_region1( tableManager );
verify_norne_oil_pvt_region2( tableManager );
Opm::EclipseStateConstPtr eclState(new EclipseState(deck, parseMode));
verify_norne_oil_pvt_region1( deck, eclState );
verify_norne_oil_pvt_region2( deck, eclState );
}