diff --git a/opm/material/constraintsolvers/PTFlash.hpp b/opm/material/constraintsolvers/PTFlash.hpp index 2c47cc524..daf8765b3 100644 --- a/opm/material/constraintsolvers/PTFlash.hpp +++ b/opm/material/constraintsolvers/PTFlash.hpp @@ -65,13 +65,13 @@ class PTFlash enum { numComponents = FluidSystem::numComponents }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx}; enum { gasPhaseIdx = FluidSystem::gasPhaseIdx}; - enum { numMiscibleComponents = FluidSystem::numMiscibleComponents}; //octane, co2 // should be brine instead of brine here. + enum { numMiscibleComponents = FluidSystem::numMiscibleComponents}; enum { numMisciblePhases = FluidSystem::numMisciblePhases}; //oil, gas enum { numEq = numMisciblePhases+ numMisciblePhases*numMiscibleComponents - }; //pressure, saturation, composition + }; public: /*! @@ -82,9 +82,9 @@ public: static void solve(FluidState& fluid_state, const Dune::FieldVector& z, int spatialIdx, - int verbosity, std::string twoPhaseMethod, - Scalar tolerance = -1.) + Scalar tolerance = -1., + int verbosity = 0) { using InputEval = typename FluidState::Scalar; @@ -163,20 +163,16 @@ public: // Cell is one-phase. Use Li's phase labeling method to see if it's liquid or vapor L_scalar = li_single_phase_label_(fluid_state_scalar, z_scalar, verbosity); } + fluid_state_scalar.setLvalue(L_scalar); // Print footer if (verbosity >= 1) { std::cout << "********" << std::endl; } - // all the solution should be processed in scalar form - // now we should update the derivatives - // TODO: should be able the handle the derivatives directly, which will affect the final organization - // of the code - // TODO: Does fluid_state_scalar contain z with derivatives? - fluid_state_scalar.setLvalue(L_scalar); - fluid_state.setLvalue(L_scalar); - // ensure that things in fluid_state_scalar is transformed to fluid_state + // the flash solution process were performed in scalar form, after the flash calculation finishes, + // we transform the derivatives to the final solution + // ensure that things in fluid_state_scalar is transformed to fluid_state for (int compIdx=0; compIdx static typename Vector::field_type rachfordRice_g_(const Vector& K, typename Vector::field_type L, const Vector& z) { @@ -394,12 +390,13 @@ protected: typename Vector::field_type gLmin = rachfordRice_g_(K, Lmin, z); // Print new header - if (verbosity == 3 || verbosity == 4) { + if (verbosity >= 3) { std::cout << std::setw(10) << "Iteration" << std::setw(16) << "g(Lmid)" << std::setw(16) << "L" << std::endl; } - + + constexpr int max_it = 100; // Bisection loop - for (int iteration=1; iteration<100; ++iteration){ + for (int iteration = 0; iteration < max_it; ++iteration){ // New midpoint auto L = (Lmin + Lmax) / 2; auto gMid = rachfordRice_g_(K, L, z); @@ -423,7 +420,7 @@ protected: gLmin = gMid; } } - throw std::runtime_error(" Rachford-Rice with bisection failed!"); + throw std::runtime_error(" Rachford-Rice with bisection failed with " + std::to_string(max_it) + " iterations!"); } template @@ -599,9 +596,6 @@ protected: } } - // TODO: refactoring the template parameter later - // For the function flash_2ph, we should carry the derivatives in, and - // return with the correct and updated derivatives template static void flash_2ph(const ComponentVector& z_scalar, const std::string& flash_2p_method, @@ -653,7 +647,9 @@ protected: // Compute x and y from K, L and Z computeLiquidVapor_(fluid_state, L, K, z); - std::cout << " the current L is " << Opm::getValue(L) << std::endl; + if (verbosity >= 1) { + std::cout << " the current L is " << Opm::getValue(L) << std::endl; + } // Print initial condition if (verbosity >= 1) { @@ -732,7 +728,9 @@ protected: while (iter < max_iter) { assembleNewton_, ComponentVector, num_primary_variables, num_equations> (flash_fluid_state, z, jac, res); - std::cout << " newton residual is " << res.two_norm() << std::endl; + if (verbosity >= 1) { + std::cout << " newton residual is " << res.two_norm() << std::endl; + } converged = res.two_norm() < tolerance; if (converged) { break; @@ -763,15 +761,17 @@ protected: } } } - for (unsigned i = 0; i < num_equations; ++i) { - for (unsigned j = 0; j < num_primary_variables; ++j) { - std::cout << " " << jac[i][j] ; + if (verbosity >= 1) { + for (unsigned i = 0; i < num_equations; ++i) { + for (unsigned j = 0; j < num_primary_variables; ++j) { + std::cout << " " << jac[i][j]; + } + std::cout << std::endl; } std::cout << std::endl; } - std::cout << std::endl; if (!converged) { - throw std::runtime_error(" Newton composition update did not converge within maxIterations"); + throw std::runtime_error(" Newton composition update did not converge within maxIterations " + std::to_string(max_iter)); } // fluid_state is scalar, we need to update all the values for fluid_state here diff --git a/opm/material/eos/PengRobinsonMixture.hpp b/opm/material/eos/PengRobinsonMixture.hpp index 8a1b63574..2149e6935 100644 --- a/opm/material/eos/PengRobinsonMixture.hpp +++ b/opm/material/eos/PengRobinsonMixture.hpp @@ -86,14 +86,11 @@ public: * 4th edition, McGraw-Hill, 1987, pp. 42-44, 143-145 */ template - static LhsEval computeFugacityCoefficient(const FluidState& fs_arg, - const Params& params_arg, + static LhsEval computeFugacityCoefficient(const FluidState& fs, + const Params& params, unsigned phaseIdx, unsigned compIdx) { - auto fs = fs_arg; - auto params = params_arg; - // note that we normalize the component mole fractions, so // that their sum is 100%. This increases numerical stability // considerably if the fluid state is not physical. diff --git a/opm/material/eos/PengRobinsonParamsMixture.hpp b/opm/material/eos/PengRobinsonParamsMixture.hpp index 3243d2da8..ccac31111 100644 --- a/opm/material/eos/PengRobinsonParamsMixture.hpp +++ b/opm/material/eos/PengRobinsonParamsMixture.hpp @@ -70,7 +70,7 @@ public: /*! * \brief TODO */ - Scalar getaCache(unsigned compIIdx, unsigned compJIdx ) + Scalar getaCache(unsigned compIIdx, unsigned compJIdx ) const { return aCache_[compIIdx][compJIdx]; } diff --git a/opm/material/fluidsystems/PTFlashParameterCache.hpp b/opm/material/fluidsystems/PTFlashParameterCache.hpp index 12dcbdf6b..5056df397 100644 --- a/opm/material/fluidsystems/PTFlashParameterCache.hpp +++ b/opm/material/fluidsystems/PTFlashParameterCache.hpp @@ -2,7 +2,8 @@ // vi: set et ts=4 sw=4 sts=4: /* Copyright 2022 NORCE. - + Copyright 2022 SINTEF Digital, Mathematics and Cybernetics. + This file is part of the Open Porous Media project (OPM). OPM is free software: you can redistribute it and/or modify @@ -73,16 +74,8 @@ public: unsigned phaseIdx, int exceptQuantities = ParentType::None) { - // if (phaseIdx != oilPhaseIdx) - // return; - updateEosParams(fluidState, phaseIdx, exceptQuantities); - // if we don't need to recalculate the molar volume, we exit - // here - // if (VmUpToDate_[phaseIdx]) - // return; - // update the phase's molar volume updateMolarVolume_(fluidState, phaseIdx); } @@ -117,7 +110,7 @@ public: case gasPhaseIdx: return gasPhaseParams_.a(); default: throw std::logic_error("The a() parameter is only defined for " - "oil phases"); + "oil and gas phases"); }; } @@ -134,7 +127,7 @@ public: case gasPhaseIdx: return gasPhaseParams_.b(); default: throw std::logic_error("The b() parameter is only defined for " - "oil phase"); + "oil and gas phase"); }; } @@ -154,7 +147,7 @@ public: case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).a(); default: throw std::logic_error("The a() parameter is only defined for " - "oil phase"); + "oil and gas phase"); }; } @@ -173,7 +166,7 @@ public: case gasPhaseIdx: return gasPhaseParams_.pureParams(compIdx).b(); default: throw std::logic_error("The b() parameter is only defined for " - "oil phase"); + "oil and gas phase"); }; } @@ -186,7 +179,7 @@ public: * \param compIdx The component phase of interest * \param compJIdx Additional component index */ - Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx) + Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx) const { switch (phaseIdx) { @@ -194,7 +187,7 @@ public: case gasPhaseIdx: return gasPhaseParams_.getaCache(compIdx,compJIdx); default: throw std::logic_error("The aCache() parameter is only defined for " - "oil phase"); + "oil and gas phase"); }; } @@ -235,9 +228,6 @@ public: unsigned phaseIdx, int exceptQuantities = ParentType::None) { - // if (phaseIdx != oilPhaseIdx) - // return; - if (!(exceptQuantities & ParentType::Temperature)) { updatePure_(fluidState, phaseIdx); diff --git a/opm/material/fluidsystems/Spe5ParameterCache.hpp b/opm/material/fluidsystems/Spe5ParameterCache.hpp index 9c92fa9ce..4d27d56d7 100644 --- a/opm/material/fluidsystems/Spe5ParameterCache.hpp +++ b/opm/material/fluidsystems/Spe5ParameterCache.hpp @@ -182,7 +182,7 @@ public: * \param compIdx The component phase of interest * \param compJIdx Additional component index */ - Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx) + Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx) const { switch (phaseIdx) { diff --git a/opm/material/fluidsystems/ThreeComponentFluidSystem.hh b/opm/material/fluidsystems/ThreeComponentFluidSystem.hh index 78e800124..ba70f6949 100644 --- a/opm/material/fluidsystems/ThreeComponentFluidSystem.hh +++ b/opm/material/fluidsystems/ThreeComponentFluidSystem.hh @@ -181,7 +181,6 @@ namespace Opm { LhsEval dens; if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) { - // paramCache.updatePhase(fluidState, phaseIdx); dens = fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx); } return dens; @@ -210,11 +209,7 @@ namespace Opm { assert(phaseIdx < numPhases); assert(compIdx < numComponents); - // TODO: here the derivatives for the phi are dropped. Should we keep the derivatives against the pressure - // and temperature? - LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); - //Scalar phi = Opm::getValue( - // PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx)); + LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx); return phi; } diff --git a/tests/test_co2brine_ptflash.cpp b/tests/test_co2brine_ptflash.cpp index 03914b424..1bedfe3f2 100644 --- a/tests/test_co2brine_ptflash.cpp +++ b/tests/test_co2brine_ptflash.cpp @@ -127,7 +127,7 @@ bool testPTFlash(const std::string& flash_twophase_method) const int spatialIdx = 0; using Flash = Opm::PTFlash; - Flash::solve(fluid_state, z, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance); + Flash::solve(fluid_state, z, spatialIdx, flash_twophase_method, flash_tolerance, flash_verbosity); return result_okay(fluid_state); } @@ -229,4 +229,4 @@ int main(int argc, char **argv) } return 0; -} \ No newline at end of file +} diff --git a/tests/test_threecomponents_ptflash.cpp b/tests/test_threecomponents_ptflash.cpp index a53c7ab76..80a5cec82 100644 --- a/tests/test_threecomponents_ptflash.cpp +++ b/tests/test_threecomponents_ptflash.cpp @@ -117,7 +117,7 @@ bool testPTFlash(const std::string& flash_twophase_method) } const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional - const int flash_verbosity = 1; + const int flash_verbosity = 0; // TODO: should we set these? // Set initial K and L @@ -130,7 +130,7 @@ bool testPTFlash(const std::string& flash_twophase_method) const int spatialIdx = 0; using Flash = Opm::PTFlash; - Flash::solve(fluid_state, z, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance); + Flash::solve(fluid_state, z, spatialIdx, flash_twophase_method, flash_tolerance, flash_verbosity); return result_okay(fluid_state); } @@ -242,7 +242,9 @@ int main(int argc, char **argv) } if (!test_passed) { - throw std::runtime_error(" PTFlash tests failed"); + throw std::runtime_error(" test_threecomponents_ptflash tests FAILED"); + } else { + std::cout << "test_threecomponents_ptflash testing is SUCCESSFUL"; } return 0;