diff --git a/opm/core/fluid/BlackoilFluid.hpp b/opm/core/fluid/BlackoilFluid.hpp index d1cd6df2..6aa3142d 100644 --- a/opm/core/fluid/BlackoilFluid.hpp +++ b/opm/core/fluid/BlackoilFluid.hpp @@ -21,11 +21,10 @@ #define OPM_BLACKOILFLUID_HEADER_INCLUDED -#include -#include -#include -#include -#include +#include +#include +#include +#include #include @@ -174,26 +173,34 @@ namespace Opm 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& dB = states.formation_volume_factor_deriv[i]; const PhaseVec& R = states.solution_factor[i]; - const PhaseVec& dR = states.solution_factor_deriv[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 } } @@ -276,6 +283,13 @@ namespace Opm 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); @@ -284,19 +298,17 @@ namespace Opm dR[Aqua] = 0.0; dR[Vapour] = pvt_.dRdp(p[Vapour], z, Vapour); dR[Liquid] = pvt_.dRdp(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_; 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_; @@ -307,6 +319,7 @@ namespace Opm +#ifdef COMPUTE_OLD_TERMS static void computeSingleEquilibrium(const PhaseVec& B, const PhaseVec& dB, const PhaseVec& R, @@ -319,6 +332,15 @@ namespace Opm 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 @@ -345,6 +367,7 @@ namespace Opm 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}) @@ -372,6 +395,7 @@ namespace Opm dAt.mtv(tmp1, tmp2); Ait.mtv(tmp2, tmp3); exp_term = tmp3[Aqua] + tmp3[Liquid] + tmp3[Gas]; +#endif } diff --git a/opm/core/fluid/blackoil/BlackoilComponent.hpp b/opm/core/fluid/blackoil/BlackoilComponent.hpp deleted file mode 100644 index 61294a68..00000000 --- a/opm/core/fluid/blackoil/BlackoilComponent.hpp +++ /dev/null @@ -1,165 +0,0 @@ -/* - Copyright 2010 SINTEF ICT, Applied Mathematics. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_BLACKOILCOMPONENT_HEADER_INCLUDED -#define OPM_BLACKOILCOMPONENT_HEADER_INCLUDED - - -#include - -namespace Dumux -{ - -/*! - * \ingroup Components - - * \brief - * A component class for the black oil model, intended to be used - * for all three components. - * - * \tparam Scalar The type used for scalar values - */ -template -class BlackoilComponent -{ -public: - /*! - * \brief A human readable name for the component. - */ - static const char *name() - { - return "BlackoilComponent"; - } - - /*! - * \brief The molar mass in [kg] of the component. - */ - static Scalar molarMass() - { DUNE_THROW(Dune::NotImplemented, "Component::molarMass()"); } - - /*! - * \brief Returns the critical temperature in [K] of the component. - */ - static Scalar criticalTemperature() - { DUNE_THROW(Dune::NotImplemented, "Component::criticalTemperature()"); } - - /*! - * \brief Returns the critical pressure in [Pa] of the component. - */ - static Scalar criticalPressure() - { DUNE_THROW(Dune::NotImplemented, "Component::criticalPressure()"); } - - /*! - * \brief Returns the temperature in [K] at the component's triple point. - */ - static Scalar tripleTemperature() - { DUNE_THROW(Dune::NotImplemented, "Component::tripleTemperature()"); } - - /*! - * \brief Returns the pressure in [Pa] at the component's triple point. - */ - static Scalar triplePressure() - { DUNE_THROW(Dune::NotImplemented, "Component::triplePressure()"); } - - /*! - * \brief The vapor pressure in [Pa] of the component at a given - * temperature in [K]. - * - * \param T temperature of the component in [K] - */ - static Scalar vaporPressure(Scalar T) - { DUNE_THROW(Dune::NotImplemented, "Component::vaporPressure()"); } - - /*! - * \brief The density in [kg/m^3] of the component at a given pressure in [Pa] and temperature in [K]. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static Scalar gasDensity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::density()"); } - - /*! - * \brief The density [kg/m^3] of the liquid component at a given pressure in [Pa] and temperature in [K]. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static Scalar liquidDensity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::density()"); } - - /*! - * \brief Specific enthalpy [J/kg] of the pure component in gas. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static const Scalar gasEnthalpy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasEnthalpy()"); } - - /*! - * \brief Specific enthalpy [J/kg] of the pure component in liquid. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static const Scalar liquidEnthalpy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidEnthalpy()"); } - - /*! - * \brief Specific internal energy [J/kg] of the pure component in gas. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static const Scalar gasInternalEnergy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasInternalEnergy()"); } - - /*! - * \brief Specific internal energy [J/kg] of pure the pure component in liquid. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static const Scalar liquidInternalEnergy(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidInternalEnergy()"); } - - /*! - * \brief The dynamic viscosity [Pa*s] of the pure component at a given pressure in [Pa] and temperature in [K]. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static Scalar gasViscosity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::gasViscosity()"); } - - /*! - * \brief The dynamic liquid viscosity [Pa*s] of the pure component. - * - * \param temperature temperature of component in [K] - * \param pressure pressure of component in [Pa] - */ - static Scalar liquidViscosity(Scalar temperature, Scalar pressure) - { DUNE_THROW(Dune::NotImplemented, "Component::liquidViscosity()"); } - -}; - -} // end namepace - -#endif // OPM_BLACKOILCOMPONENT_HEADER_INCLUDED diff --git a/opm/core/fluid/blackoil/BlackoilDefs.hpp b/opm/core/fluid/blackoil/BlackoilDefs.hpp index 7615c5e6..fba7546c 100644 --- a/opm/core/fluid/blackoil/BlackoilDefs.hpp +++ b/opm/core/fluid/blackoil/BlackoilDefs.hpp @@ -22,6 +22,7 @@ #include +#include #include namespace Opm @@ -36,16 +37,63 @@ namespace Opm 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) + { data_.assign(elem); } // In C++11, assign -> fill + SmallVec& operator=(const T& elem) + { data_.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: + std::tr1::array data_; + }; + template + class SmallMat + { + public: + SmallMat() + {} + SmallMat(const T& elem) + { data_.assign(elem); } // In C++11, assign -> fill + 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 Dune::FieldVector CompVec; -// typedef Dune::FieldVector PhaseVec; - typedef std::tr1::array CompVec; - typedef std::tr1::array PhaseVec; + typedef SmallVec CompVec; + typedef SmallVec PhaseVec; BOOST_STATIC_ASSERT(int(numComponents) == int(numPhases)); -// typedef Dune::FieldMatrix PhaseToCompMatrix; -// typedef Dune::FieldMatrix PhaseJacobian; - typedef std::tr1::array PhaseToCompMatrix; - typedef std::tr1::array PhaseJacobian; + 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 diff --git a/opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp b/opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp index cd323f06..b03f4700 100644 --- a/opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp +++ b/opm/core/fluid/blackoil/FluidMatrixInteractionBlackoil.hpp @@ -20,10 +20,10 @@ #ifndef OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED #define OPM_FLUIDMATRIXINTERACTIONBLACKOIL_HEADER_INCLUDED -#include -#include -#include -#include "BlackoilDefs.hpp" +#include +#include +#include +#include #include #include diff --git a/opm/core/fluid/blackoil/FluidStateBlackoil.hpp b/opm/core/fluid/blackoil/FluidStateBlackoil.hpp index dfc8a17b..fd96edb3 100644 --- a/opm/core/fluid/blackoil/FluidStateBlackoil.hpp +++ b/opm/core/fluid/blackoil/FluidStateBlackoil.hpp @@ -40,9 +40,11 @@ namespace Opm 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_; @@ -62,9 +64,7 @@ namespace Opm // Variables from PVT functions. std::vector formation_volume_factor; // B - std::vector formation_volume_factor_deriv; // dB/dp std::vector solution_factor; // R - std::vector solution_factor_deriv; // dR/dp std::vector viscosity; // mu // Variables computed from PVT data. @@ -74,9 +74,13 @@ namespace Opm 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 diff --git a/opm/core/fluid/blackoil/MiscibilityDead.cpp b/opm/core/fluid/blackoil/MiscibilityDead.cpp index 2a490ef2..eeeddb46 100644 --- a/opm/core/fluid/blackoil/MiscibilityDead.cpp +++ b/opm/core/fluid/blackoil/MiscibilityDead.cpp @@ -82,7 +82,7 @@ namespace Opm { } - double MiscibilityDead::getViscosity(int region, double press, const surfvol_t& /*surfvol*/) const + double MiscibilityDead::getViscosity(int /*region*/, double press, const surfvol_t& /*surfvol*/) const { return viscosity_(press); } @@ -100,7 +100,7 @@ namespace Opm } } - double MiscibilityDead::B(int region, double press, const surfvol_t& /*surfvol*/) const + double MiscibilityDead::B(int /*region*/, double press, const surfvol_t& /*surfvol*/) const { // Interpolate 1/B return 1.0/one_over_B_(press); diff --git a/tests/Makefile.am b/tests/Makefile.am index 98a3822b..d7e5f9ce 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -15,6 +15,9 @@ $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) noinst_PROGRAMS = \ +bo_fluid_p_and_z_deps \ +bo_fluid_pressuredeps \ +bo_fluid_test \ monotcubicinterpolator_test \ param_test \ sparsetable_test \ @@ -23,6 +26,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 index 3d1a5aaa..af3d69bc 100644 --- a/tests/bo_fluid_p_and_z_deps.cpp +++ b/tests/bo_fluid_p_and_z_deps.cpp @@ -20,9 +20,9 @@ #include "config.h" -#include -#include -#include +#include +#include +#include int main(int argc, char** argv) @@ -49,7 +49,11 @@ int main(int argc, char** argv) 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"; @@ -67,12 +71,14 @@ int main(int argc, char** argv) 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; diff --git a/tests/bo_fluid_pressuredeps.cpp b/tests/bo_fluid_pressuredeps.cpp index 1c6d1a44..7c277976 100644 --- a/tests/bo_fluid_pressuredeps.cpp +++ b/tests/bo_fluid_pressuredeps.cpp @@ -19,10 +19,51 @@ #include "config.h" -#include -#include -#include +#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) { @@ -43,7 +84,9 @@ int main(int argc, char** argv) 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; @@ -52,40 +95,14 @@ int main(int argc, char** argv) for (int phase = 0; phase < Opm::BlackoilFluid::numPhases; ++phase) { totmob += state.mobility_[phase]; } +#ifdef COMPUTE_OLD_TERMS if (print_compr) { - 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'; + printCompressibilityTerms(p, state); } else { - 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'; + printSatsEtc(p, state); } +#else + printSatsEtc(p, state); +#endif } } diff --git a/tests/bo_fluid_test.cpp b/tests/bo_fluid_test.cpp index d42a1660..019b58dc 100644 --- a/tests/bo_fluid_test.cpp +++ b/tests/bo_fluid_test.cpp @@ -19,9 +19,9 @@ #include "config.h" -#include -#include -#include +#include +#include +#include int main(int argc, char** argv) @@ -37,7 +37,7 @@ int main(int argc, char** argv) Opm::FluidMatrixInteractionBlackoilParams fluid_params; fluid_params.init(parser); typedef Opm::FluidMatrixInteractionBlackoil Law; - Dune::FieldVector s, kr; + Law::PhaseVec s, kr; const double temp = 300; // [K] int num = 41; for (int i = 0; i < num; ++i) { @@ -47,7 +47,7 @@ int main(int argc, char** argv) Law::kr(kr, fluid_params, s, temp); std::cout.width(6); std::cout.fill(' '); - std::cout << s[Law::Liquid] << " " << kr << '\n'; + std::cout << s[Law::Liquid] << " " << kr[0] << ' ' << kr[1] << ' ' << kr[2] << '\n'; } }