Merge from upstream.
This commit is contained in:
Bård Skaflestad
commit
403cf40dd0
@ -37,3 +37,5 @@ tests/test_cfs_tpfa
|
||||
tests/test_jacsys
|
||||
tests/test_readvector
|
||||
tests/test_sf2p
|
||||
tests/bo_fluid_p_and_z_deps
|
||||
tests/bo_fluid_pressuredeps
|
||||
|
||||
1
.hgtags
Normal file
1
.hgtags
Normal file
@ -0,0 +1 @@
|
||||
b3cc85c2b5ab3f5b283244624d06a1afab9f8d67 Version-InitialImport-svn935
|
||||
60
Makefile.am
60
Makefile.am
@ -3,12 +3,25 @@ ACLOCAL_AMFLAGS = -I m4
|
||||
# Recurse to build examples after libraries
|
||||
SUBDIRS = . tests examples
|
||||
|
||||
CPPFLAGS += $(BOOST_CPPFLAGS)
|
||||
LIBS += $(BOOST_LDFLAGS) \
|
||||
$(BOOST_FILESYSTEM_LIB) \
|
||||
$(BOOST_SYSTEM_LIB) \
|
||||
$(BOOST_DATE_TIME_LIB) \
|
||||
$(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
|
||||
$(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
|
||||
|
||||
# Declare libraries
|
||||
lib_LTLIBRARIES = libopmcore.la
|
||||
|
||||
libopmcore_la_SOURCES = \
|
||||
libopmcore_la_SOURCES = \
|
||||
opm/core/eclipse/EclipseGridInspector.cpp \
|
||||
opm/core/eclipse/EclipseGridParser.cpp \
|
||||
opm/core/fluid/blackoil/BlackoilPVT.cpp \
|
||||
opm/core/fluid/blackoil/MiscibilityDead.cpp \
|
||||
opm/core/fluid/blackoil/MiscibilityLiveGas.cpp \
|
||||
opm/core/fluid/blackoil/MiscibilityLiveOil.cpp \
|
||||
opm/core/fluid/blackoil/MiscibilityProps.cpp \
|
||||
opm/core/utility/MonotCubicInterpolator.cpp \
|
||||
opm/core/utility/parameters/Parameter.cpp \
|
||||
opm/core/utility/parameters/ParameterGroup.cpp \
|
||||
@ -18,21 +31,21 @@ opm/core/utility/parameters/tinyxml/tinystr.cpp \
|
||||
opm/core/utility/parameters/tinyxml/tinyxml.cpp \
|
||||
opm/core/utility/parameters/tinyxml/tinyxmlerror.cpp \
|
||||
opm/core/utility/parameters/tinyxml/tinyxmlparser.cpp \
|
||||
opm/core/utility/cart_grid.c \
|
||||
opm/core/utility/cpgpreprocess/geometry.c \
|
||||
opm/core/utility/cpgpreprocess/preprocess.c \
|
||||
opm/core/utility/cpgpreprocess/readvector.cpp \
|
||||
opm/core/utility/cpgpreprocess/cgridinterface.c \
|
||||
opm/core/utility/cpgpreprocess/sparsetable.c \
|
||||
opm/core/utility/cpgpreprocess/facetopology.c \
|
||||
opm/core/utility/cpgpreprocess/uniquepoints.c \
|
||||
opm/core/utility/cart_grid.c \
|
||||
opm/core/utility/cpgpreprocess/geometry.c \
|
||||
opm/core/utility/cpgpreprocess/preprocess.c \
|
||||
opm/core/utility/cpgpreprocess/readvector.cpp \
|
||||
opm/core/utility/cpgpreprocess/cgridinterface.c \
|
||||
opm/core/utility/cpgpreprocess/sparsetable.c \
|
||||
opm/core/utility/cpgpreprocess/facetopology.c \
|
||||
opm/core/utility/cpgpreprocess/uniquepoints.c \
|
||||
opm/core/utility/StopWatch.cpp \
|
||||
opm/core/linalg/sparse_sys.c \
|
||||
opm/core/pressure/cfsh.c \
|
||||
opm/core/pressure/flow_bc.c \
|
||||
opm/core/pressure/well.c \
|
||||
opm/core/pressure/fsh_common_impl.c \
|
||||
opm/core/pressure/fsh.c \
|
||||
opm/core/linalg/sparse_sys.c \
|
||||
opm/core/pressure/cfsh.c \
|
||||
opm/core/pressure/flow_bc.c \
|
||||
opm/core/pressure/well.c \
|
||||
opm/core/pressure/fsh_common_impl.c \
|
||||
opm/core/pressure/fsh.c \
|
||||
opm/core/pressure/tpfa/ifs_tpfa.c \
|
||||
opm/core/pressure/tpfa/compr_source.c \
|
||||
opm/core/pressure/tpfa/cfs_tpfa.c \
|
||||
@ -62,10 +75,20 @@ opm/core/eclipse/EclipseGridParser.hpp \
|
||||
opm/core/eclipse/EclipseGridParserHelpers.hpp \
|
||||
opm/core/eclipse/EclipseUnits.hpp \
|
||||
opm/core/eclipse/SpecialEclipseFields.hpp \
|
||||
opm/core/fluid/BlackoilFluid.hpp \
|
||||
opm/core/fluid/SimpleFluid2p.hpp \
|
||||
opm/core/fluid/blackoil/BlackoilDefs.hpp \
|
||||
opm/core/fluid/blackoil/BlackoilPVT.hpp \
|
||||
opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp \
|
||||
opm/core/fluid/blackoil/FluidStateBlackoil.hpp \
|
||||
opm/core/fluid/blackoil/MiscibilityDead.hpp \
|
||||
opm/core/fluid/blackoil/MiscibilityLiveGas.hpp \
|
||||
opm/core/fluid/blackoil/MiscibilityLiveOil.hpp \
|
||||
opm/core/fluid/blackoil/MiscibilityProps.hpp \
|
||||
opm/core/fluid/blackoil/MiscibilityWater.hpp \
|
||||
opm/core/utility/Average.hpp \
|
||||
opm/core/utility/ErrorMacros.hpp \
|
||||
opm/core/utility/Factory.hpp \
|
||||
opm/core/utility/linInt.hpp \
|
||||
opm/core/utility/MonotCubicInterpolator.hpp \
|
||||
opm/core/utility/parameters/Parameter.hpp \
|
||||
opm/core/utility/parameters/ParameterGroup.hpp \
|
||||
@ -81,7 +104,11 @@ opm/core/utility/RootFinders.hpp \
|
||||
opm/core/utility/SparseTable.hpp \
|
||||
opm/core/utility/SparseVector.hpp \
|
||||
opm/core/utility/StopWatch.hpp \
|
||||
opm/core/utility/UniformTableLinear.hpp \
|
||||
opm/core/utility/Units.hpp \
|
||||
opm/core/utility/buildUniformMonotoneTable.hpp \
|
||||
opm/core/utility/linInt.hpp \
|
||||
opm/core/utility/linearInterpolation.hpp \
|
||||
opm/core/utility/cpgpreprocess/readvector.hpp \
|
||||
opm/core/utility/cpgpreprocess/uniquepoints.h \
|
||||
opm/core/utility/cpgpreprocess/preprocess.h \
|
||||
@ -95,7 +122,6 @@ opm/core/utility/cart_grid.h \
|
||||
opm/core/well.h \
|
||||
opm/core/linalg/sparse_sys.h \
|
||||
opm/core/linalg/blas_lapack.h \
|
||||
opm/core/fluid/SimpleFluid2p.hpp \
|
||||
opm/core/grid.h \
|
||||
opm/core/GridAdapter.hpp \
|
||||
opm/core/pressure/fsh.h \
|
||||
|
||||
@ -4,13 +4,7 @@ $(BOOST_CPPFLAGS)
|
||||
|
||||
LDFLAGS = $(BOOST_LDFLAGS)
|
||||
|
||||
LDADD = \
|
||||
$(top_builddir)/libopmcore.la \
|
||||
$(BOOST_FILESYSTEM_LIB) \
|
||||
$(BOOST_SYSTEM_LIB) \
|
||||
$(BOOST_DATE_TIME_LIB) \
|
||||
$(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
|
||||
$(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
|
||||
LDADD = $(top_builddir)/libopmcore.la
|
||||
|
||||
noinst_PROGRAMS = \
|
||||
scaneclipsedeck
|
||||
|
||||
@ -83,7 +83,7 @@ namespace EclipseKeywords
|
||||
string("BULKMOD"), string("YOUNGMOD"), string("LAMEMOD"),
|
||||
string("SHEARMOD"), string("POISSONMOD"), string("PWAVEMOD"),
|
||||
string("MULTPV"), string("PRESSURE"), string("SGAS"),
|
||||
string("SWAT")
|
||||
string("SWAT"), string("SOIL")
|
||||
};
|
||||
const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]);
|
||||
|
||||
@ -95,7 +95,7 @@ namespace EclipseKeywords
|
||||
string("SGOF"), string("SWOF"), string("ROCK"),
|
||||
string("ROCKTAB"), string("WELSPECS"), string("COMPDAT"),
|
||||
string("WCONINJE"), string("WCONPROD"), string("WELTARG"),
|
||||
string("EQUIL"), string("PVCDO"),
|
||||
string("EQUIL"), string("PVCDO"), string("TSTEP"),
|
||||
// The following fields only have a dummy implementation
|
||||
// that allows us to ignore them.
|
||||
"SWFN",
|
||||
@ -111,7 +111,7 @@ namespace EclipseKeywords
|
||||
string("NSTACK"), string("SATNUM"),
|
||||
string("RPTRST"), string("ROIP"), string("RWIP"),
|
||||
string("RWSAT"), string("RPR"), string("WBHP"),
|
||||
string("WOIR"), string("TSTEP"), string("BOX"),
|
||||
string("WOIR"), string("BOX"),
|
||||
string("COORDSYS"), string("PBVD")
|
||||
};
|
||||
const int num_ignore_with_data = sizeof(ignore_with_data) / sizeof(ignore_with_data[0]);
|
||||
@ -249,11 +249,18 @@ void EclipseGridParser::readImpl(istream& is)
|
||||
readVectorData(is, floatmap[keyword]);
|
||||
break;
|
||||
case SpecialField: {
|
||||
std::tr1::shared_ptr<SpecialBase> sb_ptr = createSpecialField(is, keyword);
|
||||
if (sb_ptr) {
|
||||
specialmap[keyword] = sb_ptr;
|
||||
map<string, std::tr1::shared_ptr<SpecialBase> >::iterator pos =
|
||||
specialmap.find(keyword);
|
||||
if (pos == specialmap.end()) {
|
||||
std::tr1::shared_ptr<SpecialBase> sb_ptr =
|
||||
createSpecialField(is, keyword);
|
||||
if (sb_ptr) {
|
||||
specialmap[keyword] = sb_ptr;
|
||||
} else {
|
||||
THROW("Could not create field " << keyword);
|
||||
}
|
||||
} else {
|
||||
THROW("Could not create field " << keyword);
|
||||
pos->second->read(is);
|
||||
}
|
||||
break;
|
||||
}
|
||||
@ -341,7 +348,7 @@ void EclipseGridParser::convertToSI()
|
||||
} else if (key == "PORO" || key == "BULKMOD" || key == "YOUNGMOD" ||
|
||||
key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" ||
|
||||
key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" ||
|
||||
key == "SGAS" || key == "SWAT") {
|
||||
key == "SGAS" || key == "SWAT" || key == "SOIL") {
|
||||
unit = 1.0;
|
||||
} else if (key == "PRESSURE") {
|
||||
unit = units_.pressure;
|
||||
|
||||
@ -142,6 +142,7 @@ public:
|
||||
SPECIAL_FIELD(WELTARG);
|
||||
SPECIAL_FIELD(EQUIL);
|
||||
SPECIAL_FIELD(PVCDO);
|
||||
SPECIAL_FIELD(TSTEP);
|
||||
|
||||
// The following fields only have a dummy implementation
|
||||
// that allows us to ignore them.
|
||||
|
||||
@ -868,7 +868,7 @@ struct CompdatLine
|
||||
skin_factor_(0.0),
|
||||
D_factor_(-1e100),
|
||||
penetration_direct_("Z"),
|
||||
r0_(0.0)
|
||||
r0_(-1.0)
|
||||
{
|
||||
grid_ind_.resize(4);
|
||||
}
|
||||
@ -1018,13 +1018,18 @@ struct WCONINJE : public SpecialBase
|
||||
double_data[2] = wconinje_line.BHP_limit_;
|
||||
double_data[4] = wconinje_line.VFP_table_number_;
|
||||
double_data[5] = wconinje_line.concentration_;
|
||||
readDefaultedVectorData(is, double_data, 6);
|
||||
const int num_to_read = 6;
|
||||
int num_read = readDefaultedVectorData(is, double_data, num_to_read);
|
||||
wconinje_line.surface_flow_max_rate_ = double_data[0];
|
||||
wconinje_line.fluid_volume_max_rate_ = double_data[1];
|
||||
wconinje_line.BHP_limit_ = double_data[2];
|
||||
wconinje_line.THP_limit_ = double_data[3];
|
||||
wconinje_line.VFP_table_number_ = (int)double_data[4];
|
||||
wconinje_line.concentration_ = double_data[5];
|
||||
// HACK! Ignore any further items
|
||||
if (num_read == num_to_read) {
|
||||
ignoreSlashLine(is);
|
||||
}
|
||||
wconinje.push_back(wconinje_line);
|
||||
}
|
||||
}
|
||||
@ -1127,7 +1132,8 @@ struct WCONPROD : public SpecialBase
|
||||
double_data[6] = wconprod_line.THP_limit_;
|
||||
double_data[7] = wconprod_line.VFP_table_number_;
|
||||
double_data[8] = wconprod_line.artif_lift_quantity_;
|
||||
readDefaultedVectorData(is, double_data, 9);
|
||||
const int num_to_read = 9;
|
||||
int num_read = readDefaultedVectorData(is, double_data, num_to_read);
|
||||
wconprod_line.oil_max_rate_ = double_data[0];
|
||||
wconprod_line.water_max_rate_ = double_data[1];
|
||||
wconprod_line.gas_max_rate_ = double_data[2];
|
||||
@ -1138,6 +1144,11 @@ struct WCONPROD : public SpecialBase
|
||||
wconprod_line.VFP_table_number_ = (int)double_data[7];
|
||||
wconprod_line.artif_lift_quantity_ = double_data[8];
|
||||
wconprod.push_back(wconprod_line);
|
||||
// HACK! Ignore any further items
|
||||
if (num_read == num_to_read) {
|
||||
ignoreSlashLine(is);
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
@ -1408,6 +1419,38 @@ struct PVCDO : public SpecialBase
|
||||
}
|
||||
};
|
||||
|
||||
struct TSTEP : public SpecialBase
|
||||
{
|
||||
std::vector<double> tstep_;
|
||||
|
||||
virtual std::string name() const {return std::string("TSTEP");}
|
||||
|
||||
virtual void read(std::istream& is)
|
||||
{
|
||||
std::vector<double> tstep;
|
||||
readVectorData(is, tstep);
|
||||
if (!tstep.empty()) {
|
||||
tstep_.insert(tstep_.end(), tstep.begin(), tstep.end());
|
||||
}
|
||||
}
|
||||
|
||||
virtual void write(std::ostream& os) const
|
||||
{
|
||||
os << name() << '\n';
|
||||
copy(tstep_.begin(), tstep_.end(),
|
||||
std::ostream_iterator<double>(os, " "));
|
||||
os << '\n';
|
||||
}
|
||||
|
||||
virtual void convertToSI(const EclipseUnits& units)
|
||||
{
|
||||
int num_steps = tstep_.size();
|
||||
for (int i = 0; i < num_steps; ++i) {
|
||||
tstep_[i] *= units.time;
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
struct MultRec : public SpecialBase
|
||||
{
|
||||
virtual void read(std::istream& is)
|
||||
|
||||
609
opm/core/fluid/BlackoilFluid.hpp
Normal file
609
opm/core/fluid/BlackoilFluid.hpp
Normal file
@ -0,0 +1,609 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_BLACKOILFLUID_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILFLUID_HEADER_INCLUDED
|
||||
|
||||
|
||||
#include <opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp>
|
||||
#include <opm/core/fluid/blackoil/FluidStateBlackoil.hpp>
|
||||
#include <opm/core/fluid/blackoil/BlackoilPVT.hpp>
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <vector>
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
/// Forward declaration for typedef i BlackoilFluid.
|
||||
class BlackoilFluidData;
|
||||
|
||||
|
||||
/// Class responsible for computing all fluid properties from
|
||||
/// face pressures and composition.
|
||||
class BlackoilFluid : public BlackoilDefs
|
||||
{
|
||||
public:
|
||||
typedef FluidStateBlackoil FluidState;
|
||||
typedef BlackoilFluidData FluidData;
|
||||
|
||||
void init(const Dune::EclipseGridParser& parser)
|
||||
{
|
||||
fmi_params_.init(parser);
|
||||
// FluidSystemBlackoil<>::init(parser);
|
||||
pvt_.init(parser);
|
||||
const std::vector<double>& dens = parser.getDENSITY().densities_[0];
|
||||
surface_densities_[Oil] = dens[0];
|
||||
surface_densities_[Water] = dens[1];
|
||||
surface_densities_[Gas] = dens[2];
|
||||
}
|
||||
|
||||
FluidState computeState(PhaseVec phase_pressure, CompVec z) const
|
||||
{
|
||||
FluidState state;
|
||||
state.temperature_ = 300;
|
||||
state.phase_pressure_ = phase_pressure;
|
||||
state.surface_volume_ = z;
|
||||
// FluidSystemBlackoil<>::computeEquilibrium(state); // Sets everything but relperm and mobility.
|
||||
computeEquilibrium(state);
|
||||
FluidMatrixInteractionBlackoil<double>::kr(state.relperm_, fmi_params_, state.saturation_, state.temperature_);
|
||||
FluidMatrixInteractionBlackoil<double>::dkr(state.drelperm_, fmi_params_, state.saturation_, state.temperature_);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
state.mobility_[phase] = state.relperm_[phase]/state.viscosity_[phase];
|
||||
for (int p2 = 0; p2 < numPhases; ++p2) {
|
||||
// Ignoring pressure variation in viscosity for this one.
|
||||
state.dmobility_[phase][p2] = state.drelperm_[phase][p2]/state.viscosity_[phase];
|
||||
}
|
||||
}
|
||||
return state;
|
||||
}
|
||||
|
||||
|
||||
const CompVec& surfaceDensities() const
|
||||
{
|
||||
return surface_densities_;
|
||||
}
|
||||
|
||||
/// \param[in] A state matrix in fortran ordering
|
||||
PhaseVec phaseDensities(const double* A) const
|
||||
{
|
||||
PhaseVec phase_dens(0.0);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
for (int comp = 0; comp < numComponents; ++comp) {
|
||||
phase_dens[phase] += A[numComponents*phase + comp]*surface_densities_[comp];
|
||||
}
|
||||
}
|
||||
return phase_dens;
|
||||
}
|
||||
|
||||
|
||||
// ---- New interface ----
|
||||
|
||||
/// Input: p, z
|
||||
/// Output: B, R
|
||||
template <class States>
|
||||
void computeBAndR(States& states) const
|
||||
{
|
||||
const std::vector<PhaseVec>& p = states.phase_pressure;
|
||||
const std::vector<CompVec>& z = states.surface_volume_density;
|
||||
std::vector<PhaseVec>& B = states.formation_volume_factor;
|
||||
std::vector<PhaseVec>& R = states.solution_factor;
|
||||
pvt_.B(p, z, B);
|
||||
pvt_.R(p, z, R);
|
||||
}
|
||||
|
||||
/// Input: p, z
|
||||
/// Output: B, R, mu
|
||||
template <class States>
|
||||
void computePvtNoDerivs(States& states) const
|
||||
{
|
||||
computeBAndR(states);
|
||||
const std::vector<PhaseVec>& p = states.phase_pressure;
|
||||
const std::vector<CompVec>& z = states.surface_volume_density;
|
||||
std::vector<PhaseVec>& mu = states.viscosity;
|
||||
pvt_.getViscosity(p, z, mu);
|
||||
}
|
||||
|
||||
/// Input: p, z
|
||||
/// Output: B, dB/dp, R, dR/dp, mu
|
||||
template <class States>
|
||||
void computePvt(States& states) const
|
||||
{
|
||||
const std::vector<PhaseVec>& p = states.phase_pressure;
|
||||
const std::vector<CompVec>& z = states.surface_volume_density;
|
||||
std::vector<PhaseVec>& B = states.formation_volume_factor;
|
||||
std::vector<PhaseVec>& dB = states.formation_volume_factor_deriv;
|
||||
std::vector<PhaseVec>& R = states.solution_factor;
|
||||
std::vector<PhaseVec>& dR = states.solution_factor_deriv;
|
||||
std::vector<PhaseVec>& mu = states.viscosity;
|
||||
pvt_.dBdp(p, z, B, dB);
|
||||
pvt_.dRdp(p, z, R, dR);
|
||||
pvt_.getViscosity(p, z, mu);
|
||||
}
|
||||
|
||||
/// Input: B, R
|
||||
/// Output: A
|
||||
template <class States>
|
||||
void computeStateMatrix(States& states) const
|
||||
{
|
||||
int num = states.formation_volume_factor.size();
|
||||
states.state_matrix.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
const PhaseVec& B = states.formation_volume_factor[i];
|
||||
const PhaseVec& R = states.solution_factor[i];
|
||||
PhaseToCompMatrix& At = states.state_matrix[i];
|
||||
// Set the A matrix (A = RB^{-1})
|
||||
// Using A transposed (At) since we really want Fortran ordering:
|
||||
// ultimately that is what the opmpressure C library expects.
|
||||
At = 0.0;
|
||||
At[Aqua][Water] = 1.0/B[Aqua];
|
||||
At[Vapour][Gas] = 1.0/B[Vapour];
|
||||
At[Liquid][Gas] = R[Liquid]/B[Liquid];
|
||||
At[Vapour][Oil] = R[Vapour]/B[Vapour];
|
||||
At[Liquid][Oil] = 1.0/B[Liquid];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/// Input: z, B, dB/dp, R, dR/dp
|
||||
/// Output: A, u, sum(u), s, c, cT, ex
|
||||
template <class States>
|
||||
void computePvtDepending(States& states) const
|
||||
{
|
||||
int num = states.formation_volume_factor.size();
|
||||
states.state_matrix.resize(num);
|
||||
states.phase_volume_density.resize(num);
|
||||
states.total_phase_volume_density.resize(num);
|
||||
states.saturation.resize(num);
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
states.phase_compressibility.resize(num);
|
||||
states.total_compressibility.resize(num);
|
||||
states.experimental_term.resize(num);
|
||||
#endif
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
const CompVec& z = states.surface_volume_density[i];
|
||||
const PhaseVec& B = states.formation_volume_factor[i];
|
||||
const PhaseVec& R = states.solution_factor[i];
|
||||
PhaseToCompMatrix& At = states.state_matrix[i];
|
||||
PhaseVec& u = states.phase_volume_density[i];
|
||||
double& tot_phase_vol_dens = states.total_phase_volume_density[i];
|
||||
PhaseVec& s = states.saturation[i];
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
const PhaseVec& dB = states.formation_volume_factor_deriv[i];
|
||||
const PhaseVec& dR = states.solution_factor_deriv[i];
|
||||
PhaseVec& cp = states.phase_compressibility[i];
|
||||
double& tot_comp = states.total_compressibility[i];
|
||||
double& exp_term = states.experimental_term[i];
|
||||
computeSingleEquilibrium(B, dB, R, dR, z,
|
||||
At, u, tot_phase_vol_dens,
|
||||
s, cp, tot_comp, exp_term);
|
||||
#else
|
||||
computeSingleEquilibrium(B, R, z,
|
||||
At, u, tot_phase_vol_dens,
|
||||
s);
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
/// Input: s, mu
|
||||
/// Output: kr, lambda
|
||||
template <class States>
|
||||
void computeMobilitiesNoDerivs(States& states) const
|
||||
{
|
||||
int num = states.saturation.size();
|
||||
states.relperm.resize(num);
|
||||
states.mobility.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
const CompVec& s = states.saturation[i];
|
||||
const PhaseVec& mu = states.viscosity[i];
|
||||
PhaseVec& kr = states.relperm[i];
|
||||
PhaseVec& lambda = states.mobility[i];
|
||||
FluidMatrixInteractionBlackoil<double>::kr(kr, fmi_params_, s, 300.0);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
lambda[phase] = kr[phase]/mu[phase];
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
/// Input: s, mu
|
||||
/// Output: kr, dkr/ds, lambda, dlambda/ds
|
||||
template <class States>
|
||||
void computeMobilities(States& states) const
|
||||
{
|
||||
int num = states.saturation.size();
|
||||
states.relperm.resize(num);
|
||||
states.relperm_deriv.resize(num);
|
||||
states.mobility.resize(num);
|
||||
states.mobility_deriv.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
const CompVec& s = states.saturation[i];
|
||||
const PhaseVec& mu = states.viscosity[i];
|
||||
PhaseVec& kr = states.relperm[i];
|
||||
PhaseJacobian& dkr = states.relperm_deriv[i];
|
||||
PhaseVec& lambda = states.mobility[i];
|
||||
PhaseJacobian& dlambda = states.mobility_deriv[i];
|
||||
FluidMatrixInteractionBlackoil<double>::kr(kr, fmi_params_, s, 300.0);
|
||||
FluidMatrixInteractionBlackoil<double>::dkr(dkr, fmi_params_, s, 300.0);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
lambda[phase] = kr[phase]/mu[phase];
|
||||
for (int p2 = 0; p2 < numPhases; ++p2) {
|
||||
// Ignoring pressure variation in viscosity for this one.
|
||||
dlambda[phase][p2] = dkr[phase][p2]/mu[phase];
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
BlackoilPVT pvt_;
|
||||
FluidMatrixInteractionBlackoilParams<double> fmi_params_;
|
||||
CompVec surface_densities_;
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Assuming the surface volumes and the pressures of all
|
||||
* phases are known, compute everything except relperm and
|
||||
* mobility.
|
||||
*/
|
||||
template <class FluidState>
|
||||
void computeEquilibrium(FluidState& fluid_state) const
|
||||
{
|
||||
// Get B and R factors.
|
||||
const PhaseVec& p = fluid_state.phase_pressure_;
|
||||
const CompVec& z = fluid_state.surface_volume_;
|
||||
PhaseVec& B = fluid_state.formation_volume_factor_;
|
||||
B[Aqua] = pvt_.B(p[Aqua], z, Aqua);
|
||||
B[Vapour] = pvt_.B(p[Vapour], z, Vapour);
|
||||
B[Liquid] = pvt_.B(p[Liquid], z, Liquid);
|
||||
PhaseVec& R = fluid_state.solution_factor_;
|
||||
R[Aqua] = 0.0;
|
||||
R[Vapour] = pvt_.R(p[Vapour], z, Vapour);
|
||||
R[Liquid] = pvt_.R(p[Liquid], z, Liquid);
|
||||
|
||||
// Convenience vars.
|
||||
PhaseToCompMatrix& At = fluid_state.phase_to_comp_;
|
||||
PhaseVec& u = fluid_state.phase_volume_density_;
|
||||
double& tot_phase_vol_dens = fluid_state.total_phase_volume_density_;
|
||||
PhaseVec& s = fluid_state.saturation_;
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
PhaseVec dB;
|
||||
dB[Aqua] = pvt_.dBdp(p[Aqua], z, Aqua);
|
||||
dB[Vapour] = pvt_.dBdp(p[Vapour], z, Vapour);
|
||||
dB[Liquid] = pvt_.dBdp(p[Liquid], z, Liquid);
|
||||
PhaseVec dR;
|
||||
dR[Aqua] = 0.0;
|
||||
dR[Vapour] = pvt_.dRdp(p[Vapour], z, Vapour);
|
||||
dR[Liquid] = pvt_.dRdp(p[Liquid], z, Liquid);
|
||||
PhaseVec& cp = fluid_state.phase_compressibility_;
|
||||
double& tot_comp = fluid_state.total_compressibility_;
|
||||
double& exp_term = fluid_state.experimental_term_;
|
||||
computeSingleEquilibrium(B, dB, R, dR, z,
|
||||
At, u, tot_phase_vol_dens,
|
||||
s, cp, tot_comp, exp_term);
|
||||
#else
|
||||
computeSingleEquilibrium(B, R, z,
|
||||
At, u, tot_phase_vol_dens,
|
||||
s);
|
||||
#endif
|
||||
|
||||
// Compute viscosities.
|
||||
PhaseVec& mu = fluid_state.viscosity_;
|
||||
mu[Aqua] = pvt_.getViscosity(p[Aqua], z, Aqua);
|
||||
mu[Vapour] = pvt_.getViscosity(p[Vapour], z, Vapour);
|
||||
mu[Liquid] = pvt_.getViscosity(p[Liquid], z, Liquid);
|
||||
}
|
||||
|
||||
|
||||
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
static void computeSingleEquilibrium(const PhaseVec& B,
|
||||
const PhaseVec& dB,
|
||||
const PhaseVec& R,
|
||||
const PhaseVec& dR,
|
||||
const CompVec& z,
|
||||
PhaseToCompMatrix& At,
|
||||
PhaseVec& u,
|
||||
double& tot_phase_vol_dens,
|
||||
PhaseVec& s,
|
||||
PhaseVec& cp,
|
||||
double& tot_comp,
|
||||
double& exp_term)
|
||||
#else
|
||||
static void computeSingleEquilibrium(const PhaseVec& B,
|
||||
const PhaseVec& R,
|
||||
const CompVec& z,
|
||||
PhaseToCompMatrix& At,
|
||||
PhaseVec& u,
|
||||
double& tot_phase_vol_dens,
|
||||
PhaseVec& s)
|
||||
#endif
|
||||
{
|
||||
// Set the A matrix (A = RB^{-1})
|
||||
// Using At since we really want Fortran ordering
|
||||
// (since ultimately that is what the opmpressure
|
||||
// C library expects).
|
||||
// PhaseToCompMatrix& At = fluid_state.phase_to_comp_[i];
|
||||
At = 0.0;
|
||||
At[Aqua][Water] = 1.0/B[Aqua];
|
||||
At[Vapour][Gas] = 1.0/B[Vapour];
|
||||
At[Liquid][Gas] = R[Liquid]/B[Liquid];
|
||||
At[Vapour][Oil] = R[Vapour]/B[Vapour];
|
||||
At[Liquid][Oil] = 1.0/B[Liquid];
|
||||
|
||||
// Update phase volumes. This is the same as multiplying with A^{-1}
|
||||
// PhaseVec& u = fluid_state.phase_volume_density_[i];
|
||||
double detR = 1.0 - R[Vapour]*R[Liquid];
|
||||
u[Aqua] = B[Aqua]*z[Water];
|
||||
u[Vapour] = B[Vapour]*(z[Gas] - R[Liquid]*z[Oil])/detR;
|
||||
u[Liquid] = B[Liquid]*(z[Oil] - R[Vapour]*z[Gas])/detR;
|
||||
tot_phase_vol_dens = u[Aqua] + u[Vapour] + u[Liquid];
|
||||
|
||||
// PhaseVec& s = fluid_state.saturation_[i];
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
s[phase] = u[phase]/tot_phase_vol_dens;
|
||||
}
|
||||
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
// Phase compressibilities.
|
||||
// PhaseVec& cp = fluid_state.phase_compressibility_[i];
|
||||
// Set the derivative of the A matrix (A = RB^{-1})
|
||||
PhaseToCompMatrix dAt(0.0);
|
||||
dAt[Aqua][Water] = -dB[Aqua]/(B[Aqua]*B[Aqua]);
|
||||
dAt[Vapour][Gas] = -dB[Vapour]/(B[Vapour]*B[Vapour]);
|
||||
dAt[Liquid][Oil] = -dB[Liquid]/(B[Liquid]*B[Liquid]); // Different order than above.
|
||||
dAt[Liquid][Gas] = dAt[Liquid][Oil]*R[Liquid] + dR[Liquid]/B[Liquid];
|
||||
dAt[Vapour][Oil] = dAt[Vapour][Gas]*R[Vapour] + dR[Vapour]/B[Vapour];
|
||||
|
||||
PhaseToCompMatrix Ait;
|
||||
Dune::FMatrixHelp::invertMatrix(At, Ait);
|
||||
|
||||
PhaseToCompMatrix Ct;
|
||||
Dune::FMatrixHelp::multMatrix(dAt, Ait, Ct);
|
||||
|
||||
cp[Aqua] = Ct[Aqua][Water];
|
||||
cp[Liquid] = Ct[Liquid][Oil] + Ct[Liquid][Gas];
|
||||
cp[Vapour] = Ct[Vapour][Gas] + Ct[Vapour][Oil];
|
||||
tot_comp = cp*s;
|
||||
|
||||
// Experimental term.
|
||||
PhaseVec tmp1, tmp2, tmp3;
|
||||
Ait.mtv(z, tmp1);
|
||||
dAt.mtv(tmp1, tmp2);
|
||||
Ait.mtv(tmp2, tmp3);
|
||||
exp_term = tmp3[Aqua] + tmp3[Liquid] + tmp3[Gas];
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
|
||||
struct FaceFluidData : public BlackoilDefs
|
||||
{
|
||||
// Canonical state variables.
|
||||
std::vector<CompVec> surface_volume_density; // z
|
||||
std::vector<PhaseVec> phase_pressure; // p
|
||||
|
||||
// Variables from PVT functions.
|
||||
std::vector<PhaseVec> formation_volume_factor; // B
|
||||
std::vector<PhaseVec> solution_factor; // R
|
||||
|
||||
// Variables computed from PVT data.
|
||||
// The A matrices are all in Fortran order (or, equivalently,
|
||||
// we store the transposes).
|
||||
std::vector<PhaseToCompMatrix> state_matrix; // A' = (RB^{-1})'
|
||||
|
||||
// Variables computed from saturation.
|
||||
std::vector<PhaseVec> mobility; // lambda
|
||||
std::vector<PhaseJacobian> mobility_deriv; // dlambda/ds
|
||||
|
||||
// Gravity and/or capillary pressure potential differences.
|
||||
std::vector<PhaseVec> gravity_potential; // (\rho g \delta z)-ish contribution per face
|
||||
};
|
||||
|
||||
|
||||
struct AllFluidData : public BlackoilDefs
|
||||
{
|
||||
// Per-cell data
|
||||
AllFluidStates cell_data;
|
||||
std::vector<double> voldiscr;
|
||||
std::vector<double> relvoldiscr;
|
||||
|
||||
// Per-face data.
|
||||
FaceFluidData face_data;
|
||||
|
||||
public:
|
||||
template <class Grid, class Rock>
|
||||
void computeNew(const Grid& grid,
|
||||
const Rock& rock,
|
||||
const BlackoilFluid& fluid,
|
||||
const typename Grid::Vector gravity,
|
||||
const std::vector<PhaseVec>& cell_pressure,
|
||||
const std::vector<PhaseVec>& face_pressure,
|
||||
const std::vector<CompVec>& cell_z,
|
||||
const CompVec& bdy_z,
|
||||
const double dt)
|
||||
{
|
||||
int num_cells = cell_z.size();
|
||||
ASSERT(num_cells == grid.numCells());
|
||||
const int np = numPhases;
|
||||
const int nc = numComponents;
|
||||
BOOST_STATIC_ASSERT(np == nc);
|
||||
BOOST_STATIC_ASSERT(np == 3);
|
||||
|
||||
// p, z -> B, dB, R, dR, mu, A, dA, u, sum(u), s, c, cT, ex, kr, dkr, lambda, dlambda.
|
||||
cell_data.phase_pressure = cell_pressure;
|
||||
cell_data.surface_volume_density = cell_z;
|
||||
fluid.computePvt(cell_data);
|
||||
fluid.computePvtDepending(cell_data);
|
||||
fluid.computeMobilities(cell_data);
|
||||
|
||||
// Compute volume discrepancies.
|
||||
voldiscr.resize(num_cells);
|
||||
relvoldiscr.resize(num_cells);
|
||||
#pragma omp parallel for
|
||||
for (int cell = 0; cell < num_cells; ++cell) {
|
||||
double pv = rock.porosity(cell)*grid.cellVolume(cell);
|
||||
voldiscr[cell] = (cell_data.total_phase_volume_density[cell] - 1.0)*pv/dt;
|
||||
relvoldiscr[cell] = std::fabs(cell_data.total_phase_volume_density[cell] - 1.0);
|
||||
}
|
||||
|
||||
|
||||
// Compute upwinded face properties, including z.
|
||||
computeUpwindProperties(grid, fluid, gravity,
|
||||
cell_pressure, face_pressure,
|
||||
cell_z, bdy_z);
|
||||
|
||||
// Compute state matrices for faces.
|
||||
// p, z -> B, R, A
|
||||
face_data.phase_pressure = face_pressure;
|
||||
fluid.computeBAndR(face_data);
|
||||
fluid.computeStateMatrix(face_data);
|
||||
}
|
||||
|
||||
|
||||
template <class Grid>
|
||||
void computeUpwindProperties(const Grid& grid,
|
||||
const BlackoilFluid& fluid,
|
||||
const typename Grid::Vector gravity,
|
||||
const std::vector<PhaseVec>& cell_pressure,
|
||||
const std::vector<PhaseVec>& face_pressure,
|
||||
const std::vector<CompVec>& cell_z,
|
||||
const CompVec& bdy_z)
|
||||
{
|
||||
int num_faces = face_pressure.size();
|
||||
ASSERT(num_faces == grid.numFaces());
|
||||
bool nonzero_gravity = gravity.two_norm() > 0.0;
|
||||
face_data.state_matrix.resize(num_faces);
|
||||
face_data.mobility.resize(num_faces);
|
||||
face_data.mobility_deriv.resize(num_faces);
|
||||
face_data.gravity_potential.resize(num_faces);
|
||||
face_data.surface_volume_density.resize(num_faces);
|
||||
#pragma omp parallel for
|
||||
for (int face = 0; face < num_faces; ++face) {
|
||||
// Obtain properties from both sides of the face.
|
||||
typedef typename Grid::Vector Vec;
|
||||
Vec fc = grid.faceCentroid(face);
|
||||
int c[2] = { grid.faceCell(face, 0), grid.faceCell(face, 1) };
|
||||
|
||||
// Get pressures and compute gravity contributions,
|
||||
// to decide upwind directions.
|
||||
PhaseVec phase_p[2];
|
||||
PhaseVec gravcontrib[2];
|
||||
for (int j = 0; j < 2; ++j) {
|
||||
if (c[j] >= 0) {
|
||||
// Pressures
|
||||
phase_p[j] = cell_pressure[c[j]];
|
||||
// Gravity contribution.
|
||||
if (nonzero_gravity) {
|
||||
Vec cdiff = fc;
|
||||
cdiff -= grid.cellCentroid(c[j]);
|
||||
gravcontrib[j] = fluid.phaseDensities(&cell_data.state_matrix[c[j]][0][0]);
|
||||
gravcontrib[j] *= (cdiff*gravity);
|
||||
} else {
|
||||
gravcontrib[j] = 0.0;
|
||||
}
|
||||
} else {
|
||||
// Pressures
|
||||
phase_p[j] = face_pressure[face];
|
||||
// Gravity contribution.
|
||||
gravcontrib[j] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// Gravity contribution:
|
||||
// gravcapf = rho_1*g*(z_12 - z_1) - rho_2*g*(z_12 - z_2)
|
||||
// where _1 and _2 refers to two neigbour cells, z is the
|
||||
// z coordinate of the centroid, and z_12 is the face centroid.
|
||||
// Also compute the potentials.
|
||||
PhaseVec pot[2];
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
face_data.gravity_potential[face][phase] = gravcontrib[0][phase] - gravcontrib[1][phase];
|
||||
pot[0][phase] = phase_p[0][phase] + face_data.gravity_potential[face][phase];
|
||||
pot[1][phase] = phase_p[1][phase];
|
||||
}
|
||||
|
||||
// Now we can easily find the upwind direction for every phase,
|
||||
// we can also tell which boundary faces are inflow bdys.
|
||||
|
||||
// Compute face_z, which is averaged from the cells, unless on outflow or noflow bdy.
|
||||
// Get mobilities and derivatives.
|
||||
CompVec face_z(0.0);
|
||||
double face_z_factor = 0.5;
|
||||
PhaseVec phase_mob[2];
|
||||
PhaseJacobian phasemob_deriv[2];
|
||||
for (int j = 0; j < 2; ++j) {
|
||||
if (c[j] >= 0) {
|
||||
face_z += cell_z[c[j]];
|
||||
phase_mob[j] = cell_data.mobility[c[j]];
|
||||
phasemob_deriv[j] = cell_data.mobility_deriv[c[j]];
|
||||
} else if (pot[j][Liquid] > pot[(j+1)%2][Liquid]) {
|
||||
// Inflow boundary.
|
||||
face_z += bdy_z;
|
||||
FluidStateBlackoil bdy_state = fluid.computeState(face_pressure[face], bdy_z);
|
||||
phase_mob[j] = bdy_state.mobility_;
|
||||
phasemob_deriv[j] = bdy_state.dmobility_;
|
||||
} else {
|
||||
// For outflow or noflow boundaries, only cell z is used.
|
||||
face_z_factor = 1.0;
|
||||
// Also, make sure the boundary data are not used for mobilities.
|
||||
pot[j] = -1e100;
|
||||
}
|
||||
}
|
||||
face_z *= face_z_factor;
|
||||
face_data.surface_volume_density[face] = face_z;
|
||||
|
||||
// Computing upwind mobilities and derivatives
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
if (pot[0][phase] == pot[1][phase]) {
|
||||
// Average.
|
||||
double aver = 0.5*(phase_mob[0][phase] + phase_mob[1][phase]);
|
||||
face_data.mobility[face][phase] = aver;
|
||||
for (int p2 = 0; p2 < numPhases; ++p2) {
|
||||
face_data.mobility_deriv[face][phase][p2] = phasemob_deriv[0][phase][p2]
|
||||
+ phasemob_deriv[1][phase][p2];
|
||||
}
|
||||
} else {
|
||||
// Upwind.
|
||||
int upwind = pot[0][phase] > pot[1][phase] ? 0 : 1;
|
||||
face_data.mobility[face][phase] = phase_mob[upwind][phase];
|
||||
for (int p2 = 0; p2 < numPhases; ++p2) {
|
||||
face_data.mobility_deriv[face][phase][p2] = phasemob_deriv[upwind][phase][p2];
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_BLACKOILFLUID_HEADER_INCLUDED
|
||||
101
opm/core/fluid/blackoil/BlackoilDefs.hpp
Normal file
101
opm/core/fluid/blackoil/BlackoilDefs.hpp
Normal file
@ -0,0 +1,101 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_BLACKOILDEFS_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILDEFS_HEADER_INCLUDED
|
||||
|
||||
|
||||
#include <tr1/array>
|
||||
#include <iostream>
|
||||
#include <boost/static_assert.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
class BlackoilDefs
|
||||
{
|
||||
public:
|
||||
enum { numComponents = 3 };
|
||||
enum { numPhases = 3 };
|
||||
|
||||
enum ComponentIndex { Water = 0, Oil = 1, Gas = 2 };
|
||||
enum PhaseIndex { Aqua = 0, Liquid = 1, Vapour = 2 };
|
||||
|
||||
// We need types with operator= and constructor taking scalars
|
||||
// for the small vectors and matrices, to save us from having to
|
||||
// rewrite a large amount of code.
|
||||
template <typename T, int N>
|
||||
class SmallVec
|
||||
{
|
||||
public:
|
||||
SmallVec()
|
||||
{}
|
||||
SmallVec(const T& elem)
|
||||
{ assign(elem); }
|
||||
SmallVec& operator=(const T& elem)
|
||||
{ assign(elem); return *this; }
|
||||
const T& operator[](int i) const
|
||||
{ return data_[i]; }
|
||||
T& operator[](int i)
|
||||
{ return data_[i]; }
|
||||
template <typename U>
|
||||
void assign(const U& elem)
|
||||
{
|
||||
for (int i = 0; i < N; ++i) {
|
||||
data_[i] = elem;
|
||||
}
|
||||
}
|
||||
private:
|
||||
T data_[N];
|
||||
};
|
||||
template <typename T, int Rows, int Cols>
|
||||
class SmallMat
|
||||
{
|
||||
public:
|
||||
SmallMat()
|
||||
{}
|
||||
SmallMat(const T& elem)
|
||||
{ data_.assign(elem); }
|
||||
SmallMat& operator=(const T& elem)
|
||||
{ data_.assign(elem); return *this; }
|
||||
typedef SmallVec<T, Cols> RowType;
|
||||
const RowType& operator[](int i) const
|
||||
{ return data_[i]; }
|
||||
RowType& operator[](int i)
|
||||
{ return data_[i]; }
|
||||
private:
|
||||
SmallVec<RowType, Rows> data_;
|
||||
};
|
||||
|
||||
typedef double Scalar;
|
||||
typedef SmallVec<Scalar, numComponents> CompVec;
|
||||
typedef SmallVec<Scalar, numPhases> PhaseVec;
|
||||
BOOST_STATIC_ASSERT(int(numComponents) == int(numPhases));
|
||||
typedef SmallMat<Scalar, numComponents, numPhases> PhaseToCompMatrix;
|
||||
typedef SmallMat<Scalar, numPhases, numPhases> PhaseJacobian;
|
||||
// Attempting to guard against alignment issues.
|
||||
BOOST_STATIC_ASSERT(sizeof(CompVec) == numComponents*sizeof(Scalar));
|
||||
BOOST_STATIC_ASSERT(sizeof(PhaseVec) == numPhases*sizeof(Scalar));
|
||||
BOOST_STATIC_ASSERT(sizeof(PhaseToCompMatrix) == numComponents*numPhases*sizeof(Scalar));
|
||||
BOOST_STATIC_ASSERT(sizeof(PhaseJacobian) == numPhases*numPhases*sizeof(Scalar));
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // OPM_BLACKOILDEFS_HEADER_INCLUDED
|
||||
202
opm/core/fluid/blackoil/BlackoilPVT.cpp
Normal file
202
opm/core/fluid/blackoil/BlackoilPVT.cpp
Normal file
@ -0,0 +1,202 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
|
||||
#include <opm/core/fluid/blackoil/BlackoilPVT.hpp>
|
||||
#include <opm/core/fluid/blackoil/MiscibilityDead.hpp>
|
||||
#include <opm/core/fluid/blackoil/MiscibilityLiveOil.hpp>
|
||||
#include <opm/core/fluid/blackoil/MiscibilityLiveGas.hpp>
|
||||
#include <opm/core/fluid/blackoil/MiscibilityWater.hpp>
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/utility/Units.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linInt.hpp>
|
||||
|
||||
using namespace Dune;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
void BlackoilPVT::init(const Dune::EclipseGridParser& parser)
|
||||
{
|
||||
typedef std::vector<std::vector<std::vector<double> > > table_t;
|
||||
region_number_ = 0;
|
||||
|
||||
// Surface densities. Accounting for different orders in eclipse and our code.
|
||||
if (parser.hasField("DENSITY")) {
|
||||
const int region_number = 0;
|
||||
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
|
||||
const std::vector<double>& d = parser.getDENSITY().densities_[region_number];
|
||||
densities_[Aqua] = d[ECL_water];
|
||||
densities_[Vapour] = d[ECL_gas];
|
||||
densities_[Liquid] = d[ECL_oil];
|
||||
} else {
|
||||
THROW("Input is missing DENSITY\n");
|
||||
}
|
||||
|
||||
// Water PVT
|
||||
if (parser.hasField("PVTW")) {
|
||||
water_props_.reset(new MiscibilityWater(parser.getPVTW().pvtw_));
|
||||
} else {
|
||||
water_props_.reset(new MiscibilityWater(0.5*Dune::prefix::centi*Dune::unit::Poise)); // Eclipse 100 default
|
||||
}
|
||||
|
||||
// Oil PVT
|
||||
if (parser.hasField("PVDO")) {
|
||||
oil_props_.reset(new MiscibilityDead(parser.getPVDO().pvdo_));
|
||||
} else if (parser.hasField("PVTO")) {
|
||||
oil_props_.reset(new MiscibilityLiveOil(parser.getPVTO().pvto_));
|
||||
} else if (parser.hasField("PVCDO")) {
|
||||
oil_props_.reset(new MiscibilityWater(parser.getPVCDO().pvcdo_));
|
||||
} else {
|
||||
THROW("Input is missing PVDO and PVTO\n");
|
||||
}
|
||||
|
||||
// Gas PVT
|
||||
if (parser.hasField("PVDG")) {
|
||||
gas_props_.reset(new MiscibilityDead(parser.getPVDG().pvdg_));
|
||||
} else if (parser.hasField("PVTG")) {
|
||||
gas_props_.reset(new MiscibilityLiveGas(parser.getPVTG().pvtg_));
|
||||
} else {
|
||||
THROW("Input is missing PVDG and PVTG\n");
|
||||
}
|
||||
}
|
||||
|
||||
BlackoilPVT::CompVec BlackoilPVT::surfaceDensities() const
|
||||
{
|
||||
return densities_;
|
||||
}
|
||||
|
||||
double BlackoilPVT::getViscosity(double press, const CompVec& surfvol, PhaseIndex phase) const
|
||||
{
|
||||
return propsForPhase(phase).getViscosity(region_number_, press, surfvol);
|
||||
}
|
||||
|
||||
double BlackoilPVT::B(double press, const CompVec& surfvol, PhaseIndex phase) const
|
||||
{
|
||||
return propsForPhase(phase).B(region_number_, press, surfvol);
|
||||
}
|
||||
|
||||
double BlackoilPVT::dBdp(double press, const CompVec& surfvol, PhaseIndex phase) const
|
||||
{
|
||||
return propsForPhase(phase).dBdp(region_number_, press, surfvol);
|
||||
}
|
||||
|
||||
double BlackoilPVT::R(double press, const CompVec& surfvol, PhaseIndex phase) const
|
||||
{
|
||||
return propsForPhase(phase).R(region_number_, press, surfvol);
|
||||
}
|
||||
|
||||
double BlackoilPVT::dRdp(double press, const CompVec& surfvol, PhaseIndex phase) const
|
||||
{
|
||||
return propsForPhase(phase).dRdp(region_number_, press, surfvol);
|
||||
}
|
||||
|
||||
const MiscibilityProps& BlackoilPVT::propsForPhase(PhaseIndex phase) const
|
||||
{
|
||||
switch (phase) {
|
||||
case Aqua:
|
||||
return *water_props_;
|
||||
case Liquid:
|
||||
return *oil_props_;
|
||||
case Vapour:
|
||||
return *gas_props_;
|
||||
default:
|
||||
THROW("Unknown phase accessed: " << phase);
|
||||
}
|
||||
}
|
||||
|
||||
void BlackoilPVT::getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
propsForPhase(PhaseIndex(phase)).getViscosity(pressures, surfvol, phase, data1_);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i][phase] = data1_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void BlackoilPVT::B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
propsForPhase(PhaseIndex(phase)).B(pressures, surfvol, phase, data1_);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i][phase] = data1_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void BlackoilPVT::dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output_B,
|
||||
std::vector<PhaseVec>& output_dBdp) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output_B.resize(num);
|
||||
output_dBdp.resize(num);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
propsForPhase(PhaseIndex(phase)).dBdp(pressures, surfvol, phase, data1_, data2_);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output_B[i][phase] = data1_[i];
|
||||
output_dBdp[i][phase] = data2_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void BlackoilPVT::R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
propsForPhase(PhaseIndex(phase)).R(pressures, surfvol, phase, data1_);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i][phase] = data1_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void BlackoilPVT::dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output_R,
|
||||
std::vector<PhaseVec>& output_dRdp) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output_R.resize(num);
|
||||
output_dRdp.resize(num);
|
||||
for (int phase = 0; phase < numPhases; ++phase) {
|
||||
propsForPhase(PhaseIndex(phase)).dRdp(pressures, surfvol, phase, data1_, data2_);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output_R[i][phase] = data1_[i];
|
||||
output_dRdp[i][phase] = data2_[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
88
opm/core/fluid/blackoil/BlackoilPVT.hpp
Normal file
88
opm/core/fluid/blackoil/BlackoilPVT.hpp
Normal file
@ -0,0 +1,88 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_BLACKOILPVT_HEADER_INCLUDED
|
||||
#define OPM_BLACKOILPVT_HEADER_INCLUDED
|
||||
|
||||
|
||||
#include <opm/core/fluid/blackoil/MiscibilityProps.hpp>
|
||||
#include <opm/core/fluid/blackoil/BlackoilDefs.hpp>
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <boost/scoped_ptr.hpp>
|
||||
#include <string>
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
class BlackoilPVT : public BlackoilDefs
|
||||
{
|
||||
public:
|
||||
void init(const Dune::EclipseGridParser& ep);
|
||||
|
||||
double getViscosity(double press,
|
||||
const CompVec& surfvol,
|
||||
PhaseIndex phase) const;
|
||||
CompVec surfaceDensities() const;
|
||||
double B (double press,
|
||||
const CompVec& surfvol,
|
||||
PhaseIndex phase) const;
|
||||
double dBdp(double press,
|
||||
const CompVec& surfvol,
|
||||
PhaseIndex phase) const;
|
||||
double R (double press,
|
||||
const CompVec& surfvol,
|
||||
PhaseIndex phase) const;
|
||||
double dRdp(double press,
|
||||
const CompVec& surfvol,
|
||||
PhaseIndex phase) const;
|
||||
|
||||
void getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output) const;
|
||||
void B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output) const;
|
||||
void dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output_B,
|
||||
std::vector<PhaseVec>& output_dBdp) const;
|
||||
void R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output) const;
|
||||
void dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
std::vector<PhaseVec>& output_R,
|
||||
std::vector<PhaseVec>& output_dRdp) const;
|
||||
|
||||
private:
|
||||
int region_number_;
|
||||
const MiscibilityProps& propsForPhase(PhaseIndex phase) const;
|
||||
|
||||
boost::scoped_ptr<MiscibilityProps> water_props_;
|
||||
boost::scoped_ptr<MiscibilityProps> oil_props_;
|
||||
boost::scoped_ptr<MiscibilityProps> gas_props_;
|
||||
CompVec densities_;
|
||||
mutable std::vector<double> data1_;
|
||||
mutable std::vector<double> data2_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // OPM_BLACKOILPVT_HEADER_INCLUDED
|
||||
208
opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp
Normal file
208
opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp
Normal file
@ -0,0 +1,208 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
|
||||
#define OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/utility/UniformTableLinear.hpp>
|
||||
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
|
||||
#include <opm/core/fluid/blackoil/BlackoilDefs.hpp>
|
||||
#include <iostream>
|
||||
#include <stdexcept>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
// Forward declaration needed by associated parameters class.
|
||||
template <class ScalarT, class ParamsT>
|
||||
class FluidMatrixInteractionBlackoil;
|
||||
|
||||
template <class ScalarT>
|
||||
class FluidMatrixInteractionBlackoilParams
|
||||
{
|
||||
public:
|
||||
typedef ScalarT Scalar;
|
||||
void init(const Dune::EclipseGridParser& ep)
|
||||
{
|
||||
// Extract input data.
|
||||
const Dune::SGOF::table_t& sgof_table = ep.getSGOF().sgof_;
|
||||
const Dune::SWOF::table_t& swof_table = ep.getSWOF().swof_;
|
||||
if (sgof_table.size() != 1 || swof_table.size() != 1) {
|
||||
std::cerr << "We must have exactly one SWOF and one SGOF table (at the moment).\n";
|
||||
throw std::logic_error("Not implemented");
|
||||
}
|
||||
const std::vector<double>& sw = swof_table[0][0];
|
||||
const std::vector<double>& krw = swof_table[0][1];
|
||||
const std::vector<double>& krow = swof_table[0][2];
|
||||
const std::vector<double>& pcow = swof_table[0][3];
|
||||
const std::vector<double>& sg = sgof_table[0][0];
|
||||
const std::vector<double>& krg = sgof_table[0][1];
|
||||
const std::vector<double>& krog = sgof_table[0][2];
|
||||
const std::vector<double>& pcog = sgof_table[0][3];
|
||||
|
||||
// Create tables for krw, krow, krg and krog.
|
||||
int samples = 200;
|
||||
buildUniformMonotoneTable(sw, krw, samples, krw_);
|
||||
buildUniformMonotoneTable(sw, krow, samples, krow_);
|
||||
buildUniformMonotoneTable(sg, krg, samples, krg_);
|
||||
buildUniformMonotoneTable(sg, krog, samples, krog_);
|
||||
krocw_ = krow[0]; // At connate water -> ecl. SWOF
|
||||
|
||||
// Create tables for pcow and pcog.
|
||||
buildUniformMonotoneTable(sw, pcow, samples, pcow_);
|
||||
buildUniformMonotoneTable(sg, pcog, samples, pcog_);
|
||||
}
|
||||
|
||||
private:
|
||||
template <class S, class P>
|
||||
friend class FluidMatrixInteractionBlackoil;
|
||||
|
||||
Dune::utils::UniformTableLinear<Scalar> krw_;
|
||||
Dune::utils::UniformTableLinear<Scalar> krow_;
|
||||
Dune::utils::UniformTableLinear<Scalar> pcow_;
|
||||
Dune::utils::UniformTableLinear<Scalar> krg_;
|
||||
Dune::utils::UniformTableLinear<Scalar> krog_;
|
||||
Dune::utils::UniformTableLinear<Scalar> pcog_;
|
||||
Scalar krocw_; // = krow_(s_wc)
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
* \ingroup material
|
||||
*
|
||||
* \brief Capillary pressures and relative permeabilities for a black oil system.
|
||||
*/
|
||||
template <class ScalarT, class ParamsT = FluidMatrixInteractionBlackoilParams<ScalarT> >
|
||||
class FluidMatrixInteractionBlackoil : public BlackoilDefs
|
||||
{
|
||||
public:
|
||||
typedef ParamsT Params;
|
||||
typedef typename Params::Scalar Scalar;
|
||||
|
||||
/*!
|
||||
* \brief The linear capillary pressure-saturation curve.
|
||||
*
|
||||
* This material law is linear:
|
||||
* \f[
|
||||
p_C = (1 - \overline{S}_w) (p_{C,max} - p_{C,entry}) + p_{C,entry}
|
||||
\f]
|
||||
*
|
||||
* \param Swe Effective saturation of of the wetting phase \f$\overline{S}_w\f$
|
||||
*/
|
||||
template <class pcContainerT, class SatContainerT>
|
||||
static void pC(pcContainerT &pc,
|
||||
const Params ¶ms,
|
||||
const SatContainerT &saturations,
|
||||
Scalar /*temperature*/)
|
||||
{
|
||||
Scalar sw = saturations[Aqua];
|
||||
Scalar sg = saturations[Vapour];
|
||||
pc[Liquid] = 0.0;
|
||||
pc[Aqua] = params.pcow_(sw);
|
||||
pc[Vapour] = params.pcog_(sg);
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief The saturation-capillary pressure curve.
|
||||
*
|
||||
* This is the inverse of the capillary pressure-saturation curve:
|
||||
* \f[
|
||||
S_w = 1 - \frac{p_C - p_{C,entry}}{p_{C,max} - p_{C,entry}}
|
||||
\f]
|
||||
*
|
||||
* \param pC Capillary pressure \f$\p_C\f$
|
||||
* \return The effective saturaion of the wetting phase \f$\overline{S}_w\f$
|
||||
*/
|
||||
template <class SatContainerT, class pcContainerT>
|
||||
static void S(SatContainerT &saturations,
|
||||
const Params ¶ms,
|
||||
const pcContainerT &pc,
|
||||
Scalar /*temperature*/)
|
||||
{
|
||||
std::cerr << "FluidMatrixInteractionBlackoil::S() is not implemented yet\n";
|
||||
throw std::logic_error("Not implemented");
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \brief The relative permeability of all phases.
|
||||
*/
|
||||
template <class krContainerT, class SatContainerT>
|
||||
static void kr(krContainerT &kr,
|
||||
const Params ¶ms,
|
||||
const SatContainerT &saturations,
|
||||
Scalar /*temperature*/)
|
||||
{
|
||||
// Stone-II relative permeability model.
|
||||
Scalar sw = saturations[Aqua];
|
||||
Scalar sg = saturations[Vapour];
|
||||
Scalar krw = params.krw_(sw);
|
||||
Scalar krg = params.krg_(sg);
|
||||
Scalar krow = params.krow_(sw);
|
||||
Scalar krog = params.krog_(sg);
|
||||
Scalar krocw = params.krocw_;
|
||||
kr[Aqua] = krw;
|
||||
kr[Vapour] = krg;
|
||||
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
|
||||
if (kr[Liquid] < 0.0) {
|
||||
kr[Liquid] = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
* \brief The saturation derivatives of relative permeability of all phases.
|
||||
*/
|
||||
template <class krContainerT, class SatContainerT>
|
||||
static void dkr(krContainerT &dkr,
|
||||
const Params ¶ms,
|
||||
const SatContainerT &saturations,
|
||||
Scalar /*temperature*/)
|
||||
{
|
||||
for (int p1 = 0; p1 < numPhases; ++p1) {
|
||||
for (int p2 = 0; p2 < numPhases; ++p2) {
|
||||
dkr[p1][p2] = 0.0;
|
||||
}
|
||||
}
|
||||
// Stone-II relative permeability model.
|
||||
Scalar sw = saturations[Aqua];
|
||||
Scalar sg = saturations[Vapour];
|
||||
Scalar krw = params.krw_(sw);
|
||||
Scalar dkrww = params.krw_.derivative(sw);
|
||||
Scalar krg = params.krg_(sg);
|
||||
Scalar dkrgg = params.krg_.derivative(sg);
|
||||
Scalar krow = params.krow_(sw);
|
||||
Scalar dkrow = params.krow_.derivative(sw);
|
||||
Scalar krog = params.krog_(sg);
|
||||
Scalar dkrog = params.krog_.derivative(sg);
|
||||
Scalar krocw = params.krocw_;
|
||||
dkr[Aqua][Aqua] = dkrww;
|
||||
dkr[Vapour][Vapour] = dkrgg;
|
||||
dkr[Liquid][Aqua] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
|
||||
dkr[Liquid][Vapour] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg);
|
||||
}
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
|
||||
|
||||
#endif // OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED
|
||||
95
opm/core/fluid/blackoil/FluidStateBlackoil.hpp
Normal file
95
opm/core/fluid/blackoil/FluidStateBlackoil.hpp
Normal file
@ -0,0 +1,95 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
|
||||
#define OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
|
||||
|
||||
|
||||
#include "BlackoilDefs.hpp"
|
||||
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
/*!
|
||||
* \brief Fluid state for a black oil model.
|
||||
*/
|
||||
struct FluidStateBlackoil : public BlackoilDefs
|
||||
{
|
||||
Scalar temperature_;
|
||||
CompVec surface_volume_;
|
||||
PhaseVec phase_pressure_;
|
||||
PhaseVec phase_volume_density_;
|
||||
Scalar total_phase_volume_density_;
|
||||
PhaseVec formation_volume_factor_;
|
||||
PhaseVec solution_factor_;
|
||||
PhaseToCompMatrix phase_to_comp_; // RB^{-1} in Fortran ordering
|
||||
PhaseVec saturation_;
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
PhaseVec phase_compressibility_;
|
||||
Scalar total_compressibility_;
|
||||
Scalar experimental_term_;
|
||||
#endif
|
||||
PhaseVec viscosity_;
|
||||
PhaseVec relperm_;
|
||||
PhaseJacobian drelperm_;
|
||||
PhaseVec mobility_;
|
||||
PhaseJacobian dmobility_;
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
* \brief Multiple fluid states for a black oil model.
|
||||
*/
|
||||
struct AllFluidStates : public BlackoilDefs
|
||||
{
|
||||
// Canonical state variables.
|
||||
std::vector<CompVec> surface_volume_density; // z
|
||||
std::vector<PhaseVec> phase_pressure; // p
|
||||
|
||||
// Variables from PVT functions.
|
||||
std::vector<PhaseVec> formation_volume_factor; // B
|
||||
std::vector<PhaseVec> solution_factor; // R
|
||||
std::vector<PhaseVec> viscosity; // mu
|
||||
|
||||
// Variables computed from PVT data.
|
||||
// The A matrices are all in Fortran order (or, equivalently,
|
||||
// we store the transposes).
|
||||
std::vector<PhaseToCompMatrix> state_matrix; // A' = (RB^{-1})'
|
||||
std::vector<PhaseVec> phase_volume_density; // u
|
||||
std::vector<Scalar> total_phase_volume_density; // sum(u)
|
||||
std::vector<PhaseVec> saturation; // s = u/sum(u)
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
std::vector<PhaseVec> formation_volume_factor_deriv; // dB/dp
|
||||
std::vector<PhaseVec> solution_factor_deriv; // dR/dp
|
||||
std::vector<PhaseVec> phase_compressibility; // c
|
||||
std::vector<Scalar> total_compressibility; // cT
|
||||
std::vector<Scalar> experimental_term; // ex = sum(Ai*dA*Ai*z)
|
||||
#endif
|
||||
|
||||
// Variables computed from saturation.
|
||||
std::vector<PhaseVec> relperm; // kr
|
||||
std::vector<PhaseJacobian> relperm_deriv; // dkr/ds
|
||||
std::vector<PhaseVec> mobility; // lambda
|
||||
std::vector<PhaseJacobian> mobility_deriv; // dlambda/ds
|
||||
};
|
||||
|
||||
} // end namespace Opm
|
||||
|
||||
|
||||
#endif // OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED
|
||||
11
opm/core/fluid/blackoil/Makefile.am
Normal file
11
opm/core/fluid/blackoil/Makefile.am
Normal file
@ -0,0 +1,11 @@
|
||||
noinst_LTLIBRARIES = libblackoilfluid_noinst.la
|
||||
|
||||
libblackoilfluid_noinst_la_SOURCES = \
|
||||
BlackoilPVT.cpp \
|
||||
MiscibilityDead.cpp \
|
||||
MiscibilityLiveGas.cpp \
|
||||
MiscibilityLiveOil.cpp \
|
||||
MiscibilityProps.cpp
|
||||
|
||||
|
||||
include $(top_srcdir)/am/global-rules
|
||||
179
opm/core/fluid/blackoil/MiscibilityDead.cpp
Normal file
179
opm/core/fluid/blackoil/MiscibilityDead.cpp
Normal file
@ -0,0 +1,179 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityDead.cpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:06:04 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <opm/core/fluid/blackoil/MiscibilityDead.hpp>
|
||||
#include <algorithm>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linInt.hpp>
|
||||
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
|
||||
#include <boost/lexical_cast.hpp>
|
||||
#include <string>
|
||||
#include <fstream>
|
||||
|
||||
using namespace std;
|
||||
using namespace Dune;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
|
||||
/// Constructor
|
||||
MiscibilityDead::MiscibilityDead(const table_t& pvd_table)
|
||||
{
|
||||
const int region_number = 0;
|
||||
if (pvd_table.size() != 1) {
|
||||
THROW("More than one PVT-region");
|
||||
}
|
||||
|
||||
// Copy data
|
||||
const int sz = pvd_table[region_number][0].size();
|
||||
std::vector<double> press(sz);
|
||||
std::vector<double> B_inv(sz);
|
||||
std::vector<double> visc(sz);
|
||||
for (int i = 0; i < sz; ++i) {
|
||||
press[i] = pvd_table[region_number][0][i];
|
||||
B_inv[i] = 1.0 / pvd_table[region_number][1][i];
|
||||
visc[i] = pvd_table[region_number][2][i];
|
||||
}
|
||||
int samples = 1025;
|
||||
buildUniformMonotoneTable(press, B_inv, samples, one_over_B_);
|
||||
buildUniformMonotoneTable(press, visc, samples, viscosity_);
|
||||
|
||||
// Dumping the created tables.
|
||||
// static int count = 0;
|
||||
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
|
||||
// os.precision(15);
|
||||
// os << "1/B\n\n" << one_over_B_
|
||||
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
|
||||
}
|
||||
|
||||
// Destructor
|
||||
MiscibilityDead::~MiscibilityDead()
|
||||
{
|
||||
}
|
||||
|
||||
double MiscibilityDead::getViscosity(int /*region*/, double press, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
return viscosity_(press);
|
||||
}
|
||||
|
||||
void MiscibilityDead::getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = viscosity_(pressures[i][phase]);
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityDead::B(int /*region*/, double press, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
// Interpolate 1/B
|
||||
return 1.0/one_over_B_(press);
|
||||
}
|
||||
|
||||
void MiscibilityDead::B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = 1.0/one_over_B_(pressures[i][phase]);
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityDead::dBdp(int region, double press, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
// Interpolate 1/B
|
||||
surfvol_t dummy_surfvol;
|
||||
double Bg = B(region, press, dummy_surfvol);
|
||||
return -Bg*Bg*one_over_B_.derivative(press);
|
||||
}
|
||||
|
||||
void MiscibilityDead::dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvols,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const
|
||||
{
|
||||
B(pressures, surfvols, phase, output_B);
|
||||
int num = pressures.size();
|
||||
output_dBdp.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
double Bg = output_B[i];
|
||||
output_dBdp[i] = -Bg*Bg*one_over_B_.derivative(pressures[i][phase]);
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityDead::R(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
void MiscibilityDead::R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.clear();
|
||||
output.resize(num, 0.0);
|
||||
}
|
||||
|
||||
double MiscibilityDead::dRdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
void MiscibilityDead::dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output_R.clear();
|
||||
output_R.resize(num, 0.0);
|
||||
output_dRdp.clear();
|
||||
output_dRdp.resize(num, 0.0);
|
||||
}
|
||||
|
||||
}
|
||||
89
opm/core/fluid/blackoil/MiscibilityDead.hpp
Normal file
89
opm/core/fluid/blackoil/MiscibilityDead.hpp
Normal file
@ -0,0 +1,89 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityDead.hpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:05:47 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef SINTEF_MISCIBILITYDEAD_HEADER
|
||||
#define SINTEF_MISCIBILITYDEAD_HEADER
|
||||
|
||||
/** Class for immiscible dead oil and dry gas.
|
||||
* Detailed description.
|
||||
*/
|
||||
|
||||
#include <opm/core/fluid/blackoil/MiscibilityProps.hpp>
|
||||
#include <opm/core/utility/UniformTableLinear.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
class MiscibilityDead : public MiscibilityProps
|
||||
{
|
||||
public:
|
||||
typedef std::vector<std::vector<std::vector<double> > > table_t;
|
||||
|
||||
MiscibilityDead(const table_t& pvd_table);
|
||||
virtual ~MiscibilityDead();
|
||||
|
||||
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double B(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double R(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
|
||||
|
||||
virtual void getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const;
|
||||
virtual void R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const;
|
||||
|
||||
private:
|
||||
// PVT properties of dry gas or dead oil
|
||||
Dune::utils::UniformTableLinear<double> one_over_B_;
|
||||
Dune::utils::UniformTableLinear<double> viscosity_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // SINTEF_MISCIBILITYDEAD_HEADER
|
||||
|
||||
286
opm/core/fluid/blackoil/MiscibilityLiveGas.cpp
Normal file
286
opm/core/fluid/blackoil/MiscibilityLiveGas.cpp
Normal file
@ -0,0 +1,286 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// 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.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
|
||||
#include <opm/core/fluid/blackoil/MiscibilityLiveGas.hpp>
|
||||
#include <algorithm>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linInt.hpp>
|
||||
|
||||
using namespace std;
|
||||
using namespace Dune;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
|
||||
/// Constructor
|
||||
MiscibilityLiveGas::MiscibilityLiveGas(const table_t& pvtg)
|
||||
{
|
||||
// GAS, PVTG
|
||||
const int region_number = 0;
|
||||
if (pvtg.size() != 1) {
|
||||
THROW("More than one PVD-region");
|
||||
}
|
||||
saturated_gas_table_.resize(4);
|
||||
const int sz = pvtg[region_number].size();
|
||||
for (int k=0; k<4; ++k) {
|
||||
saturated_gas_table_[k].resize(sz);
|
||||
}
|
||||
|
||||
for (int i=0; i<sz; ++i) {
|
||||
saturated_gas_table_[0][i] = pvtg[region_number][i][0]; // p
|
||||
saturated_gas_table_[1][i] = pvtg[region_number][i][2]; // Bg
|
||||
saturated_gas_table_[2][i] = pvtg[region_number][i][3]; // mu_g
|
||||
saturated_gas_table_[3][i] = pvtg[region_number][i][1]; // Rv
|
||||
}
|
||||
|
||||
undersat_gas_tables_.resize(sz);
|
||||
for (int i=0; i<sz; ++i) {
|
||||
undersat_gas_tables_[i].resize(3);
|
||||
int tsize = (pvtg[region_number][i].size() - 1)/3;
|
||||
undersat_gas_tables_[i][0].resize(tsize);
|
||||
undersat_gas_tables_[i][1].resize(tsize);
|
||||
undersat_gas_tables_[i][2].resize(tsize);
|
||||
for (int j=0, k=0; j<tsize; ++j) {
|
||||
undersat_gas_tables_[i][0][j] = pvtg[region_number][i][++k]; // Rv
|
||||
undersat_gas_tables_[i][1][j] = pvtg[region_number][i][++k]; // Bg
|
||||
undersat_gas_tables_[i][2][j] = pvtg[region_number][i][++k]; // mu_g
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Destructor
|
||||
MiscibilityLiveGas::~MiscibilityLiveGas()
|
||||
{
|
||||
}
|
||||
|
||||
double MiscibilityLiveGas::getViscosity(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
return miscible_gas(press, surfvol, 2, false);
|
||||
}
|
||||
|
||||
void MiscibilityLiveGas::getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = miscible_gas(pressures[i][phase], surfvol[i], 2, false);
|
||||
}
|
||||
}
|
||||
|
||||
// Vaporised oil-gas ratio.
|
||||
double MiscibilityLiveGas::R(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
if (surfvol[Liquid] == 0.0) {
|
||||
return 0.0;
|
||||
}
|
||||
double R = linearInterpolationExtrap(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3], press);
|
||||
double maxR = surfvol[Liquid]/surfvol[Vapour];
|
||||
if (R < maxR ) { // Saturated case
|
||||
return R;
|
||||
} else {
|
||||
return maxR; // Undersaturated case
|
||||
}
|
||||
}
|
||||
|
||||
void MiscibilityLiveGas::R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = R(0, pressures[i][phase], surfvol[i]);
|
||||
}
|
||||
}
|
||||
|
||||
// Vaporised oil-gas ratio derivative
|
||||
double MiscibilityLiveGas::dRdp(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
double R = linearInterpolationExtrap(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3], press);
|
||||
double maxR = surfvol[Liquid]/surfvol[Vapour];
|
||||
if (R < maxR ) { // Saturated case
|
||||
return linearInterpolDerivative(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3],
|
||||
press);
|
||||
} else {
|
||||
return 0.0; // Undersaturated case
|
||||
}
|
||||
}
|
||||
|
||||
void MiscibilityLiveGas::dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
R(pressures, surfvol, phase, output_R);
|
||||
int num = pressures.size();
|
||||
output_dRdp.resize(num);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output_dRdp[i] = dRdp(0, pressures[i][phase], surfvol[i]); // \TODO Speedup here by using already evaluated R.
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityLiveGas::B(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
if (surfvol[Vapour] == 0.0) return 1.0; // To handle no-gas case.
|
||||
return miscible_gas(press, surfvol, 1, false);
|
||||
}
|
||||
|
||||
void MiscibilityLiveGas::B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = B(0, pressures[i][phase], surfvol[i]);
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityLiveGas::dBdp(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
if (surfvol[Vapour] == 0.0) return 0.0; // To handle no-gas case.
|
||||
return miscible_gas(press, surfvol, 1, true);
|
||||
}
|
||||
|
||||
void MiscibilityLiveGas::dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
B(pressures, surfvol, phase, output_B);
|
||||
int num = pressures.size();
|
||||
output_dBdp.resize(num);
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output_dBdp[i] = dBdp(0, pressures[i][phase], surfvol[i]); // \TODO Speedup here by using already evaluated B.
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityLiveGas::miscible_gas(double press, const surfvol_t& surfvol, int item,
|
||||
bool deriv) const
|
||||
{
|
||||
int section;
|
||||
double R = linearInterpolationExtrap(saturated_gas_table_[0],
|
||||
saturated_gas_table_[3], press,
|
||||
section);
|
||||
double maxR = surfvol[Liquid]/surfvol[Vapour];
|
||||
if (deriv) {
|
||||
if (R < maxR ) { // Saturated case
|
||||
return linearInterpolDerivative(saturated_gas_table_[0],
|
||||
saturated_gas_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = section;
|
||||
if (undersat_gas_tables_[is][0].size() < 2) {
|
||||
double val = (saturated_gas_table_[item][is+1]
|
||||
- saturated_gas_table_[item][is]) /
|
||||
(saturated_gas_table_[0][is+1] -
|
||||
saturated_gas_table_[0][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolationExtrap(undersat_gas_tables_[is][0],
|
||||
undersat_gas_tables_[is][item],
|
||||
maxR);
|
||||
double val2 =
|
||||
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
|
||||
undersat_gas_tables_[is+1][item],
|
||||
maxR);
|
||||
double val = (val2 - val1)/
|
||||
(saturated_gas_table_[0][is+1] - saturated_gas_table_[0][is]);
|
||||
return val;
|
||||
}
|
||||
} else {
|
||||
if (R < maxR ) { // Saturated case
|
||||
return linearInterpolationExtrap(saturated_gas_table_[0],
|
||||
saturated_gas_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = section;
|
||||
// Extrapolate from first table section
|
||||
if (is == 0 && press < saturated_gas_table_[0][0]) {
|
||||
return linearInterpolationExtrap(undersat_gas_tables_[0][0],
|
||||
undersat_gas_tables_[0][item],
|
||||
maxR);
|
||||
}
|
||||
|
||||
// Extrapolate from last table section
|
||||
int ltp = saturated_gas_table_[0].size() - 1;
|
||||
if (is+1 == ltp && press > saturated_gas_table_[0][ltp]) {
|
||||
return linearInterpolationExtrap(undersat_gas_tables_[ltp][0],
|
||||
undersat_gas_tables_[ltp][item],
|
||||
maxR);
|
||||
}
|
||||
|
||||
// Interpolate between table sections
|
||||
double w = (press - saturated_gas_table_[0][is]) /
|
||||
(saturated_gas_table_[0][is+1] -
|
||||
saturated_gas_table_[0][is]);
|
||||
if (undersat_gas_tables_[is][0].size() < 2) {
|
||||
double val = saturated_gas_table_[item][is] +
|
||||
w*(saturated_gas_table_[item][is+1] -
|
||||
saturated_gas_table_[item][is]);
|
||||
return val;
|
||||
}
|
||||
double val1 =
|
||||
linearInterpolationExtrap(undersat_gas_tables_[is][0],
|
||||
undersat_gas_tables_[is][item],
|
||||
maxR);
|
||||
double val2 =
|
||||
linearInterpolationExtrap(undersat_gas_tables_[is+1][0],
|
||||
undersat_gas_tables_[is+1][item],
|
||||
maxR);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
92
opm/core/fluid/blackoil/MiscibilityLiveGas.hpp
Normal file
92
opm/core/fluid/blackoil/MiscibilityLiveGas.hpp
Normal file
@ -0,0 +1,92 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityLiveGas.hpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:21:26 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef SINTEF_MISCIBILITYLIVEGAS_HEADER
|
||||
#define SINTEF_MISCIBILITYLIVEGAS_HEADER
|
||||
|
||||
/** Class for miscible wet gas.
|
||||
* Detailed description.
|
||||
*/
|
||||
|
||||
#include "MiscibilityProps.hpp"
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
class MiscibilityLiveGas : public MiscibilityProps
|
||||
{
|
||||
public:
|
||||
typedef std::vector<std::vector<std::vector<double> > > table_t;
|
||||
|
||||
MiscibilityLiveGas(const table_t& pvto);
|
||||
virtual ~MiscibilityLiveGas();
|
||||
|
||||
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double R(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double B(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
|
||||
|
||||
virtual void getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const;
|
||||
virtual void R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const;
|
||||
|
||||
protected:
|
||||
// item: 1=B 2=mu;
|
||||
double miscible_gas(double press, const surfvol_t& surfvol, int item,
|
||||
bool deriv = false) const;
|
||||
// PVT properties of wet gas (with vaporised oil)
|
||||
std::vector<std::vector<double> > saturated_gas_table_;
|
||||
std::vector<std::vector<std::vector<double> > > undersat_gas_tables_;
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // SINTEF_MISCIBILITYLIVEGAS_HEADER
|
||||
|
||||
395
opm/core/fluid/blackoil/MiscibilityLiveOil.cpp
Normal file
395
opm/core/fluid/blackoil/MiscibilityLiveOil.cpp
Normal file
@ -0,0 +1,395 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibiltyLiveOil.cpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:08:25 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include <opm/core/fluid/blackoil/MiscibilityLiveOil.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linInt.hpp>
|
||||
#include <algorithm>
|
||||
|
||||
using namespace std;
|
||||
using namespace Dune;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
|
||||
/// Constructor
|
||||
MiscibilityLiveOil::MiscibilityLiveOil(const table_t& pvto)
|
||||
{
|
||||
// OIL, PVTO
|
||||
const int region_number = 0;
|
||||
if (pvto.size() != 1) {
|
||||
THROW("More than one PVD-region");
|
||||
}
|
||||
saturated_oil_table_.resize(4);
|
||||
const int sz = pvto[region_number].size();
|
||||
for (int k=0; k<4; ++k) {
|
||||
saturated_oil_table_[k].resize(sz);
|
||||
}
|
||||
for (int i=0; i<sz; ++i) {
|
||||
saturated_oil_table_[0][i] = pvto[region_number][i][1]; // p
|
||||
saturated_oil_table_[1][i] = 1.0/pvto[region_number][i][2]; // 1/Bo
|
||||
saturated_oil_table_[2][i] = pvto[region_number][i][3]; // mu_o
|
||||
saturated_oil_table_[3][i] = pvto[region_number][i][0]; // Rs
|
||||
}
|
||||
|
||||
undersat_oil_tables_.resize(sz);
|
||||
for (int i=0; i<sz; ++i) {
|
||||
undersat_oil_tables_[i].resize(3);
|
||||
int tsize = (pvto[region_number][i].size() - 1)/3;
|
||||
undersat_oil_tables_[i][0].resize(tsize);
|
||||
undersat_oil_tables_[i][1].resize(tsize);
|
||||
undersat_oil_tables_[i][2].resize(tsize);
|
||||
for (int j=0, k=0; j<tsize; ++j) {
|
||||
undersat_oil_tables_[i][0][j] = pvto[region_number][i][++k]; // p
|
||||
undersat_oil_tables_[i][1][j] = 1.0/pvto[region_number][i][++k]; // 1/Bo
|
||||
undersat_oil_tables_[i][2][j] = pvto[region_number][i][++k]; // mu_o
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Fill in additional entries in undersaturated tables by interpolating/extrapolating 1/Bo and mu_o ...
|
||||
int iPrev = -1;
|
||||
int iNext = 1;
|
||||
while (undersat_oil_tables_[iNext][0].size() < 2) {
|
||||
++iNext;
|
||||
}
|
||||
ASSERT(iNext < sz);
|
||||
for (int i=0; i<sz; ++i) {
|
||||
if (undersat_oil_tables_[i][0].size() > 1) {
|
||||
iPrev = i;
|
||||
continue;
|
||||
}
|
||||
|
||||
bool flagPrev = (iPrev >= 0);
|
||||
bool flagNext = true;
|
||||
if (iNext < i) {
|
||||
iPrev = iNext;
|
||||
flagPrev = true;
|
||||
iNext = i+1;
|
||||
while (undersat_oil_tables_[iNext][0].size() < 2) {
|
||||
++iNext;
|
||||
}
|
||||
}
|
||||
double slopePrevBinv = 0.0;
|
||||
double slopePrevVisc = 0.0;
|
||||
double slopeNextBinv = 0.0;
|
||||
double slopeNextVisc = 0.0;
|
||||
while (flagPrev || flagNext) {
|
||||
double pressure0 = undersat_oil_tables_[i][0].back();
|
||||
double pressure = 1.0e47;
|
||||
if (flagPrev) {
|
||||
std::vector<double>::iterator itPrev = upper_bound(undersat_oil_tables_[iPrev][0].begin(),
|
||||
undersat_oil_tables_[iPrev][0].end(),pressure0+1.);
|
||||
if (itPrev == undersat_oil_tables_[iPrev][0].end()) {
|
||||
--itPrev; // Extrapolation ...
|
||||
} else if (itPrev == undersat_oil_tables_[iPrev][0].begin()) {
|
||||
++itPrev;
|
||||
}
|
||||
if (itPrev == undersat_oil_tables_[iPrev][0].end()-1) {
|
||||
flagPrev = false; // Last data set for "prev" ...
|
||||
}
|
||||
double dPPrev = *itPrev - *(itPrev-1);
|
||||
pressure = *itPrev;
|
||||
int index = int(itPrev - undersat_oil_tables_[iPrev][0].begin());
|
||||
slopePrevBinv = (undersat_oil_tables_[iPrev][1][index] - undersat_oil_tables_[iPrev][1][index-1])/dPPrev;
|
||||
slopePrevVisc = (undersat_oil_tables_[iPrev][2][index] - undersat_oil_tables_[iPrev][2][index-1])/dPPrev;
|
||||
}
|
||||
if (flagNext) {
|
||||
std::vector<double>::iterator itNext = upper_bound(undersat_oil_tables_[iNext][0].begin(),
|
||||
undersat_oil_tables_[iNext][0].end(),pressure0+1.);
|
||||
if (itNext == undersat_oil_tables_[iNext][0].end()) {
|
||||
--itNext; // Extrapolation ...
|
||||
} else if (itNext == undersat_oil_tables_[iNext][0].begin()) {
|
||||
++itNext;
|
||||
}
|
||||
if (itNext == undersat_oil_tables_[iNext][0].end()-1) {
|
||||
flagNext = false; // Last data set for "next" ...
|
||||
}
|
||||
double dPNext = *itNext - *(itNext-1);
|
||||
if (flagPrev) {
|
||||
pressure = std::min(pressure,*itNext);
|
||||
} else {
|
||||
pressure = *itNext;
|
||||
}
|
||||
int index = int(itNext - undersat_oil_tables_[iNext][0].begin());
|
||||
slopeNextBinv = (undersat_oil_tables_[iNext][1][index] - undersat_oil_tables_[iNext][1][index-1])/dPNext;
|
||||
slopeNextVisc = (undersat_oil_tables_[iNext][2][index] - undersat_oil_tables_[iNext][2][index-1])/dPNext;
|
||||
}
|
||||
double dP = pressure - pressure0;
|
||||
if (iPrev >= 0) {
|
||||
double w = (saturated_oil_table_[3][i] - saturated_oil_table_[3][iPrev]) /
|
||||
(saturated_oil_table_[3][iNext] - saturated_oil_table_[3][iPrev]);
|
||||
undersat_oil_tables_[i][0].push_back(pressure0+dP);
|
||||
undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back() +
|
||||
dP*(slopePrevBinv+w*(slopeNextBinv-slopePrevBinv)));
|
||||
undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back() +
|
||||
dP*(slopePrevVisc+w*(slopeNextVisc-slopePrevVisc)));
|
||||
} else {
|
||||
undersat_oil_tables_[i][0].push_back(pressure0+dP);
|
||||
undersat_oil_tables_[i][1].push_back(undersat_oil_tables_[i][1].back()+dP*slopeNextBinv);
|
||||
undersat_oil_tables_[i][2].push_back(undersat_oil_tables_[i][2].back()+dP*slopeNextVisc);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
// Destructor
|
||||
MiscibilityLiveOil::~MiscibilityLiveOil()
|
||||
{
|
||||
}
|
||||
|
||||
double MiscibilityLiveOil::getViscosity(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
return miscible_oil(press, surfvol, 2, false);
|
||||
}
|
||||
|
||||
void MiscibilityLiveOil::getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = miscible_oil(pressures[i][phase], surfvol[i], 2, false);
|
||||
}
|
||||
}
|
||||
|
||||
// Dissolved gas-oil ratio
|
||||
double MiscibilityLiveOil::R(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
return evalR(press, surfvol);
|
||||
}
|
||||
|
||||
void MiscibilityLiveOil::R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = evalR(pressures[i][phase], surfvol[i]);
|
||||
}
|
||||
}
|
||||
|
||||
// Dissolved gas-oil ratio derivative
|
||||
double MiscibilityLiveOil::dRdp(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
double R = linearInterpolationExtrap(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3], press);
|
||||
double maxR = surfvol[Vapour]/surfvol[Liquid];
|
||||
if (R < maxR ) { // Saturated case
|
||||
return linearInterpolDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],
|
||||
press);
|
||||
} else {
|
||||
return 0.0; // Undersaturated case
|
||||
}
|
||||
}
|
||||
|
||||
void MiscibilityLiveOil::dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
int num = pressures.size();
|
||||
output_R.resize(num);
|
||||
output_dRdp.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
evalRDeriv(pressures[i][phase], surfvol[i], output_R[i], output_dRdp[i]);
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityLiveOil::B(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
return evalB(press, surfvol);
|
||||
}
|
||||
|
||||
void MiscibilityLiveOil::B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
int num = pressures.size();
|
||||
output.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output[i] = evalB(pressures[i][phase], surfvol[i]);
|
||||
}
|
||||
}
|
||||
|
||||
double MiscibilityLiveOil::dBdp(int /*region*/, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
// if (surfvol[Liquid] == 0.0) return 0.0; // To handle no-oil case.
|
||||
double Bo = evalB(press, surfvol); // \TODO check if we incur virtual call overhead here.
|
||||
return -Bo*Bo*miscible_oil(press, surfvol, 1, true);
|
||||
}
|
||||
|
||||
void MiscibilityLiveOil::dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const
|
||||
{
|
||||
ASSERT(pressures.size() == surfvol.size());
|
||||
B(pressures, surfvol, phase, output_B);
|
||||
int num = pressures.size();
|
||||
output_dBdp.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output_dBdp[i] = dBdp(0, pressures[i][phase], surfvol[i]); // \TODO Speedup here by using already evaluated B.
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double MiscibilityLiveOil::evalR(double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
if (surfvol[Vapour] == 0.0) {
|
||||
return 0.0;
|
||||
}
|
||||
double R = linearInterpolationExtrap(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3], press);
|
||||
double maxR = surfvol[Vapour]/surfvol[Liquid];
|
||||
if (R < maxR ) { // Saturated case
|
||||
return R;
|
||||
} else {
|
||||
return maxR; // Undersaturated case
|
||||
}
|
||||
}
|
||||
|
||||
void MiscibilityLiveOil::evalRDeriv(const double press, const surfvol_t& surfvol,
|
||||
double& R, double& dRdp) const
|
||||
{
|
||||
if (surfvol[Vapour] == 0.0) {
|
||||
R = 0.0;
|
||||
dRdp = 0.0;
|
||||
return;
|
||||
}
|
||||
R = linearInterpolationExtrap(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3], press);
|
||||
double maxR = surfvol[Vapour]/surfvol[Liquid];
|
||||
if (R < maxR ) {
|
||||
// Saturated case
|
||||
dRdp = linearInterpolDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],
|
||||
press);
|
||||
} else {
|
||||
// Undersaturated case
|
||||
R = maxR;
|
||||
dRdp = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
double MiscibilityLiveOil::evalB(double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
// if (surfvol[Liquid] == 0.0) return 1.0; // To handle no-oil case.
|
||||
return 1.0/miscible_oil(press, surfvol, 1, false);
|
||||
}
|
||||
|
||||
|
||||
void MiscibilityLiveOil::evalBDeriv(const double press, const surfvol_t& surfvol,
|
||||
double& B, double& dBdp) const
|
||||
{
|
||||
B = evalB(press, surfvol);
|
||||
dBdp = -B*B*miscible_oil(press, surfvol, 1, true);
|
||||
}
|
||||
|
||||
|
||||
double MiscibilityLiveOil::miscible_oil(double press, const surfvol_t& surfvol,
|
||||
int item, bool deriv) const
|
||||
{
|
||||
int section;
|
||||
double R = linearInterpolationExtrap(saturated_oil_table_[0],
|
||||
saturated_oil_table_[3],
|
||||
press, section);
|
||||
double maxR = (surfvol[Liquid] == 0.0) ? 0.0 : surfvol[Vapour]/surfvol[Liquid];
|
||||
if (deriv) {
|
||||
if (R < maxR ) { // Saturated case
|
||||
return linearInterpolDerivative(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
int is = tableIndex(saturated_oil_table_[3], maxR);
|
||||
double w = (maxR - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
ASSERT(undersat_oil_tables_[is][0].size() >= 2);
|
||||
ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolDerivative(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolDerivative(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
} else {
|
||||
if (R < maxR ) { // Saturated case
|
||||
return linearInterpolationExtrap(saturated_oil_table_[0],
|
||||
saturated_oil_table_[item],
|
||||
press);
|
||||
} else { // Undersaturated case
|
||||
// Interpolate between table sections
|
||||
int is = tableIndex(saturated_oil_table_[3], maxR);
|
||||
double w = (maxR - saturated_oil_table_[3][is]) /
|
||||
(saturated_oil_table_[3][is+1] - saturated_oil_table_[3][is]);
|
||||
ASSERT(undersat_oil_tables_[is][0].size() >= 2);
|
||||
ASSERT(undersat_oil_tables_[is+1][0].size() >= 2);
|
||||
double val1 =
|
||||
linearInterpolationExtrap(undersat_oil_tables_[is][0],
|
||||
undersat_oil_tables_[is][item],
|
||||
press);
|
||||
double val2 =
|
||||
linearInterpolationExtrap(undersat_oil_tables_[is+1][0],
|
||||
undersat_oil_tables_[is+1][item],
|
||||
press);
|
||||
double val = val1 + w*(val2 - val1);
|
||||
return val;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
97
opm/core/fluid/blackoil/MiscibilityLiveOil.hpp
Normal file
97
opm/core/fluid/blackoil/MiscibilityLiveOil.hpp
Normal file
@ -0,0 +1,97 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityLiveOil.hpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:08:09 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef SINTEF_MISCIBILITYLIVEOIL_HEADER
|
||||
#define SINTEF_MISCIBILITYLIVEOIL_HEADER
|
||||
|
||||
/** Class for miscible live oil.
|
||||
* Detailed description.
|
||||
*/
|
||||
|
||||
#include "MiscibilityProps.hpp"
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
class MiscibilityLiveOil : public MiscibilityProps
|
||||
{
|
||||
public:
|
||||
typedef std::vector<std::vector<std::vector<double> > > table_t;
|
||||
|
||||
MiscibilityLiveOil(const table_t& pvto);
|
||||
virtual ~MiscibilityLiveOil();
|
||||
|
||||
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double R (int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double B (int region, double press, const surfvol_t& surfvol) const;
|
||||
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const;
|
||||
|
||||
virtual void getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const;
|
||||
virtual void R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const;
|
||||
virtual void dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const;
|
||||
|
||||
protected:
|
||||
double evalR(double press, const surfvol_t& surfvol) const;
|
||||
void evalRDeriv(double press, const surfvol_t& surfvol, double& R, double& dRdp) const;
|
||||
double evalB(double press, const surfvol_t& surfvol) const;
|
||||
void evalBDeriv(double press, const surfvol_t& surfvol, double& B, double& dBdp) const;
|
||||
|
||||
// item: 1=B 2=mu;
|
||||
double miscible_oil(double press, const surfvol_t& surfvol, int item,
|
||||
bool deriv = false) const;
|
||||
|
||||
// PVT properties of live oil (with dissolved gas)
|
||||
std::vector<std::vector<double> > saturated_oil_table_;
|
||||
std::vector<std::vector<std::vector<double> > > undersat_oil_tables_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // SINTEF_MISCIBILITYLIVEOIL_HEADER
|
||||
|
||||
52
opm/core/fluid/blackoil/MiscibilityProps.cpp
Normal file
52
opm/core/fluid/blackoil/MiscibilityProps.cpp
Normal file
@ -0,0 +1,52 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityProps.cpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:05:05 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "MiscibilityProps.hpp"
|
||||
|
||||
using namespace std;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
//------------------------------------------------------------------------
|
||||
// Member functions
|
||||
//-------------------------------------------------------------------------
|
||||
|
||||
/// Constructor
|
||||
MiscibilityProps::MiscibilityProps()
|
||||
{
|
||||
}
|
||||
|
||||
MiscibilityProps::~MiscibilityProps()
|
||||
{
|
||||
}
|
||||
|
||||
} // namespace Opm
|
||||
87
opm/core/fluid/blackoil/MiscibilityProps.hpp
Normal file
87
opm/core/fluid/blackoil/MiscibilityProps.hpp
Normal file
@ -0,0 +1,87 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityProps.hpp
|
||||
//
|
||||
// Created: Wed Feb 10 09:04:35 2010
|
||||
//
|
||||
// Author: Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// Revision: $Id$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef SINTEF_MISCIBILITYPROPS_HEADER
|
||||
#define SINTEF_MISCIBILITYPROPS_HEADER
|
||||
|
||||
/** Base class for properties of fluids and rocks.
|
||||
* Detailed description.
|
||||
*/
|
||||
|
||||
#include "BlackoilDefs.hpp"
|
||||
#include <vector>
|
||||
#include <tr1/array>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
|
||||
|
||||
class MiscibilityProps : public BlackoilDefs
|
||||
{
|
||||
public:
|
||||
typedef CompVec surfvol_t;
|
||||
|
||||
MiscibilityProps();
|
||||
virtual ~MiscibilityProps();
|
||||
|
||||
virtual double getViscosity(int region, double press, const surfvol_t& surfvol) const = 0;
|
||||
virtual double B (int region, double press, const surfvol_t& surfvol) const = 0;
|
||||
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const = 0;
|
||||
virtual double R (int region, double press, const surfvol_t& surfvol) const = 0;
|
||||
virtual double dRdp(int region, double press, const surfvol_t& surfvol) const = 0;
|
||||
|
||||
virtual void getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const = 0;
|
||||
virtual void B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const = 0;
|
||||
virtual void dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const = 0;
|
||||
virtual void R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output) const = 0;
|
||||
virtual void dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvol,
|
||||
int phase,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const = 0;
|
||||
};
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif // SINTEF_MISCIBILITYPROPS_HEADER
|
||||
|
||||
194
opm/core/fluid/blackoil/MiscibilityWater.hpp
Normal file
194
opm/core/fluid/blackoil/MiscibilityWater.hpp
Normal file
@ -0,0 +1,194 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: MiscibilityWater.hpp
|
||||
//
|
||||
// Created: Tue May 18 10:26:13 2010
|
||||
//
|
||||
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
|
||||
// Bjørn Spjelkavik <bsp@sintef.no>
|
||||
//
|
||||
// $Date$
|
||||
//
|
||||
// $Revision$
|
||||
//
|
||||
//===========================================================================
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPENRS_MISCIBILITYWATER_HEADER
|
||||
#define OPENRS_MISCIBILITYWATER_HEADER
|
||||
|
||||
#include <opm/core/fluid/blackoil/MiscibilityProps.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
// Forward declaration.
|
||||
class PVTW;
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
class MiscibilityWater : public MiscibilityProps
|
||||
{
|
||||
public:
|
||||
typedef std::vector<std::vector<double> > table_t;
|
||||
|
||||
MiscibilityWater(const table_t& pvtw)
|
||||
{
|
||||
const int region_number = 0;
|
||||
if (pvtw.size() != 1) {
|
||||
THROW("More than one PVD-region");
|
||||
}
|
||||
ref_press_ = pvtw[region_number][0];
|
||||
ref_B_ = pvtw[region_number][1];
|
||||
comp_ = pvtw[region_number][2];
|
||||
viscosity_ = pvtw[region_number][3];
|
||||
if (pvtw[region_number].size() > 4 && pvtw[region_number][4] != 0.0) {
|
||||
THROW("MiscibilityWater does not support 'viscosibility'.");
|
||||
}
|
||||
}
|
||||
|
||||
MiscibilityWater(double visc)
|
||||
: ref_press_(0.0),
|
||||
ref_B_(1.0),
|
||||
comp_(0.0),
|
||||
viscosity_(visc)
|
||||
{
|
||||
}
|
||||
|
||||
virtual ~MiscibilityWater()
|
||||
{
|
||||
}
|
||||
|
||||
virtual double getViscosity(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
return viscosity_;
|
||||
}
|
||||
|
||||
virtual void getViscosity(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.clear();
|
||||
output.resize(num, viscosity_);
|
||||
}
|
||||
|
||||
virtual double B(int /*region*/, double press, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
if (comp_) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = comp_*(press - ref_press_);
|
||||
return ref_B_/(1.0 + x + 0.5*x*x);
|
||||
} else {
|
||||
return ref_B_;
|
||||
}
|
||||
}
|
||||
|
||||
virtual void B(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int phase,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
if (comp_) {
|
||||
output.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
// Computing a polynomial approximation to the exponential.
|
||||
double x = comp_*(pressures[i][phase] - ref_press_);
|
||||
output[i] = ref_B_/(1.0 + x + 0.5*x*x);
|
||||
}
|
||||
} else {
|
||||
output.clear();
|
||||
output.resize(num, ref_B_);
|
||||
}
|
||||
}
|
||||
|
||||
virtual double dBdp(int region, double press, const surfvol_t& surfvol) const
|
||||
{
|
||||
if (comp_) {
|
||||
return -comp_*B(region, press, surfvol);
|
||||
} else {
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
virtual void dBdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>& surfvols,
|
||||
int phase,
|
||||
std::vector<double>& output_B,
|
||||
std::vector<double>& output_dBdp) const
|
||||
{
|
||||
B(pressures, surfvols, phase, output_B);
|
||||
int num = pressures.size();
|
||||
if (comp_) {
|
||||
output_dBdp.resize(num);
|
||||
#pragma omp parallel for
|
||||
for (int i = 0; i < num; ++i) {
|
||||
output_dBdp[i] = -comp_*output_B[i];
|
||||
}
|
||||
} else {
|
||||
output_dBdp.clear();
|
||||
output_dBdp.resize(num, 0.0);
|
||||
}
|
||||
}
|
||||
|
||||
virtual double R(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
virtual void R(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int,
|
||||
std::vector<double>& output) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output.clear();
|
||||
output.resize(num, 0.0);
|
||||
}
|
||||
|
||||
virtual double dRdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
virtual void dRdp(const std::vector<PhaseVec>& pressures,
|
||||
const std::vector<CompVec>&,
|
||||
int,
|
||||
std::vector<double>& output_R,
|
||||
std::vector<double>& output_dRdp) const
|
||||
{
|
||||
int num = pressures.size();
|
||||
output_R.clear();
|
||||
output_R.resize(num, 0.0);
|
||||
output_dRdp.clear();
|
||||
output_dRdp.resize(num, 0.0);
|
||||
}
|
||||
|
||||
private:
|
||||
double ref_press_;
|
||||
double ref_B_;
|
||||
double comp_;
|
||||
double viscosity_;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
#endif // OPENRS_MISCIBILITYWATER_HEADER
|
||||
260
opm/core/utility/UniformTableLinear.hpp
Normal file
260
opm/core/utility/UniformTableLinear.hpp
Normal file
@ -0,0 +1,260 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED
|
||||
#define OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED
|
||||
|
||||
#include <cmath>
|
||||
#include <exception>
|
||||
#include <vector>
|
||||
#include <utility>
|
||||
#include <iostream>
|
||||
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
#include <opm/core/utility/linearInterpolation.hpp>
|
||||
|
||||
namespace Dune {
|
||||
namespace utils {
|
||||
|
||||
|
||||
/// @brief This class uses linear interpolation to compute the value
|
||||
/// (and its derivative) of a function f sampled at uniform points.
|
||||
/// @tparam T the range type of the function (should be an algebraic ring type)
|
||||
template<typename T>
|
||||
class UniformTableLinear
|
||||
{
|
||||
public:
|
||||
/// @brief Default constructor.
|
||||
UniformTableLinear();
|
||||
|
||||
/// @brief Construct from vector of y-values.
|
||||
/// @param xmin the x value corresponding to the first y value.
|
||||
/// @param xmax the x value corresponding to the last y value.
|
||||
/// @param y_values vector of range values.
|
||||
UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const std::vector<T>& y_values);
|
||||
|
||||
/// @brief Construct from array of y-values.
|
||||
/// @param xmin the x value corresponding to the first y value.
|
||||
/// @param xmax the x value corresponding to the last y value.
|
||||
/// @param y_values array of range values.
|
||||
/// @param num_y_values the number of values in y_values.
|
||||
UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const T* y_values,
|
||||
int num_y_values);
|
||||
|
||||
/// @brief Get the domain.
|
||||
/// @return the domain as a pair of doubles.
|
||||
std::pair<double, double> domain();
|
||||
|
||||
/// @brief Rescale the domain.
|
||||
/// @param new_domain the new domain as a pair of doubles.
|
||||
void rescaleDomain(std::pair<double, double> new_domain);
|
||||
|
||||
/// @brief Evaluate the value at x.
|
||||
/// @param x a domain value
|
||||
/// @return f(x)
|
||||
double operator()(const double x) const;
|
||||
|
||||
/// @brief Evaluate the derivative at x.
|
||||
/// @param x a domain value
|
||||
/// @return f'(x)
|
||||
double derivative(const double x) const;
|
||||
|
||||
/// @brief Equality operator.
|
||||
/// @param other another UniformTableLinear.
|
||||
/// @return true if they are represented exactly alike.
|
||||
bool operator==(const UniformTableLinear& other) const;
|
||||
|
||||
/// @brief Policies for how to behave when trying to evaluate outside the domain.
|
||||
enum RangePolicy {Throw = 0, ClosestValue = 1, Extrapolate = 2};
|
||||
|
||||
/// @brief Sets the behavioural policy for evaluation to the left of the domain.
|
||||
/// @param rp the policy
|
||||
void setLeftPolicy(RangePolicy rp);
|
||||
|
||||
/// @brief Sets the behavioural policy for evaluation to the right of the domain.
|
||||
/// @param rp the policy
|
||||
void setRightPolicy(RangePolicy rp);
|
||||
|
||||
protected:
|
||||
double xmin_;
|
||||
double xmax_;
|
||||
double xdelta_;
|
||||
std::vector<T> y_values_;
|
||||
RangePolicy left_;
|
||||
RangePolicy right_;
|
||||
template <typename U>
|
||||
friend std::ostream& operator<<(std::ostream& os, const UniformTableLinear<U>& t);
|
||||
};
|
||||
|
||||
|
||||
// Member implementations.
|
||||
|
||||
template<typename T>
|
||||
inline
|
||||
UniformTableLinear<T>
|
||||
::UniformTableLinear()
|
||||
: left_(ClosestValue), right_(ClosestValue)
|
||||
{
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline
|
||||
UniformTableLinear<T>
|
||||
::UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const std::vector<T>& y_values)
|
||||
: xmin_(xmin), xmax_(xmax), y_values_(y_values),
|
||||
left_(ClosestValue), right_(ClosestValue)
|
||||
{
|
||||
ASSERT(xmax > xmin);
|
||||
ASSERT(y_values.size() > 1);
|
||||
xdelta_ = (xmax - xmin)/(y_values.size() - 1);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline
|
||||
UniformTableLinear<T>
|
||||
::UniformTableLinear(double xmin,
|
||||
double xmax,
|
||||
const T* y_values,
|
||||
int num_y_values)
|
||||
: xmin_(xmin), xmax_(xmax),
|
||||
y_values_(y_values, y_values + num_y_values),
|
||||
left_(ClosestValue), right_(ClosestValue)
|
||||
{
|
||||
ASSERT(xmax > xmin);
|
||||
ASSERT(y_values_.size() > 1);
|
||||
xdelta_ = (xmax - xmin)/(y_values_.size() - 1);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline std::pair<double, double>
|
||||
UniformTableLinear<T>
|
||||
::domain()
|
||||
{
|
||||
return std::make_pair(xmin_, xmax_);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline void
|
||||
UniformTableLinear<T>
|
||||
::rescaleDomain(std::pair<double, double> new_domain)
|
||||
{
|
||||
xmin_ = new_domain.first;
|
||||
xmax_ = new_domain.second;
|
||||
xdelta_ = (xmax_ - xmin_)/(y_values_.size() - 1);
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline double
|
||||
UniformTableLinear<T>
|
||||
::operator()(const double xparam) const
|
||||
{
|
||||
// Implements ClosestValue policy.
|
||||
double x = std::min(xparam, xmax_);
|
||||
x = std::max(x, xmin_);
|
||||
|
||||
// Lookup is easy since we are uniform in x.
|
||||
double pos = (x - xmin_)/xdelta_;
|
||||
double posi = std::floor(pos);
|
||||
int left = int(posi);
|
||||
if (left == int(y_values_.size()) - 1) {
|
||||
// We are at xmax_
|
||||
return y_values_.back();
|
||||
}
|
||||
double w = pos - posi;
|
||||
return (1.0 - w)*y_values_[left] + w*y_values_[left + 1];
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline double
|
||||
UniformTableLinear<T>
|
||||
::derivative(const double xparam) const
|
||||
{
|
||||
// Implements ClosestValue policy.
|
||||
double x = std::min(xparam, xmax_);
|
||||
x = std::max(x, xmin_);
|
||||
|
||||
// Lookup is easy since we are uniform in x.
|
||||
double pos = (x - xmin_)/xdelta_;
|
||||
double posi = std::floor(pos);
|
||||
int left = int(posi);
|
||||
if (left == int(y_values_.size()) - 1) {
|
||||
// We are at xmax_
|
||||
--left;
|
||||
}
|
||||
return (y_values_[left + 1] - y_values_[left])/xdelta_;
|
||||
}
|
||||
|
||||
|
||||
template<typename T>
|
||||
inline bool
|
||||
UniformTableLinear<T>
|
||||
::operator==(const UniformTableLinear<T>& other) const
|
||||
{
|
||||
return xmin_ == other.xmin_
|
||||
&& xdelta_ == other.xdelta_
|
||||
&& y_values_ == other.y_values_
|
||||
&& left_ == other.left_
|
||||
&& right_ == other.right_;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline void
|
||||
UniformTableLinear<T>
|
||||
::setLeftPolicy(RangePolicy rp)
|
||||
{
|
||||
if (rp != ClosestValue) {
|
||||
THROW("Only ClosestValue RangePolicy implemented.");
|
||||
}
|
||||
left_ = rp;
|
||||
}
|
||||
|
||||
template<typename T>
|
||||
inline void
|
||||
UniformTableLinear<T>
|
||||
::setRightPolicy(RangePolicy rp)
|
||||
{
|
||||
if (rp != ClosestValue) {
|
||||
THROW("Only ClosestValue RangePolicy implemented.");
|
||||
}
|
||||
right_ = rp;
|
||||
}
|
||||
|
||||
|
||||
template <typename T>
|
||||
inline std::ostream& operator<<(std::ostream& os, const UniformTableLinear<T>& t)
|
||||
{
|
||||
int n = t.y_values_.size();
|
||||
for (int i = 0; i < n; ++i) {
|
||||
double f = double(i)/double(n - 1);
|
||||
os << (1.0 - f)*t.xmin_ + f*t.xmax_
|
||||
<< " " << t.y_values_[i] << '\n';
|
||||
}
|
||||
return os;
|
||||
}
|
||||
|
||||
} // namespace utils
|
||||
} // namespace Dune
|
||||
|
||||
#endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED
|
||||
52
opm/core/utility/buildUniformMonotoneTable.hpp
Normal file
52
opm/core/utility/buildUniformMonotoneTable.hpp
Normal file
@ -0,0 +1,52 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED
|
||||
#define OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED
|
||||
|
||||
#include <opm/core/utility/MonotCubicInterpolator.hpp>
|
||||
#include <opm/core/utility/UniformTableLinear.hpp>
|
||||
|
||||
namespace Dune {
|
||||
namespace utils {
|
||||
|
||||
template <typename T>
|
||||
void buildUniformMonotoneTable(const std::vector<double>& xv,
|
||||
const std::vector<T>& yv,
|
||||
const int samples,
|
||||
UniformTableLinear<T>& table)
|
||||
{
|
||||
MonotCubicInterpolator interp(xv, yv);
|
||||
std::vector<T> uniform_yv(samples);
|
||||
double xmin = xv[0];
|
||||
double xmax = xv.back();
|
||||
for (int i = 0; i < samples; ++i) {
|
||||
double w = double(i)/double(samples - 1);
|
||||
double x = (1.0 - w)*xmin + w*xmax;
|
||||
uniform_yv[i] = interp(x);
|
||||
}
|
||||
table = UniformTableLinear<T>(xmin, xmax, uniform_yv);
|
||||
}
|
||||
|
||||
} // namespace utils
|
||||
} // namespace Dune
|
||||
|
||||
|
||||
|
||||
#endif // OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED
|
||||
90
opm/core/utility/linearInterpolation.hpp
Normal file
90
opm/core/utility/linearInterpolation.hpp
Normal file
@ -0,0 +1,90 @@
|
||||
//===========================================================================
|
||||
//
|
||||
// File: linearInterpolation.hpp
|
||||
//
|
||||
// Created: Tue Sep 9 12:49:39 2008
|
||||
//
|
||||
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
|
||||
//
|
||||
// $Date$
|
||||
//
|
||||
// $Revision$
|
||||
//
|
||||
//===========================================================================
|
||||
|
||||
/*
|
||||
Copyright 2009, 2010 SINTEF ICT, Applied Mathematics.
|
||||
Copyright 2009, 2010 Statoil ASA.
|
||||
|
||||
This file is part of The Open Reservoir Simulator Project (OpenRS).
|
||||
|
||||
OpenRS 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.
|
||||
|
||||
OpenRS 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 OpenRS. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#ifndef OPENRS_LINEARINTERPOLATION_HEADER
|
||||
#define OPENRS_LINEARINTERPOLATION_HEADER
|
||||
|
||||
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
|
||||
namespace Dune
|
||||
{
|
||||
|
||||
/** Linear interpolation.
|
||||
* Given an increasing vector xv of parameter values and
|
||||
* a vector yv of point values of the same size,
|
||||
* the function returns ...
|
||||
*/
|
||||
template <typename T>
|
||||
T linearInterpolation(const std::vector<double>& xv,
|
||||
const std::vector<T>& yv,
|
||||
double x)
|
||||
{
|
||||
std::vector<double>::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x);
|
||||
int lb_ix = lb - xv.begin();
|
||||
if (lb_ix == 0) {
|
||||
return yv[0];
|
||||
} else if (lb_ix == int(xv.size())) {
|
||||
return yv.back();
|
||||
} else {
|
||||
double w = (x - xv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]);
|
||||
return (1.0 - w)*yv[lb_ix - 1] + w*yv[lb_ix];
|
||||
}
|
||||
}
|
||||
|
||||
/// @brief
|
||||
/// @todo Doc me!
|
||||
/// @tparam
|
||||
/// @param
|
||||
/// @return
|
||||
template <typename T>
|
||||
T linearInterpolationDerivative(const std::vector<double>& xv,
|
||||
const std::vector<T>& yv,
|
||||
double x)
|
||||
{
|
||||
double epsilon = 1e-4; // @@ Ad hoc, should choose based on xv.
|
||||
double x_low = std::max(xv[0], x - epsilon);
|
||||
double x_high = std::min(xv.back(), x + epsilon);
|
||||
T low = linearInterpolation(xv, yv, x_low);
|
||||
T high = linearInterpolation(xv, yv, x_high);
|
||||
return (high - low)/(x_high - x_low);
|
||||
}
|
||||
|
||||
|
||||
} // namespace Dune
|
||||
|
||||
|
||||
|
||||
#endif // OPENRS_LINEARINTERPOLATION_HEADER
|
||||
@ -5,16 +5,13 @@ $(BOOST_CPPFLAGS)
|
||||
LDFLAGS = $(BOOST_LDFLAGS)
|
||||
|
||||
|
||||
LDADD = \
|
||||
$(top_builddir)/libopmcore.la \
|
||||
$(BOOST_FILESYSTEM_LIB) \
|
||||
$(BOOST_SYSTEM_LIB) \
|
||||
$(BOOST_DATE_TIME_LIB) \
|
||||
$(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
|
||||
$(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)
|
||||
LDADD = $(top_builddir)/libopmcore.la
|
||||
|
||||
|
||||
noinst_PROGRAMS = \
|
||||
bo_fluid_p_and_z_deps \
|
||||
bo_fluid_pressuredeps \
|
||||
bo_fluid_test \
|
||||
monotcubicinterpolator_test \
|
||||
param_test \
|
||||
sparsetable_test \
|
||||
@ -23,6 +20,9 @@ unit_test \
|
||||
test_readvector \
|
||||
test_sf2p
|
||||
|
||||
bo_fluid_p_and_z_deps_SOURCES = bo_fluid_p_and_z_deps.cpp
|
||||
bo_fluid_pressuredeps_SOURCES = bo_fluid_pressuredeps.cpp
|
||||
bo_fluid_test_SOURCES = bo_fluid_test.cpp
|
||||
monotcubicinterpolator_test_SOURCES = monotcubicinterpolator_test.cpp
|
||||
param_test_SOURCES = param_test.cpp
|
||||
sparsetable_test_SOURCES = sparsetable_test.cpp
|
||||
|
||||
132
tests/bo_fluid_p_and_z_deps.cpp
Normal file
132
tests/bo_fluid_p_and_z_deps.cpp
Normal file
@ -0,0 +1,132 @@
|
||||
/*
|
||||
Copyright 2011 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/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/fluid/BlackoilFluid.hpp>
|
||||
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
std::cout << "%{\n";
|
||||
|
||||
// Parameters.
|
||||
Dune::parameter::ParameterGroup param(argc, argv);
|
||||
|
||||
// Parser.
|
||||
std::string ecl_file = param.get<std::string>("filename");
|
||||
Dune::EclipseGridParser parser(ecl_file);
|
||||
// Look at the BlackoilFluid behaviour
|
||||
Opm::BlackoilFluid fluid;
|
||||
fluid.init(parser);
|
||||
Opm::BlackoilFluid::CompVec z0(0.0);
|
||||
z0[Opm::BlackoilFluid::Water] = param.getDefault("z_w", 0.0);
|
||||
z0[Opm::BlackoilFluid::Oil] = param.getDefault("z_o", 1.0);
|
||||
z0[Opm::BlackoilFluid::Gas] = param.getDefault("z_g", 0.0);
|
||||
int num_pts_p = param.getDefault("num_pts_p", 41);
|
||||
int num_pts_z = param.getDefault("num_pts_z", 51);
|
||||
double min_press = param.getDefault("min_press", 1e7);
|
||||
double max_press = param.getDefault("max_press", 3e7);
|
||||
int changing_component = param.getDefault("changing_component", int(Opm::BlackoilFluid::Gas));
|
||||
double min_z = param.getDefault("min_z", 0.0);
|
||||
double max_z = param.getDefault("max_z", 500.0);
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
int variable = param.getDefault("variable", 0);
|
||||
#else
|
||||
int variable = param.getDefault("variable", 2);
|
||||
#endif
|
||||
Opm::BlackoilFluid::CompVec z = z0;
|
||||
std::cout << "%}\n"
|
||||
<< "data = [\n";
|
||||
|
||||
for (int i = 0; i < num_pts_p; ++i) {
|
||||
double pfactor = num_pts_p < 2 ? 0.0 : double(i)/double(num_pts_p - 1);
|
||||
double p = (1.0 - pfactor)*min_press + pfactor*max_press;
|
||||
for (int j = 0; j < num_pts_z; ++j) {
|
||||
double zfactor = num_pts_z < 2 ? 0.0 : double(j)/double(num_pts_z - 1);
|
||||
z[changing_component] = (1.0 - zfactor)*min_z + zfactor*max_z;
|
||||
// std::cout << p << ' ' << z << '\n';
|
||||
Opm::BlackoilFluid::FluidState state = fluid.computeState(Opm::BlackoilFluid::PhaseVec(p), z);
|
||||
std::cout.precision(6);
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
double var = 0.0;
|
||||
switch (variable) {
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
case 0:
|
||||
var = state.total_compressibility_;
|
||||
break;
|
||||
case 1:
|
||||
var = state.experimental_term_;
|
||||
break;
|
||||
#endif
|
||||
case 2:
|
||||
var = state.saturation_[0];
|
||||
break;
|
||||
case 3:
|
||||
var = state.saturation_[1];
|
||||
break;
|
||||
case 4:
|
||||
var = state.saturation_[2];
|
||||
break;
|
||||
case 5:
|
||||
var = state.formation_volume_factor_[0];
|
||||
break;
|
||||
case 6:
|
||||
var = state.formation_volume_factor_[1];
|
||||
break;
|
||||
case 7:
|
||||
var = state.formation_volume_factor_[2];
|
||||
break;
|
||||
case 8:
|
||||
var = state.solution_factor_[0];
|
||||
break;
|
||||
case 9:
|
||||
var = state.solution_factor_[1];
|
||||
break;
|
||||
case 10:
|
||||
var = state.solution_factor_[2];
|
||||
break;
|
||||
default:
|
||||
THROW("Unknown varable specification: " << variable);
|
||||
break;
|
||||
}
|
||||
std::cout << var << ' ';
|
||||
}
|
||||
std::cout << '\n';
|
||||
}
|
||||
std::cout << "];\n\n"
|
||||
<< "paxis = [\n";
|
||||
for (int i = 0; i < num_pts_p; ++i) {
|
||||
double pfactor = double(i)/double(num_pts_p - 1);
|
||||
double p = (1.0 - pfactor)*min_press + pfactor*max_press;
|
||||
std::cout << p << '\n';
|
||||
}
|
||||
std::cout << "];\n\n"
|
||||
<< "zaxis = [\n";
|
||||
for (int j = 0; j < num_pts_z; ++j) {
|
||||
double zfactor = double(j)/double(num_pts_z - 1);
|
||||
std::cout << (1.0 - zfactor)*min_z + zfactor*max_z << '\n';
|
||||
}
|
||||
std::cout << "];\n";
|
||||
}
|
||||
|
||||
108
tests/bo_fluid_pressuredeps.cpp
Normal file
108
tests/bo_fluid_pressuredeps.cpp
Normal file
@ -0,0 +1,108 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/fluid/BlackoilFluid.hpp>
|
||||
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
void printCompressibilityTerms(double p, const Opm::BlackoilFluid::FluidState& state)
|
||||
{
|
||||
std::cout.precision(6);
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << p << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << state.total_compressibility_ << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << totmob << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << state.total_phase_volume_density_ << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << state.experimental_term_ << '\n';
|
||||
}
|
||||
#endif
|
||||
|
||||
void printSatsEtc(double p, const Opm::BlackoilFluid::FluidState& state)
|
||||
{
|
||||
std::cout.precision(6);
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << p << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << state.saturation_[0] << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << state.saturation_[1] << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << state.saturation_[2] << " ";
|
||||
std::cout.width(15);
|
||||
std::cout.fill(' ');
|
||||
std::cout << state.total_phase_volume_density_ << '\n';
|
||||
}
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
// Parameters.
|
||||
Dune::parameter::ParameterGroup param(argc, argv);
|
||||
|
||||
// Parser.
|
||||
std::string ecl_file = param.get<std::string>("filename");
|
||||
Dune::EclipseGridParser parser(ecl_file);
|
||||
|
||||
// Look at the BlackoilFluid behaviour
|
||||
Opm::BlackoilFluid fluid;
|
||||
fluid.init(parser);
|
||||
Opm::BlackoilFluid::CompVec z(0.0);
|
||||
z[Opm::BlackoilFluid::Water] = param.getDefault("z_w", 0.0);
|
||||
z[Opm::BlackoilFluid::Oil] = param.getDefault("z_o", 1.0);
|
||||
z[Opm::BlackoilFluid::Gas] = param.getDefault("z_g", 0.0);
|
||||
int num_pts = param.getDefault("num_pts", 41);
|
||||
double min_press = param.getDefault("min_press", 1e7);
|
||||
double max_press = param.getDefault("max_press", 3e7);
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
bool print_compr = param.getDefault("print_compr", true);
|
||||
#endif
|
||||
for (int i = 0; i < num_pts; ++i) {
|
||||
double factor = double(i)/double(num_pts - 1);
|
||||
double p = (1.0 - factor)*min_press + factor*max_press;
|
||||
Opm::BlackoilFluid::FluidState state = fluid.computeState(Opm::BlackoilFluid::PhaseVec(p), z);
|
||||
double totmob = 0.0;
|
||||
for (int phase = 0; phase < Opm::BlackoilFluid::numPhases; ++phase) {
|
||||
totmob += state.mobility_[phase];
|
||||
}
|
||||
#ifdef COMPUTE_OLD_TERMS
|
||||
if (print_compr) {
|
||||
printCompressibilityTerms(p, state);
|
||||
} else {
|
||||
printSatsEtc(p, state);
|
||||
}
|
||||
#else
|
||||
printSatsEtc(p, state);
|
||||
#endif
|
||||
}
|
||||
}
|
||||
53
tests/bo_fluid_test.cpp
Normal file
53
tests/bo_fluid_test.cpp
Normal file
@ -0,0 +1,53 @@
|
||||
/*
|
||||
Copyright 2010 SINTEF ICT, Applied Mathematics.
|
||||
|
||||
This file is part of the Open Porous Media project (OPM).
|
||||
|
||||
OPM is free software: you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
the Free Software Foundation, either version 3 of the License, or
|
||||
(at your option) any later version.
|
||||
|
||||
OPM is distributed in the hope that it will be useful,
|
||||
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
GNU General Public License for more details.
|
||||
|
||||
You should have received a copy of the GNU General Public License
|
||||
along with OPM. If not, see <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "config.h"
|
||||
|
||||
#include <opm/core/utility/parameters/ParameterGroup.hpp>
|
||||
#include <opm/core/eclipse/EclipseGridParser.hpp>
|
||||
#include <opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp>
|
||||
|
||||
|
||||
int main(int argc, char** argv)
|
||||
{
|
||||
// Parameters.
|
||||
Dune::parameter::ParameterGroup param(argc, argv);
|
||||
|
||||
// Parser.
|
||||
std::string ecl_file = param.get<std::string>("filename");
|
||||
Dune::EclipseGridParser parser(ecl_file);
|
||||
|
||||
// Test the FluidMatrixInteractionBlackoil class.
|
||||
Opm::FluidMatrixInteractionBlackoilParams<double> fluid_params;
|
||||
fluid_params.init(parser);
|
||||
typedef Opm::FluidMatrixInteractionBlackoil<double> Law;
|
||||
Law::PhaseVec s, kr;
|
||||
const double temp = 300; // [K]
|
||||
int num = 41;
|
||||
for (int i = 0; i < num; ++i) {
|
||||
s[Law::Aqua] = 0.0;
|
||||
s[Law::Liquid] = double(i)/double(num - 1);
|
||||
s[Law::Vapour] = 1.0 - s[Law::Aqua] - s[Law::Liquid];
|
||||
Law::kr(kr, fluid_params, s, temp);
|
||||
std::cout.width(6);
|
||||
std::cout.fill(' ');
|
||||
std::cout << s[Law::Liquid] << " " << kr[0] << ' ' << kr[1] << ' ' << kr[2] << '\n';
|
||||
}
|
||||
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user