From 82c6bfe61025198458542619d4e7dac4d774d291 Mon Sep 17 00:00:00 2001 From: Trine Mykkeltvedt Date: Mon, 20 Jun 2022 09:54:15 +0200 Subject: [PATCH] added Kais startup cleaning, no functionality change --- opm/material/constraintsolvers/ChiFlash.hpp | 53 ++++++------- tests/test_chiflash.cpp | 88 ++++++++++----------- 2 files changed, 64 insertions(+), 77 deletions(-) diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp index e8ba707c5..ae12cba7b 100644 --- a/opm/material/constraintsolvers/ChiFlash.hpp +++ b/opm/material/constraintsolvers/ChiFlash.hpp @@ -58,7 +58,6 @@ namespace Opm { template class ChiFlash { - //using Problem = GetPropType; enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; enum { oilPhaseIdx = FluidSystem::oilPhaseIdx}; @@ -69,9 +68,7 @@ class ChiFlash numEq = numMisciblePhases+ numMisciblePhases*numMiscibleComponents - };//pressure, saturation, composition - - + }; //pressure, saturation, composition public: /*! @@ -84,7 +81,7 @@ public: int spatialIdx, int verbosity, std::string twoPhaseMethod, - Scalar tolerance) + Scalar tolerance = -1.) { using InputEval = typename FluidState::Scalar; @@ -93,11 +90,11 @@ public: #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7) Dune::FMatrixPrecision::set_singular_limit(1e-35); #endif + if (tolerance <= 0) { + tolerance = std::min(1e-3, 1e8 * std::numeric_limits::epsilon()); + } - if (tolerance <= 0) - tolerance = std::min(1e-3, 1e8*std::numeric_limits::epsilon()); - - //K and L from previous timestep (wilson and -1 initially) + // K and L from previous timestep (wilson and -1 initially) ComponentVector K; for(int compIdx = 0; compIdx < numComponents; ++compIdx) { K[compIdx] = fluid_state.K(compIdx); @@ -146,7 +143,6 @@ public: std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl; } phaseStabilityTest_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity); - } if (verbosity >= 1) { std::cout << "Inputs after stability test are K = [" << K_scalar << "], L = [" << L_scalar << "], z = [" << z_scalar << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl; @@ -164,11 +160,9 @@ public: } else { // 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); - single = true; + single = true; } - - // Print footer if (verbosity >= 1) { std::cout << "********" << std::endl; @@ -193,9 +187,7 @@ public: fluid_state.setKvalue(compIdx, K_scalar[compIdx]); fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]); } - updateDerivatives_(fluid_state_scalar, z, fluid_state, single); - }//end solve /*! @@ -433,7 +425,6 @@ protected: throw std::runtime_error(" Rachford-Rice with bisection failed!"); } - template static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& globalComposition, int verbosity) { @@ -475,6 +466,7 @@ protected: } } } + template static void checkStability_(const FlashFluidState& fluidState, bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc, typename FlashFluidState::Scalar& S_loc, const ComponentVector& globalComposition, bool isGas, int verbosity) @@ -537,6 +529,7 @@ protected: fluidState_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake); fluidState_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal); } + ComponentVector R; for (int compIdx=0; compIdx static void updateCurrentSol_(DefectVector& x, DefectVector& d) @@ -1021,7 +1021,6 @@ protected: for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { primary_z[comp_idx] = Opm::getValue(z[comp_idx]); } - for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx); primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii); @@ -1106,24 +1105,20 @@ protected: deri[idx] += pz * zi.derivative(idx); } } - for (unsigned idx = 0; idx < num_deri; ++idx) { x[compIdx].setDerivative(idx, deri[idx]); - } - + } // handling y for (unsigned idx = 0; idx < num_deri; ++idx) { deri[idx] = - xx[compIdx + numComponents][0]* p_v.derivative(idx); - } - + } for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { const double pz = -xx[compIdx + numComponents][cIdx + 1]; const auto& zi = z[cIdx]; for (unsigned idx = 0; idx < num_deri; ++idx) { deri[idx] += pz * zi.derivative(idx); } - } - + } for (unsigned idx = 0; idx < num_deri; ++idx) { y[compIdx].setDerivative(idx, deri[idx]); } @@ -1133,7 +1128,6 @@ protected: for (unsigned idx = 0; idx < num_deri; ++idx) { deriL[idx] = - xx[2*numComponents][0] * p_v.derivative(idx); } - for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) { const double pz = -xx[2*numComponents][cIdx + 1]; const auto& zi = z[cIdx]; @@ -1169,6 +1163,7 @@ protected: fluidState.setMoleFraction(gasPhaseIdx, compIdx, x[compIdx + numMiscibleComponents]); } + // Compute fugacities using ValueType = typename FluidState::Scalar; using ParamCache = typename FluidSystem::template ParameterCache; @@ -1198,7 +1193,7 @@ protected: // sum(x) - sum(y) = 0 b[numMiscibleComponents*numMisciblePhases] += -x[compIdx] + x[compIdx + numMiscibleComponents]; } - }//end evalDefect + }//end valDefect template static void evalJacobian_(DefectMatrix& A, @@ -1349,4 +1344,4 @@ protected: } // namespace Opm -#endif +#endif \ No newline at end of file diff --git a/tests/test_chiflash.cpp b/tests/test_chiflash.cpp index 68cce489a..cbb66fcac 100644 --- a/tests/test_chiflash.cpp +++ b/tests/test_chiflash.cpp @@ -30,21 +30,17 @@ #include #include - #include #include #include #include #include -// the following include should be removed later -// #include void testChiFlash() { - - using Scalar = double; + // TODO: the name Julia should not be there, remaining to be changed using FluidSystem = Opm::JuliaThreeComponentFluidSystem; constexpr auto numComponents = FluidSystem::numComponents; @@ -52,89 +48,85 @@ void testChiFlash() typedef Dune::FieldVector ComponentVector; typedef Opm::CompositionalFluidState FluidState; - // input +// It is a three component system + // Initial: the primary variables are, pressure, molar fractions of the first and second component Evaluation p_init = Evaluation::createVariable(10e5, 0); // 10 bar ComponentVector comp; comp[0] = Evaluation::createVariable(0.5, 1); comp[1] = Evaluation::createVariable(0.3, 2); comp[2] = 1. - comp[0] - comp[1]; - ComponentVector sat; - sat[0] = 1.0; sat[1] = 1.0-sat[0]; - // TODO: should we put the derivative against the temperature here? - Scalar temp = 300.0; - // From co2-compositional branch, it uses - // typedef typename FluidSystem::template ParameterCache ParameterCache; - - FluidState fs; - // TODO: no capillary pressure for now - fs.setPressure(FluidSystem::oilPhaseIdx, p_init); - fs.setPressure(FluidSystem::gasPhaseIdx, p_init); + // TODO: not sure whether the saturation matter here. + ComponentVector sat; + // We assume that currently everything is in the oil phase + sat[0] = 1.0; sat[1] = 1.0-sat[0]; + Scalar temp = 300.0; - fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]); - fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]); - fs.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp2Idx, comp[2]); + + // FluidState will be the input for the flash calculation + FluidState fluid_state; + fluid_state.setPressure(FluidSystem::oilPhaseIdx, p_init); + fluid_state.setPressure(FluidSystem::gasPhaseIdx, p_init); - fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]); - fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]); - fs.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]); + fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp0Idx, comp[0]); + fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp1Idx, comp[1]); + fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, FluidSystem::Comp2Idx, comp[2]); + + fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp0Idx, comp[0]); + fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp1Idx, comp[1]); + fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, FluidSystem::Comp2Idx, comp[2]); // It is used here only for calculate the z - fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]); - fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]); + fluid_state.setSaturation(FluidSystem::oilPhaseIdx, sat[0]); + fluid_state.setSaturation(FluidSystem::gasPhaseIdx, sat[1]); - fs.setTemperature(temp); + fluid_state.setTemperature(temp); // ParameterCache paramCache; { typename FluidSystem::template ParameterCache paramCache; - paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx); - paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx); - fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx)); - fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx)); + paramCache.updatePhase(fluid_state, FluidSystem::oilPhaseIdx); + paramCache.updatePhase(fluid_state, FluidSystem::gasPhaseIdx); + fluid_state.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx)); + fluid_state.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx)); } - ComponentVector zInit(0.); // TODO; zInit needs to be normalized. + ComponentVector z(0.); // TODO; z needs to be normalized. { Scalar sumMoles = 0.0; for (unsigned phaseIdx = 0; phaseIdx < FluidSystem::numPhases; ++phaseIdx) { for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - Scalar tmp = Opm::getValue(fs.molarity(phaseIdx, compIdx) * fs.saturation(phaseIdx)); - zInit[compIdx] += Opm::max(tmp, 1e-8); + Scalar tmp = Opm::getValue(fluid_state.molarity(phaseIdx, compIdx) * fluid_state.saturation(phaseIdx)); + z[compIdx] += Opm::max(tmp, 1e-8); sumMoles += tmp; } } - zInit /= sumMoles; - // initialize the derivatives - // TODO: the derivative eventually should be from the reservoir flow equations + z /= sumMoles; + // p And z is the primary variables Evaluation z_last = 1.; for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) { - zInit[compIdx] = Evaluation::createVariable(Opm::getValue(zInit[compIdx]), compIdx + 1); - z_last -= zInit[compIdx]; + z[compIdx] = Evaluation::createVariable(Opm::getValue(z[compIdx]), compIdx + 1); + z_last -= z[compIdx]; } - zInit[numComponents - 1] = z_last; + z[numComponents - 1] = z_last; } - // TODO: only, p, z need the derivatives. const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional const int flash_verbosity = 1; - //const std::string flash_twophase_method = "ssi"; // "ssi" const std::string flash_twophase_method = "newton"; - // const std::string flash_twophase_method = "ssi+newton"; // TODO: should we set these? // Set initial K and L for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { - const Evaluation Ktmp = fs.wilsonK_(compIdx); - fs.setKvalue(compIdx, Ktmp); + const Evaluation Ktmp = fluid_state.wilsonK_(compIdx); + fluid_state.setKvalue(compIdx, Ktmp); } const Evaluation Ltmp = 1.; - fs.setLvalue(Ltmp); + fluid_state.setLvalue(Ltmp); const int spatialIdx = 0; using Flash = Opm::ChiFlash; - // TODO: here the zInit does not have the proper derivatives - Flash::solve(fs, zInit, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance); + Flash::solve(fluid_state, z, spatialIdx, flash_verbosity, flash_twophase_method, flash_tolerance); } @@ -145,4 +137,4 @@ int main(int argc, char **argv) testChiFlash(); return 0; -} +} \ No newline at end of file