diff --git a/.hgignore b/.hgignore index 53e2d2a8..e684c047 100644 --- a/.hgignore +++ b/.hgignore @@ -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 diff --git a/.hgtags b/.hgtags new file mode 100644 index 00000000..df0776bb --- /dev/null +++ b/.hgtags @@ -0,0 +1 @@ +b3cc85c2b5ab3f5b283244624d06a1afab9f8d67 Version-InitialImport-svn935 diff --git a/Makefile.am b/Makefile.am index e7cbba9e..dd5f74ce 100644 --- a/Makefile.am +++ b/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 \ diff --git a/examples/Makefile.am b/examples/Makefile.am index 204e23df..c41ea05b 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -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 diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index 12a36ff7..bafd23bc 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -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 sb_ptr = createSpecialField(is, keyword); - if (sb_ptr) { - specialmap[keyword] = sb_ptr; + map >::iterator pos = + specialmap.find(keyword); + if (pos == specialmap.end()) { + std::tr1::shared_ptr 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; diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index 59899633..c7ff48e4 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -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. diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index 6cc76fc4..5cd51488 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -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 tstep_; + + virtual std::string name() const {return std::string("TSTEP");} + + virtual void read(std::istream& is) + { + std::vector 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(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) diff --git a/opm/core/fluid/BlackoilFluid.hpp b/opm/core/fluid/BlackoilFluid.hpp new file mode 100644 index 00000000..6aa3142d --- /dev/null +++ b/opm/core/fluid/BlackoilFluid.hpp @@ -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 . +*/ + +#ifndef OPM_BLACKOILFLUID_HEADER_INCLUDED +#define OPM_BLACKOILFLUID_HEADER_INCLUDED + + +#include +#include +#include +#include +#include + + +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& 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::kr(state.relperm_, fmi_params_, state.saturation_, state.temperature_); + FluidMatrixInteractionBlackoil::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 + void computeBAndR(States& states) const + { + const std::vector& p = states.phase_pressure; + const std::vector& z = states.surface_volume_density; + std::vector& B = states.formation_volume_factor; + std::vector& R = states.solution_factor; + pvt_.B(p, z, B); + pvt_.R(p, z, R); + } + + /// Input: p, z + /// Output: B, R, mu + template + void computePvtNoDerivs(States& states) const + { + computeBAndR(states); + const std::vector& p = states.phase_pressure; + const std::vector& z = states.surface_volume_density; + std::vector& mu = states.viscosity; + pvt_.getViscosity(p, z, mu); + } + + /// Input: p, z + /// Output: B, dB/dp, R, dR/dp, mu + template + void computePvt(States& states) const + { + const std::vector& p = states.phase_pressure; + const std::vector& z = states.surface_volume_density; + std::vector& B = states.formation_volume_factor; + std::vector& dB = states.formation_volume_factor_deriv; + std::vector& R = states.solution_factor; + std::vector& dR = states.solution_factor_deriv; + std::vector& 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 + 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 + 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 + 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::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 + 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::kr(kr, fmi_params_, s, 300.0); + FluidMatrixInteractionBlackoil::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 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 + 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 surface_volume_density; // z + std::vector phase_pressure; // p + + // Variables from PVT functions. + std::vector formation_volume_factor; // B + std::vector solution_factor; // R + + // Variables computed from PVT data. + // The A matrices are all in Fortran order (or, equivalently, + // we store the transposes). + std::vector state_matrix; // A' = (RB^{-1})' + + // Variables computed from saturation. + std::vector mobility; // lambda + std::vector mobility_deriv; // dlambda/ds + + // Gravity and/or capillary pressure potential differences. + std::vector gravity_potential; // (\rho g \delta z)-ish contribution per face + }; + + + struct AllFluidData : public BlackoilDefs + { + // Per-cell data + AllFluidStates cell_data; + std::vector voldiscr; + std::vector relvoldiscr; + + // Per-face data. + FaceFluidData face_data; + + public: + template + void computeNew(const Grid& grid, + const Rock& rock, + const BlackoilFluid& fluid, + const typename Grid::Vector gravity, + const std::vector& cell_pressure, + const std::vector& face_pressure, + const std::vector& 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 + void computeUpwindProperties(const Grid& grid, + const BlackoilFluid& fluid, + const typename Grid::Vector gravity, + const std::vector& cell_pressure, + const std::vector& face_pressure, + const std::vector& 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 diff --git a/opm/core/fluid/blackoil/BlackoilDefs.hpp b/opm/core/fluid/blackoil/BlackoilDefs.hpp new file mode 100644 index 00000000..08e44cb8 --- /dev/null +++ b/opm/core/fluid/blackoil/BlackoilDefs.hpp @@ -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 . +*/ + +#ifndef OPM_BLACKOILDEFS_HEADER_INCLUDED +#define OPM_BLACKOILDEFS_HEADER_INCLUDED + + +#include +#include +#include + +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 + 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 + void assign(const U& elem) + { + for (int i = 0; i < N; ++i) { + data_[i] = elem; + } + } + private: + T data_[N]; + }; + template + class SmallMat + { + public: + SmallMat() + {} + SmallMat(const T& elem) + { data_.assign(elem); } + SmallMat& operator=(const T& elem) + { data_.assign(elem); return *this; } + typedef SmallVec RowType; + const RowType& operator[](int i) const + { return data_[i]; } + RowType& operator[](int i) + { return data_[i]; } + private: + SmallVec data_; + }; + + typedef double Scalar; + typedef SmallVec CompVec; + typedef SmallVec PhaseVec; + BOOST_STATIC_ASSERT(int(numComponents) == int(numPhases)); + typedef SmallMat PhaseToCompMatrix; + typedef SmallMat 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 diff --git a/opm/core/fluid/blackoil/BlackoilPVT.cpp b/opm/core/fluid/blackoil/BlackoilPVT.cpp new file mode 100644 index 00000000..076116a8 --- /dev/null +++ b/opm/core/fluid/blackoil/BlackoilPVT.cpp @@ -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 . +*/ + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace Dune; + +namespace Opm +{ + + + void BlackoilPVT::init(const Dune::EclipseGridParser& parser) + { + typedef std::vector > > 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& 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& pressures, + const std::vector& surfvol, + std::vector& 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& pressures, + const std::vector& surfvol, + std::vector& 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& pressures, + const std::vector& surfvol, + std::vector& output_B, + std::vector& 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& pressures, + const std::vector& surfvol, + std::vector& 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& pressures, + const std::vector& surfvol, + std::vector& output_R, + std::vector& 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 diff --git a/opm/core/fluid/blackoil/BlackoilPVT.hpp b/opm/core/fluid/blackoil/BlackoilPVT.hpp new file mode 100644 index 00000000..92ca25c7 --- /dev/null +++ b/opm/core/fluid/blackoil/BlackoilPVT.hpp @@ -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 . +*/ + +#ifndef OPM_BLACKOILPVT_HEADER_INCLUDED +#define OPM_BLACKOILPVT_HEADER_INCLUDED + + +#include +#include +#include +#include +#include + + +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& pressures, + const std::vector& surfvol, + std::vector& output) const; + void B(const std::vector& pressures, + const std::vector& surfvol, + std::vector& output) const; + void dBdp(const std::vector& pressures, + const std::vector& surfvol, + std::vector& output_B, + std::vector& output_dBdp) const; + void R(const std::vector& pressures, + const std::vector& surfvol, + std::vector& output) const; + void dRdp(const std::vector& pressures, + const std::vector& surfvol, + std::vector& output_R, + std::vector& output_dRdp) const; + + private: + int region_number_; + const MiscibilityProps& propsForPhase(PhaseIndex phase) const; + + boost::scoped_ptr water_props_; + boost::scoped_ptr oil_props_; + boost::scoped_ptr gas_props_; + CompVec densities_; + mutable std::vector data1_; + mutable std::vector data2_; + }; + +} + + +#endif // OPM_BLACKOILPVT_HEADER_INCLUDED diff --git a/opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp b/opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp new file mode 100644 index 00000000..b03f4700 --- /dev/null +++ b/opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp @@ -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 . +*/ + +#ifndef OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED +#define OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + +// Forward declaration needed by associated parameters class. +template +class FluidMatrixInteractionBlackoil; + +template +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& sw = swof_table[0][0]; + const std::vector& krw = swof_table[0][1]; + const std::vector& krow = swof_table[0][2]; + const std::vector& pcow = swof_table[0][3]; + const std::vector& sg = sgof_table[0][0]; + const std::vector& krg = sgof_table[0][1]; + const std::vector& krog = sgof_table[0][2]; + const std::vector& 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 + friend class FluidMatrixInteractionBlackoil; + + Dune::utils::UniformTableLinear krw_; + Dune::utils::UniformTableLinear krow_; + Dune::utils::UniformTableLinear pcow_; + Dune::utils::UniformTableLinear krg_; + Dune::utils::UniformTableLinear krog_; + Dune::utils::UniformTableLinear pcog_; + Scalar krocw_; // = krow_(s_wc) +}; + + +/*! + * \ingroup material + * + * \brief Capillary pressures and relative permeabilities for a black oil system. + */ +template > +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 + 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 + 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 + 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 + 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 diff --git a/opm/core/fluid/blackoil/FluidStateBlackoil.hpp b/opm/core/fluid/blackoil/FluidStateBlackoil.hpp new file mode 100644 index 00000000..fd96edb3 --- /dev/null +++ b/opm/core/fluid/blackoil/FluidStateBlackoil.hpp @@ -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 . +*/ + +#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 surface_volume_density; // z + std::vector phase_pressure; // p + + // Variables from PVT functions. + std::vector formation_volume_factor; // B + std::vector solution_factor; // R + std::vector viscosity; // mu + + // Variables computed from PVT data. + // The A matrices are all in Fortran order (or, equivalently, + // we store the transposes). + std::vector state_matrix; // A' = (RB^{-1})' + std::vector phase_volume_density; // u + std::vector total_phase_volume_density; // sum(u) + std::vector saturation; // s = u/sum(u) +#ifdef COMPUTE_OLD_TERMS + std::vector formation_volume_factor_deriv; // dB/dp + std::vector solution_factor_deriv; // dR/dp + std::vector phase_compressibility; // c + std::vector total_compressibility; // cT + std::vector experimental_term; // ex = sum(Ai*dA*Ai*z) +#endif + + // Variables computed from saturation. + std::vector relperm; // kr + std::vector relperm_deriv; // dkr/ds + std::vector mobility; // lambda + std::vector mobility_deriv; // dlambda/ds + }; + +} // end namespace Opm + + +#endif // OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED diff --git a/opm/core/fluid/blackoil/Makefile.am b/opm/core/fluid/blackoil/Makefile.am new file mode 100644 index 00000000..efec558f --- /dev/null +++ b/opm/core/fluid/blackoil/Makefile.am @@ -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 diff --git a/opm/core/fluid/blackoil/MiscibilityDead.cpp b/opm/core/fluid/blackoil/MiscibilityDead.cpp new file mode 100644 index 00000000..eeeddb46 --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityDead.cpp @@ -0,0 +1,179 @@ +//=========================================================================== +// +// File: MiscibilityDead.cpp +// +// Created: Wed Feb 10 09:06:04 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + +#include +#include +#include +#include +#include +#include +#include +#include + +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 press(sz); + std::vector B_inv(sz); + std::vector 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(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& pressures, + const std::vector&, + int phase, + std::vector& 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& pressures, + const std::vector&, + int phase, + std::vector& 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& pressures, + const std::vector& surfvols, + int phase, + std::vector& output_B, + std::vector& 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& pressures, + const std::vector&, + int, + std::vector& 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& pressures, + const std::vector&, + int, + std::vector& output_R, + std::vector& 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); + } + +} diff --git a/opm/core/fluid/blackoil/MiscibilityDead.hpp b/opm/core/fluid/blackoil/MiscibilityDead.hpp new file mode 100644 index 00000000..f6f1a96e --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityDead.hpp @@ -0,0 +1,89 @@ +//=========================================================================== +// +// File: MiscibilityDead.hpp +// +// Created: Wed Feb 10 09:05:47 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + +#ifndef SINTEF_MISCIBILITYDEAD_HEADER +#define SINTEF_MISCIBILITYDEAD_HEADER + +/** Class for immiscible dead oil and dry gas. + * Detailed description. + */ + +#include +#include + +namespace Opm +{ + class MiscibilityDead : public MiscibilityProps + { + public: + typedef std::vector > > 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void B(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void dBdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_B, + std::vector& output_dBdp) const; + virtual void R(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void dRdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_R, + std::vector& output_dRdp) const; + + private: + // PVT properties of dry gas or dead oil + Dune::utils::UniformTableLinear one_over_B_; + Dune::utils::UniformTableLinear viscosity_; + }; + +} + +#endif // SINTEF_MISCIBILITYDEAD_HEADER + diff --git a/opm/core/fluid/blackoil/MiscibilityLiveGas.cpp b/opm/core/fluid/blackoil/MiscibilityLiveGas.cpp new file mode 100644 index 00000000..6e101277 --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityLiveGas.cpp @@ -0,0 +1,286 @@ +//=========================================================================== +// +// File: MiscibilityLiveGas.cpp +// +// Created: Wed Feb 10 09:21:53 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + + +#include +#include +#include +#include + +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& pressures, + const std::vector& surfvol, + int phase, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_R, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_B, + std::vector& 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 diff --git a/opm/core/fluid/blackoil/MiscibilityLiveGas.hpp b/opm/core/fluid/blackoil/MiscibilityLiveGas.hpp new file mode 100644 index 00000000..801649ed --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityLiveGas.hpp @@ -0,0 +1,92 @@ +//=========================================================================== +// +// File: MiscibilityLiveGas.hpp +// +// Created: Wed Feb 10 09:21:26 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + +#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 > > 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void B(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void dBdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_B, + std::vector& output_dBdp) const; + virtual void R(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void dRdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_R, + std::vector& 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 > saturated_gas_table_; + std::vector > > undersat_gas_tables_; + + }; + +} + +#endif // SINTEF_MISCIBILITYLIVEGAS_HEADER + diff --git a/opm/core/fluid/blackoil/MiscibilityLiveOil.cpp b/opm/core/fluid/blackoil/MiscibilityLiveOil.cpp new file mode 100644 index 00000000..60a67a1f --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityLiveOil.cpp @@ -0,0 +1,395 @@ +//=========================================================================== +// +// File: MiscibiltyLiveOil.cpp +// +// Created: Wed Feb 10 09:08:25 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + +#include +#include +#include +#include + +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 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::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::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& pressures, + const std::vector& surfvol, + int phase, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_R, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_B, + std::vector& 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 diff --git a/opm/core/fluid/blackoil/MiscibilityLiveOil.hpp b/opm/core/fluid/blackoil/MiscibilityLiveOil.hpp new file mode 100644 index 00000000..6fd8ea42 --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityLiveOil.hpp @@ -0,0 +1,97 @@ +//=========================================================================== +// +// File: MiscibilityLiveOil.hpp +// +// Created: Wed Feb 10 09:08:09 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + +#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 > > 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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void B(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void dBdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_B, + std::vector& output_dBdp) const; + virtual void R(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const; + virtual void dRdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_R, + std::vector& 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 > saturated_oil_table_; + std::vector > > undersat_oil_tables_; + }; + +} + +#endif // SINTEF_MISCIBILITYLIVEOIL_HEADER + diff --git a/opm/core/fluid/blackoil/MiscibilityProps.cpp b/opm/core/fluid/blackoil/MiscibilityProps.cpp new file mode 100644 index 00000000..012606f7 --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityProps.cpp @@ -0,0 +1,52 @@ +//=========================================================================== +// +// File: MiscibilityProps.cpp +// +// Created: Wed Feb 10 09:05:05 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + +#include "MiscibilityProps.hpp" + +using namespace std; + +namespace Opm +{ + + + //------------------------------------------------------------------------ + // Member functions + //------------------------------------------------------------------------- + + /// Constructor + MiscibilityProps::MiscibilityProps() + { + } + + MiscibilityProps::~MiscibilityProps() + { + } + +} // namespace Opm diff --git a/opm/core/fluid/blackoil/MiscibilityProps.hpp b/opm/core/fluid/blackoil/MiscibilityProps.hpp new file mode 100644 index 00000000..2bbacf33 --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityProps.hpp @@ -0,0 +1,87 @@ +//=========================================================================== +// +// File: MiscibilityProps.hpp +// +// Created: Wed Feb 10 09:04:35 2010 +// +// Author: Bjørn Spjelkavik +// +// 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 . +*/ + +#ifndef SINTEF_MISCIBILITYPROPS_HEADER +#define SINTEF_MISCIBILITYPROPS_HEADER + + /** Base class for properties of fluids and rocks. + * Detailed description. + */ + +#include "BlackoilDefs.hpp" +#include +#include + +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& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const = 0; + virtual void B(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const = 0; + virtual void dBdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_B, + std::vector& output_dBdp) const = 0; + virtual void R(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output) const = 0; + virtual void dRdp(const std::vector& pressures, + const std::vector& surfvol, + int phase, + std::vector& output_R, + std::vector& output_dRdp) const = 0; + }; + +} // namespace Opm + +#endif // SINTEF_MISCIBILITYPROPS_HEADER + diff --git a/opm/core/fluid/blackoil/MiscibilityWater.hpp b/opm/core/fluid/blackoil/MiscibilityWater.hpp new file mode 100644 index 00000000..0075f059 --- /dev/null +++ b/opm/core/fluid/blackoil/MiscibilityWater.hpp @@ -0,0 +1,194 @@ +//=========================================================================== +// +// File: MiscibilityWater.hpp +// +// Created: Tue May 18 10:26:13 2010 +// +// Author(s): Atgeirr F Rasmussen +// Bjørn Spjelkavik +// +// $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 . +*/ + +#ifndef OPENRS_MISCIBILITYWATER_HEADER +#define OPENRS_MISCIBILITYWATER_HEADER + +#include +#include + +// Forward declaration. +class PVTW; + +namespace Opm +{ + class MiscibilityWater : public MiscibilityProps + { + public: + typedef std::vector > 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& pressures, + const std::vector&, + int, + std::vector& 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& pressures, + const std::vector&, + int phase, + std::vector& 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& pressures, + const std::vector& surfvols, + int phase, + std::vector& output_B, + std::vector& 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& pressures, + const std::vector&, + int, + std::vector& 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& pressures, + const std::vector&, + int, + std::vector& output_R, + std::vector& 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 diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp new file mode 100644 index 00000000..c8840770 --- /dev/null +++ b/opm/core/utility/UniformTableLinear.hpp @@ -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 . +*/ + +#ifndef OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED +#define OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED + +#include +#include +#include +#include +#include + +#include +#include + +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 + 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& 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 domain(); + + /// @brief Rescale the domain. + /// @param new_domain the new domain as a pair of doubles. + void rescaleDomain(std::pair 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 y_values_; + RangePolicy left_; + RangePolicy right_; + template + friend std::ostream& operator<<(std::ostream& os, const UniformTableLinear& t); + }; + + + // Member implementations. + + template + inline + UniformTableLinear + ::UniformTableLinear() + : left_(ClosestValue), right_(ClosestValue) + { + } + + template + inline + UniformTableLinear + ::UniformTableLinear(double xmin, + double xmax, + const std::vector& 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 + inline + UniformTableLinear + ::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 + inline std::pair + UniformTableLinear + ::domain() + { + return std::make_pair(xmin_, xmax_); + } + + template + inline void + UniformTableLinear + ::rescaleDomain(std::pair new_domain) + { + xmin_ = new_domain.first; + xmax_ = new_domain.second; + xdelta_ = (xmax_ - xmin_)/(y_values_.size() - 1); + } + + template + inline double + UniformTableLinear + ::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 + inline double + UniformTableLinear + ::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 + inline bool + UniformTableLinear + ::operator==(const UniformTableLinear& other) const + { + return xmin_ == other.xmin_ + && xdelta_ == other.xdelta_ + && y_values_ == other.y_values_ + && left_ == other.left_ + && right_ == other.right_; + } + + template + inline void + UniformTableLinear + ::setLeftPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + left_ = rp; + } + + template + inline void + UniformTableLinear + ::setRightPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + right_ = rp; + } + + + template + inline std::ostream& operator<<(std::ostream& os, const UniformTableLinear& 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 diff --git a/opm/core/utility/buildUniformMonotoneTable.hpp b/opm/core/utility/buildUniformMonotoneTable.hpp new file mode 100644 index 00000000..09501f9a --- /dev/null +++ b/opm/core/utility/buildUniformMonotoneTable.hpp @@ -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 . +*/ + +#ifndef OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED +#define OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED + +#include +#include + +namespace Dune { + namespace utils { + + template + void buildUniformMonotoneTable(const std::vector& xv, + const std::vector& yv, + const int samples, + UniformTableLinear& table) + { + MonotCubicInterpolator interp(xv, yv); + std::vector 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(xmin, xmax, uniform_yv); + } + + } // namespace utils +} // namespace Dune + + + +#endif // OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED diff --git a/opm/core/utility/linearInterpolation.hpp b/opm/core/utility/linearInterpolation.hpp new file mode 100644 index 00000000..830545f8 --- /dev/null +++ b/opm/core/utility/linearInterpolation.hpp @@ -0,0 +1,90 @@ +//=========================================================================== +// +// File: linearInterpolation.hpp +// +// Created: Tue Sep 9 12:49:39 2008 +// +// Author(s): Atgeirr F Rasmussen +// +// $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 . +*/ + +#ifndef OPENRS_LINEARINTERPOLATION_HEADER +#define OPENRS_LINEARINTERPOLATION_HEADER + + +#include +#include + +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 + T linearInterpolation(const std::vector& xv, + const std::vector& yv, + double x) + { + std::vector::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 + T linearInterpolationDerivative(const std::vector& xv, + const std::vector& 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 diff --git a/tests/Makefile.am b/tests/Makefile.am index 98a3822b..b8d95ec6 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -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 diff --git a/tests/bo_fluid_p_and_z_deps.cpp b/tests/bo_fluid_p_and_z_deps.cpp new file mode 100644 index 00000000..af3d69bc --- /dev/null +++ b/tests/bo_fluid_p_and_z_deps.cpp @@ -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 . +*/ + + +#include "config.h" + +#include +#include +#include + + +int main(int argc, char** argv) +{ + std::cout << "%{\n"; + + // Parameters. + Dune::parameter::ParameterGroup param(argc, argv); + + // Parser. + std::string ecl_file = param.get("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"; +} + diff --git a/tests/bo_fluid_pressuredeps.cpp b/tests/bo_fluid_pressuredeps.cpp new file mode 100644 index 00000000..7c277976 --- /dev/null +++ b/tests/bo_fluid_pressuredeps.cpp @@ -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 . +*/ + +#include "config.h" + +#include +#include +#include + +#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("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 + } +} diff --git a/tests/bo_fluid_test.cpp b/tests/bo_fluid_test.cpp new file mode 100644 index 00000000..019b58dc --- /dev/null +++ b/tests/bo_fluid_test.cpp @@ -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 . +*/ + +#include "config.h" + +#include +#include +#include + + +int main(int argc, char** argv) +{ + // Parameters. + Dune::parameter::ParameterGroup param(argc, argv); + + // Parser. + std::string ecl_file = param.get("filename"); + Dune::EclipseGridParser parser(ecl_file); + + // Test the FluidMatrixInteractionBlackoil class. + Opm::FluidMatrixInteractionBlackoilParams fluid_params; + fluid_params.init(parser); + typedef Opm::FluidMatrixInteractionBlackoil 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'; + } + +}