diff --git a/opm/material/components/C1.hpp b/opm/material/components/C1.hpp index 9344bf193..bb3c2ddc6 100644 --- a/opm/material/components/C1.hpp +++ b/opm/material/components/C1.hpp @@ -75,9 +75,9 @@ public: { return 4.60e6; /* [N/m^2] */ } /*! - * \brief Critical volume of \f$C_1\f$ [m2/kmol]. + * \brief Critical volume of \f$C_1\f$ [m3/kmol]. */ - static Scalar criticalVolume() {return 9.863e-5; } + static Scalar criticalVolume() {return 9.863e-2; } /*! * \brief Acentric factor of \f$C_1\f$. diff --git a/opm/material/components/C10.hpp b/opm/material/components/C10.hpp index c0230eb8f..b135b33c2 100644 --- a/opm/material/components/C10.hpp +++ b/opm/material/components/C10.hpp @@ -75,9 +75,9 @@ public: { return 2.10e6; /* [N/m^2] */ } /*! - * \brief Critical volume of \f$C_10\f$ [m2/kmol]. + * \brief Critical volume of \f$C_10\f$ [m3/kmol]. */ - static Scalar criticalVolume() {return 6.098e-4; } + static Scalar criticalVolume() {return 6.098e-1; } /*! * \brief Acentric factor of \f$C_10\f$. diff --git a/opm/material/components/SimpleCO2.hpp b/opm/material/components/SimpleCO2.hpp index f451a03bd..3d2168561 100644 --- a/opm/material/components/SimpleCO2.hpp +++ b/opm/material/components/SimpleCO2.hpp @@ -89,7 +89,7 @@ public: * \brief Critical volume of \f$CO_2\f$ [m2/kmol]. */ // Critical volume [m3/kmol] - static Scalar criticalVolume() {return 9.412e-5; } + static Scalar criticalVolume() {return 9.412e-2; } /*! * \brief Returns the pressure \f$\mathrm{[Pa]}\f$ at the triple point of \f$CO_2\f$. diff --git a/opm/material/constraintsolvers/PTFlash.hpp b/opm/material/constraintsolvers/PTFlash.hpp index 7dad59458..e8442584f 100644 --- a/opm/material/constraintsolvers/PTFlash.hpp +++ b/opm/material/constraintsolvers/PTFlash.hpp @@ -79,7 +79,6 @@ public: template static void solve(FluidState& fluid_state, const Dune::FieldVector& z, - int spatialIdx, std::string twoPhaseMethod, Scalar tolerance = -1., int verbosity = 0) @@ -104,7 +103,6 @@ public: // Print header if (verbosity >= 1) { std::cout << "********" << std::endl; - std::cout << "Flash calculations on Cell " << spatialIdx << std::endl; std::cout << "Inputs are K = [" << K << "], L = [" << L << "], z = [" << z << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl; } @@ -1127,7 +1125,7 @@ protected: const bool newton_afterwards, const int verbosity) { // Determine max. iterations based on if it will be used as a standalone flash or as a pre-process to Newton (or other) method. - const int maxIterations = newton_afterwards ? 3 : 10; + const int maxIterations = newton_afterwards ? 3 : 100; // Store cout format before manipulation std::ios_base::fmtflags f(std::cout.flags()); @@ -1169,7 +1167,7 @@ protected: } // Print iteration info - if (verbosity == 2 || verbosity == 4) { + if (verbosity >= 2) { int prec = 5; int fugWidth = (prec + 3); int convWidth = prec + 9; @@ -1184,7 +1182,7 @@ protected: } // Check convergence - if (convFugRatio.two_norm() < 1e-6){ + if (convFugRatio.two_norm() < 1e-6) { // Restore cout format std::cout.flags(f); diff --git a/opm/material/eos/PengRobinson.hpp b/opm/material/eos/PengRobinson.hpp index d3ea4e1a1..8a5bb659b 100644 --- a/opm/material/eos/PengRobinson.hpp +++ b/opm/material/eos/PengRobinson.hpp @@ -54,7 +54,7 @@ namespace Opm { * R. Reid, et al.: The Properties of Gases and Liquids, * 4th edition, McGraw-Hill, 1987, pp. 42-44, 82 */ -template +template class PengRobinson { //! The ideal gas constant [Pa * m^3/mol/K] @@ -205,27 +205,25 @@ public: Evaluation VmCubic = max(1e-7, Z[0]*RT/p); Vm = VmCubic; - // find the extrema (if they are present) - Evaluation Vmin, Vmax, pmin, pmax; - if (findExtrema_(Vmin, Vmax, - pmin, pmax, - a, b, T)) - { - if (isGasPhase) - Vm = std::max(Vmax, VmCubic); - else { - if (Vmin > 0) - Vm = std::min(Vmin, VmCubic); - else - Vm = VmCubic; + if (UseLegacy) { + // find the extrema (if they are present) + Evaluation Vmin, Vmax, pmin, pmax; + if (findExtrema_(Vmin, Vmax, pmin, pmax, a, b, T)) { + if (isGasPhase) + Vm = std::max(Vmax, VmCubic); + else { + if (Vmin > 0) + Vm = std::min(Vmin, VmCubic); + else + Vm = VmCubic; + } + } else { + // the EOS does not exhibit any physically meaningful + // extrema, and the fluid is critical... + Vm = VmCubic; + handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase); } } - else { - // the EOS does not exhibit any physically meaningful - // extrema, and the fluid is critical... - Vm = VmCubic; - handleCriticalFluid_(Vm, fs, params, phaseIdx, isGasPhase); - } } Valgrind::CheckDefined(Vm); diff --git a/opm/material/eos/PengRobinsonMixture.hpp b/opm/material/eos/PengRobinsonMixture.hpp index 2149e6935..5b8668f49 100644 --- a/opm/material/eos/PengRobinsonMixture.hpp +++ b/opm/material/eos/PengRobinsonMixture.hpp @@ -40,7 +40,6 @@ template class PengRobinsonMixture { enum { numComponents = StaticParameters::numComponents }; - typedef ::Opm::PengRobinson PengRobinson; // this class cannot be instantiated! PengRobinsonMixture() = default; @@ -53,20 +52,6 @@ class PengRobinsonMixture static const Scalar w; public: - /*! - * \brief Computes molar volumes where the Peng-Robinson EOS is - * true. - * - * \return Number of solutions. - */ - template - static int computeMolarVolumes(Scalar* Vm, - const MutableParams& params, - unsigned phaseIdx, - const FluidState& fs) - { - return PengRobinson::computeMolarVolumes(Vm, params, phaseIdx, fs); - } /*! * \brief Returns the fugacity coefficient of an individual diff --git a/opm/material/fluidsystems/Co2BrineFluidSystem.hh b/opm/material/fluidsystems/Co2BrineFluidSystem.hh index 1e432686b..397652366 100644 --- a/opm/material/fluidsystems/Co2BrineFluidSystem.hh +++ b/opm/material/fluidsystems/Co2BrineFluidSystem.hh @@ -34,13 +34,12 @@ #include #include -//#include namespace Opm { /*! * \ingroup FluidSystem * - * \brief A two phase two component system, co2 brine + * \brief A two phase two component system with components co2 brine */ template @@ -64,7 +63,6 @@ namespace Opm { template using ParameterCache = Opm::PTFlashParameterCache>; using ViscosityModel = typename Opm::ViscosityModels>; - //using ViscosityModel = typename Opm::ViscosityModels>; using PengRobinsonMixture = typename Opm::PengRobinsonMixture>; diff --git a/opm/material/fluidsystems/PTFlashParameterCache.hpp b/opm/material/fluidsystems/PTFlashParameterCache.hpp index 5056df397..83be5ab0e 100644 --- a/opm/material/fluidsystems/PTFlashParameterCache.hpp +++ b/opm/material/fluidsystems/PTFlashParameterCache.hpp @@ -48,7 +48,7 @@ class PTFlashParameterCache { using ThisType = PTFlashParameterCache; using ParentType = Opm::ParameterCacheBase; - using PengRobinson = Opm::PengRobinson; + using PengRobinson = Opm::PengRobinson; // false refer to discard some old code enum { numPhases = FluidSystem::numPhases }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; diff --git a/opm/material/fluidsystems/ThreeComponentFluidSystem.hh b/opm/material/fluidsystems/ThreeComponentFluidSystem.hh index ba70f6949..1726ff493 100644 --- a/opm/material/fluidsystems/ThreeComponentFluidSystem.hh +++ b/opm/material/fluidsystems/ThreeComponentFluidSystem.hh @@ -153,7 +153,7 @@ namespace Opm { static const char* name[] = {"o", // oleic phase "g"}; // gas phase - assert(0 <= phaseIdx && phaseIdx < 2); + assert(phaseIdx < 2); return name[phaseIdx]; } @@ -166,7 +166,7 @@ namespace Opm { Comp2::name(), }; - assert(0 <= compIdx && compIdx < 3); + assert(compIdx < 3); return name[compIdx]; } @@ -196,6 +196,7 @@ namespace Opm { // Use LBC method to calculate viscosity LhsEval mu; mu = ViscosityModel::LBC(fluidState, paramCache, phaseIdx); + return mu; } diff --git a/opm/material/viscositymodels/LBC.hpp b/opm/material/viscositymodels/LBC.hpp index 9a94426d9..567e7368c 100644 --- a/opm/material/viscositymodels/LBC.hpp +++ b/opm/material/viscositymodels/LBC.hpp @@ -49,54 +49,51 @@ public: unsigned phaseIdx) { const Scalar MPa_atm = 0.101325; - const Scalar R = 8.3144598e-3;//Mj/kmol*K + const Scalar R = Opm::Constants::R; const auto& T = Opm::decay(fluidState.temperature(phaseIdx)); - const auto& rho = Opm::decay(fluidState.density(phaseIdx)); + const auto& P = Opm::decay(fluidState.pressure(phaseIdx)); + const auto& Z = Opm::decay(fluidState.compressFactor(phaseIdx)); - LhsEval sumMm = 0.0; LhsEval sumVolume = 0.0; for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { - const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa; - const Scalar& T_c = FluidSystem::criticalTemperature(compIdx); - const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol; const auto& x = Opm::decay(fluidState.moleFraction(phaseIdx, compIdx)); - const Scalar v_c = FluidSystem::criticalVolume(compIdx); // in m3/kmol - sumMm += x*Mm; + const Scalar v_c = FluidSystem::criticalVolume(compIdx) / 1000; // converting to m3/mol from m3/kmol sumVolume += x*v_c; } - LhsEval rho_pc = sumMm/sumVolume; //mixture pseudocritical density - LhsEval rho_r = rho/rho_pc; + LhsEval rho_pc = 1.0 / sumVolume; + LhsEval V = (R * T * Z)/P; + LhsEval rho = 1.0 / V; + LhsEval rho_r = rho / rho_pc; - - LhsEval xsum_T_c = 0.0; //mixture pseudocritical temperature - LhsEval xsum_Mm = 0.0; //mixture molar mass - LhsEval xsum_p_ca = 0.0; //mixture pseudocritical pressure + LhsEval xsum_T_c = 0.0; // mixture pseudocritical temperature + LhsEval xsum_Mm = 0.0; // mixture molar mass + LhsEval xsum_p_ca = 0.0; // mixture pseudocritical pressure for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { - const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa; + const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6; // converting to Mpa from pascal const Scalar& T_c = FluidSystem::criticalTemperature(compIdx); - const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol; + const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; // converting to kg/kmol from kg/mol; const auto& x = Opm::decay(fluidState.moleFraction(phaseIdx, compIdx)); Scalar p_ca = p_c / MPa_atm; - xsum_T_c += x*T_c; - xsum_Mm += x*Mm; - xsum_p_ca += x*p_ca; + xsum_T_c += x * T_c; + xsum_Mm += x * Mm; + xsum_p_ca += x * p_ca; } LhsEval zeta_tot = Opm::pow(xsum_T_c / (Opm::pow(xsum_Mm,3.0) * Opm::pow(xsum_p_ca,4.0)),1./6); LhsEval my0 = 0.0; LhsEval sumxrM = 0.0; for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) { - const Scalar& p_c = FluidSystem::criticalPressure(compIdx)/1e6; // in Mpa; + const Scalar& p_c = FluidSystem::criticalPressure(compIdx) / 1e6; // converting to Mpa from pa; const Scalar& T_c = FluidSystem::criticalTemperature(compIdx); - const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; //in kg/kmol; + const Scalar Mm = FluidSystem::molarMass(compIdx) * 1000; // converting to kg/kmol from kg/mol; const auto& x = Opm::decay(fluidState.moleFraction(phaseIdx, compIdx)); Scalar p_ca = p_c / MPa_atm; Scalar zeta = std::pow(T_c / (std::pow(Mm,3.0) * std::pow(p_ca,4.0)),1./6); LhsEval T_r = T/T_c; LhsEval xrM = x * std::pow(Mm,0.5); LhsEval mys = 0.0; - if (T_r <=1.5) { + if (T_r <= 1.5) { mys = 34.0e-5*Opm::pow(T_r,0.94)/zeta; } else { mys = 17.78e-5*Opm::pow(4.58*T_r - 1.67, 0.625)/zeta; diff --git a/tests/material/test_co2brine_ptflash.cpp b/tests/material/test_co2brine_ptflash.cpp index 7347e896e..2b51f9d6f 100644 --- a/tests/material/test_co2brine_ptflash.cpp +++ b/tests/material/test_co2brine_ptflash.cpp @@ -127,9 +127,8 @@ BOOST_DATA_TEST_CASE(PtFlash, test_methods) const Evaluation Ltmp = 1.; fluid_state.setLvalue(Ltmp); - const int spatialIdx = 0; using Flash = Opm::PTFlash; - Flash::solve(fluid_state, z, spatialIdx, sample, flash_tolerance, flash_verbosity); + Flash::solve(fluid_state, z, sample, flash_tolerance, flash_verbosity); ComponentVector x, y; const Evaluation L = fluid_state.L(); @@ -241,9 +240,8 @@ BOOST_DATA_TEST_CASE(PtFlashSingle, test_methods) const Evaluation Ltmp = 1.; fluid_state.setLvalue(Ltmp); - const int spatialIdx = 0; using Flash = Opm::PTFlash; - Flash::solve(fluid_state, z, spatialIdx, sample, flash_tolerance, flash_verbosity); + Flash::solve(fluid_state, z, sample, flash_tolerance, flash_verbosity); ComponentVector x, y; const Evaluation L = fluid_state.L(); diff --git a/tests/material/test_threecomponents_ptflash.cpp b/tests/material/test_threecomponents_ptflash.cpp index 8470162c7..8f1f3a6f8 100644 --- a/tests/material/test_threecomponents_ptflash.cpp +++ b/tests/material/test_threecomponents_ptflash.cpp @@ -127,9 +127,8 @@ BOOST_DATA_TEST_CASE(PtFlash, test_methods) const Evaluation Ltmp = 1.; fluid_state.setLvalue(Ltmp); - const int spatialIdx = 0; using Flash = Opm::PTFlash; - Flash::solve(fluid_state, z, spatialIdx, sample, flash_tolerance, flash_verbosity); + Flash::solve(fluid_state, z, sample, flash_tolerance, flash_verbosity); ComponentVector x, y; const Evaluation L = fluid_state.L(); @@ -138,6 +137,25 @@ BOOST_DATA_TEST_CASE(PtFlash, test_methods) y[comp_idx] = fluid_state.moleFraction(FluidSystem::gasPhaseIdx, comp_idx); } + if (flash_verbosity >= 1) { + for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { + std::cout << " x for component: " << comp_idx << "is " << x[comp_idx] << std::endl; + for (int i = 0; i < 3; ++i) { + std::cout << " x deriv " << i << " is: " << x[comp_idx].derivative(i) << std::endl; + } + + std::cout << " y for component: " << comp_idx << "is " << y[comp_idx] << std::endl; + for (int i = 0; i < 3; ++i) { + std::cout << " y deriv " << i << " is: " << y[comp_idx].derivative(i) << std::endl; + } + } + std::cout << " L is " << L << std::endl; + for (int i = 0; i < L.size(); ++i) { + std::cout << " L deriv " << i << " is: " << L.derivative(i) << std::endl; + } + } + + Evaluation ref_L = 1 - 0.763309246; ref_L.setDerivative(0, 4.072857907696467e-8); ref_L.setDerivative(1, -1.1606117844565438);