diff --git a/dune/porsol/blackoil/BlackoilFluid.hpp b/dune/porsol/blackoil/BlackoilFluid.hpp new file mode 100644 index 00000000..914ceb6b --- /dev/null +++ b/dune/porsol/blackoil/BlackoilFluid.hpp @@ -0,0 +1,68 @@ +/* + 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 +{ + + class BlackoilFluid + { + public: + enum { numPhases = 3 }; + enum { numComponents = 3 }; + typedef Dune::FieldVector PhaseVec; + typedef Dune::FieldVector CompVec; + + void init(const Dune::EclipseGridParser& parser) + { + fmi_params_.init(parser); + FluidSystemBlackoil<>::init(parser); + } + + FluidStateBlackoil computeState(PhaseVec phase_pressure, CompVec z) + { + FluidStateBlackoil state; + state.temperature_ = 300; + state.phase_pressure_ = phase_pressure; + state.surface_volume_ = z; + FluidSystemBlackoil<>::computeEquilibrium(state); // Sets everything but relperm and mobility. + FluidMatrixInteractionBlackoil::kr(state.relperm_, fmi_params_, state.saturation_, state.temperature_); + for (int phase = 0; phase < numPhases; ++phase) { + state.mobility_[phase] = state.relperm_[phase]/state.viscosity_[phase]; + } + return state; + } + + private: + FluidMatrixInteractionBlackoilParams fmi_params_; + }; + +} // namespace Opm + +#endif // OPM_BLACKOILFLUID_HEADER_INCLUDED diff --git a/dune/porsol/blackoil/fluid/BlackoilDefs.hpp b/dune/porsol/blackoil/fluid/BlackoilDefs.hpp index 3a38d8ae..c0cfe15b 100644 --- a/dune/porsol/blackoil/fluid/BlackoilDefs.hpp +++ b/dune/porsol/blackoil/fluid/BlackoilDefs.hpp @@ -20,6 +20,10 @@ #ifndef OPM_BLACKOILDEFS_HEADER_INCLUDED #define OPM_BLACKOILDEFS_HEADER_INCLUDED + +#include + + namespace Opm { @@ -31,6 +35,10 @@ namespace Opm enum ComponentIndex { Water = 0, Gas = 1, Oil = 2 }; enum PhaseIndex { Aqua = 0, Vapour = 1, Liquid = 2 }; + + typedef double Scalar; + typedef Dune::FieldVector CompVec; + typedef Dune::FieldVector PhaseVec; }; } // namespace Opm diff --git a/dune/porsol/blackoil/fluid/BlackoilPVT.cpp b/dune/porsol/blackoil/fluid/BlackoilPVT.cpp index e568534c..7ae03c6d 100644 --- a/dune/porsol/blackoil/fluid/BlackoilPVT.cpp +++ b/dune/porsol/blackoil/fluid/BlackoilPVT.cpp @@ -78,32 +78,32 @@ namespace Opm } } - BlackoilPVT::surfvol_t BlackoilPVT::surfaceDensities() const + BlackoilPVT::CompVec BlackoilPVT::surfaceDensities() const { return densities_; } - double BlackoilPVT::getViscosity(double press, const surfvol_t& surfvol, PhaseIndex phase) const + 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 surfvol_t& surfvol, PhaseIndex phase) const + 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 surfvol_t& surfvol, PhaseIndex phase) const + 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 surfvol_t& surfvol, PhaseIndex phase) const + 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 surfvol_t& surfvol, PhaseIndex phase) const + double BlackoilPVT::dRdp(double press, const CompVec& surfvol, PhaseIndex phase) const { return propsForPhase(phase).dRdp(region_number_, press, surfvol); } diff --git a/dune/porsol/blackoil/fluid/BlackoilPVT.hpp b/dune/porsol/blackoil/fluid/BlackoilPVT.hpp index 1f5156d4..0563daab 100644 --- a/dune/porsol/blackoil/fluid/BlackoilPVT.hpp +++ b/dune/porsol/blackoil/fluid/BlackoilPVT.hpp @@ -33,20 +33,23 @@ namespace Opm class BlackoilPVT : public BlackoilDefs { public: - typedef MiscibilityProps::surfvol_t surfvol_t; - void init(const Dune::EclipseGridParser& ep); - double getViscosity(double press, const surfvol_t& surfvol, + double getViscosity(double press, + const CompVec& surfvol, PhaseIndex phase) const; - surfvol_t surfaceDensities() const; - double B (double press, const surfvol_t& surfvol, + CompVec surfaceDensities() const; + double B (double press, + const CompVec& surfvol, PhaseIndex phase) const; - double dBdp(double press, const surfvol_t& surfvol, + double dBdp(double press, + const CompVec& surfvol, PhaseIndex phase) const; - double R (double press, const surfvol_t& surfvol, + double R (double press, + const CompVec& surfvol, PhaseIndex phase) const; - double dRdp(double press, const surfvol_t& surfvol, + double dRdp(double press, + const CompVec& surfvol, PhaseIndex phase) const; @@ -57,7 +60,7 @@ namespace Opm boost::scoped_ptr water_props_; boost::scoped_ptr oil_props_; boost::scoped_ptr gas_props_; - surfvol_t densities_; + CompVec densities_; }; } diff --git a/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp b/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp index 79f296de..da8e025e 100644 --- a/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp +++ b/dune/porsol/blackoil/fluid/FluidStateBlackoil.hpp @@ -21,53 +21,25 @@ #define OPM_FLUIDSTATEBLACKOIL_HEADER_INCLUDED -#include "FluidSystemBlackoil.hpp" +#include "BlackoilDefs.hpp" namespace Opm { /*! - * \brief Calculates the phase state from the primary variables in the - * blackoil model. + * \brief Fluid states for a black oil model. */ struct FluidStateBlackoil : public BlackoilDefs { - typedef FluidSystemBlackoil FluidSystem; - typedef double Scalar; - typedef std::tr1::array PrimaryVariables; // Surface volumes and phase pressures. - typedef FluidMatrixInteractionBlackoil MaterialLaw; - typedef typename MaterialLaw::Params MaterialLawParams; - Scalar temperature_; - Scalar surface_volume_[numComponents]; - Scalar phase_pressure_[numPhases]; - Scalar phase_volume_[numPhases]; - Scalar saturation_[numPhases]; - - - /*! - * \brief Update the phase state from the primary variables. - */ - void update(const PrimaryVariables& primary_vars, - const MaterialLawParams& material_params, - Scalar temperature) - { - // Set the temperature. - temperature_ = temperature; - - // Set the surface volumes. - for (int i = 0; i < numComponents; ++i) { - surface_volume_[i] = primary_vars[i]; - } - - // Set the phase pressures. - for (int i = 0; i < numPhases; ++i) { - phase_pressure_[i] = primary_vars[numComponents + i]; - } - - // Compute phase volumes by the fluid system rules. - FluidSystem::computeEquilibrium(*this); - } + CompVec surface_volume_; + PhaseVec phase_pressure_; + PhaseVec phase_volume_; + Scalar phase_to_comp_[numPhases*numComponents]; + PhaseVec saturation_; + PhaseVec viscosity_; + PhaseVec relperm_; + PhaseVec mobility_; }; } // end namespace Opm diff --git a/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp b/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp index 6077bc1c..ccbcca18 100644 --- a/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp +++ b/dune/porsol/blackoil/fluid/FluidSystemBlackoil.hpp @@ -35,6 +35,8 @@ #define OPM_FLUIDSYSTEMBLACKOIL_HEADER_INCLUDED #include "BlackoilPVT.hpp" +#include "BlackoilDefs.hpp" +#include #include #include @@ -64,19 +66,10 @@ private: * \brief A black oil fluid system. */ template -class FluidSystemBlackoil +class FluidSystemBlackoil : public BlackoilDefs { public: typedef ParamsT Params; - typedef double Scalar; - - enum { numComponents = 3 }; - enum { numPhases = 3 }; - - enum ComponentIndex { Water = 0, Gas = 1, Oil = 2 }; - enum PhaseIndex { Aqua = 0, Vapour = 1, Liquid = 2 }; - - typedef BlackoilPVT::surfvol_t FluidVec; /*! * \brief Initialize system from input. @@ -151,16 +144,16 @@ public: /*! - * \brief Assuming the surface volumes and the pressure of all - * phases are known, compute the phase volumes and - * saturations. + * \brief Assuming the surface volumes and the pressures of all + * phases are known, compute everything except relperm and + * mobility. */ template static void computeEquilibrium(FluidState& fluid_state) { // Get B and R factors. - const double* p = fluid_state.phase_pressure_; - const double* z = fluid_state.surface_volume_; + const PhaseVec& p = fluid_state.phase_pressure_; + const CompVec& z = fluid_state.surface_volume_; double B[3]; B[Aqua] = params().pvt_.B(p[Aqua], z, Aqua); B[Vapour] = params().pvt_.B(p[Vapour], z, Vapour); @@ -168,9 +161,17 @@ public: double R[3]; // Only using 2 of them, though. R[Vapour] = params().pvt_.B(p[Vapour], z, Vapour); R[Liquid] = params().pvt_.B(p[Liquid], z, Liquid); + // Set the A matrix (A = RB^{-1}) + Dune::SharedFortranMatrix A(numComponents, numPhases, fluid_state.phase_to_comp_); + zero(A); + A(Water, Aqua) = 1.0/B[Aqua]; + A(Gas, Vapour) = 1.0/B[Vapour]; + A(Gas, Liquid) = R[Liquid]/B[Liquid]; + A(Oil, Vapour) = R[Vapour]/B[Vapour]; + A(Oil, Liquid) = 1.0/B[Liquid]; - // Update phase volumes. - double* u = fluid_state.phase_volume_; + // Update phase volumes. This is the same as multiplying with A^{-1} + PhaseVec& u = fluid_state.phase_volume_; 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; @@ -178,10 +179,16 @@ public: // Update saturations. double sumu = u[Aqua] + u[Vapour] + u[Liquid]; - double* s = fluid_state.saturation_; + PhaseVec& s = fluid_state.saturation_; for (int i = 0; i < 3; ++i) { s[i] = u[i]/sumu; } + + // Compute viscosities. + PhaseVec& mu = fluid_state.viscosity_; + mu[Aqua] = params().pvt_.getViscosity(p[Aqua], z, Aqua); + mu[Vapour] = params().pvt_.getViscosity(p[Vapour], z, Vapour); + mu[Liquid] = params().pvt_.getViscosity(p[Liquid], z, Liquid); } /*! diff --git a/dune/porsol/blackoil/fluid/MiscibilityProps.hpp b/dune/porsol/blackoil/fluid/MiscibilityProps.hpp index 588c5e0f..693f87ec 100644 --- a/dune/porsol/blackoil/fluid/MiscibilityProps.hpp +++ b/dune/porsol/blackoil/fluid/MiscibilityProps.hpp @@ -46,7 +46,7 @@ namespace Opm class MiscibilityProps : public BlackoilDefs { public: - typedef std::tr1::array surfvol_t; + typedef CompVec surfvol_t; MiscibilityProps(); virtual ~MiscibilityProps();