diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp index ffd796792..8685f0cf7 100644 --- a/opm/material/constraintsolvers/ChiFlash.hpp +++ b/opm/material/constraintsolvers/ChiFlash.hpp @@ -156,7 +156,14 @@ public: if (verbosity >= 1) { std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl; } + //phaseStabilityTestMichelsen_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity); phaseStabilityTest_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity); +/* for (int compIdx = 0; compIdx= 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; @@ -167,8 +174,10 @@ public: // Rachford Rice equation to get initial L for composition solver L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity); + Scalar Vtest = 1 - L_scalar; + - const std::string twoPhaseMethod = "newton"; // "ssi" + // const std::string twoPhaseMethod = "ssi"; // "ssi" flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity); @@ -189,9 +198,12 @@ public: // of the code // TODO: Does fluid_state_scalar contain z with derivatives? fluid_state_scalar.setLvalue(L_scalar); - //std::cout << " HEIHEIHEI L " << fluid_state_scalar.L(0) << std::endl; updateDerivatives_(fluid_state_scalar, z, fluid_state); + std::cout << " ------ SUMMARY ------ " << std::endl; + std::cout << " L " << L_scalar << std::endl; + std::cout << " K " << K_scalar << std::endl; + // Update phases /* typename FluidSystem::template ParameterCache paramCache; @@ -539,7 +551,7 @@ protected: } template - static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& z, int verbosity) + static void phaseStabilityTestMichelsen_(bool& stable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& z, int verbosity) { // Declarations bool isTrivialL, isTrivialV; @@ -552,40 +564,37 @@ protected: if (verbosity == 3 || verbosity == 4) { std::cout << "Stability test for vapor phase:" << std::endl; } - //stable_vapor, i_v = michelsen_test!(vapor, f_z, f_xy, vapor.z, z, K, eos, c, forces, Val(true); kwarg...) bool stable_vapour = false; - michelsenTest_(fluidState, z, K0,stable_vapour,/*isGas*/true, verbosity); - checkStability_(fluidState, isTrivialV, K0, y, S_v, z, /*isGas=*/true, verbosity); - bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV; + michelsenTest_(fluidState, z, y, K0,stable_vapour,/*isGas*/true, verbosity); - // Check for liquids stable phase - if (verbosity == 3 || verbosity == 4) { - std::cout << "Stability test for liquid phase:" << std::endl; + bool stable_liquid = false; + michelsenTest_(fluidState, z, x, K1,stable_liquid,/*isGas*/false, verbosity); + + //bool stable = false; + stable = stable_liquid && stable_vapour; + + if (!stable) { + for (int compIdx = 0; compIdx= 1) { + std::cout << "Stability test done for - vapour - liquid - sum:" << stable_vapour << " - " << stable_liquid << " - " << stable < - static void michelsenTest_(const FlashFluidState& fluidState, const ComponentVector z, + static void michelsenTest_(const FlashFluidState& fluidState, const ComponentVector z, ComponentVector& xy_out, ComponentVector& K, bool& stable, bool isGas, int verbosity) { using FlashEval = typename FlashFluidState::Scalar; @@ -600,14 +609,34 @@ template FlashEval S_loc = 0.0; FlashEval xy_c = 0.0; bool isTrivial; + bool isConverged; + + int phaseIdx = (isGas ? static_cast(gasPhaseIdx) : static_cast(oilPhaseIdx)); + int phaseIdx2 = (isGas ? static_cast(oilPhaseIdx) : static_cast(gasPhaseIdx)); + // Setup output if (verbosity >= 3 || verbosity >= 4) { std::cout << std::setw(10) << "Iteration" << std::setw(16) << "K-Norm" << std::setw(16) << "R-Norm" << std::endl; } + //mixture fugacity + for (int compIdx=0; compIdx paramCache_zi; + paramCache_zi.updatePhase(fluidState_zi, oilPhaseIdx); + + for (int compIdx=0; compIdx } xy_loc /= S_loc; + if (isGas) + xy_out = z; + else + xy_out = xy_loc; + //to get out fugacities each phase for (int compIdx=0; compIdx(gasPhaseIdx) : static_cast(oilPhaseIdx)); - int phaseIdx2 = (isGas ? static_cast(oilPhaseIdx) : static_cast(gasPhaseIdx)); - - for (int compIdx=0; compIdx paramCache_xy; paramCache_xy.updatePhase(fluidState_xy, phaseIdx); - typename FluidSystem::template ParameterCache paramCache_zi; - paramCache_zi.updatePhase(fluidState_zi, phaseIdx2); - for (int compIdx=0; compIdx(gasPhaseIdx) : static_cast(oilPhaseIdx)); - int phaseIdx2 = (isGas ? static_cast(oilPhaseIdx) : static_cast(gasPhaseIdx)); - - - for (int compIdx=0; compIdx paramCache_xy; - paramCache_xy.updatePhase(fluidState_xy, phaseIdx); - - typename FluidSystem::template ParameterCache paramCache_zi; - paramCache_zi.updatePhase(fluidState_zi, phaseIdx2); - - //fugacity for fake phases each component - for (int compIdx=0; compIdx // Check convergence isTrivial = (K_norm < 1e-5); - if (isTrivial || R_norm < 1e-10) + isConverged = (R_norm < 1e-10); + + bool ok = isTrivial || isConverged; + bool done = ok || i == maxIter; + + if (done && !ok) { + isTrivial = true; + throw std::runtime_error(" Stability test did not converge"); + //@warn "Stability test failed to converge in $maxiter iterations. Assuming stability." cond xy K_norm R_norm K + } + if (ok) { + stable = isTrivial || S_loc <= 1 + 1e-5; return; + } //todo: make sure that no mole fraction is smaller than 1e-8 ? - //todo: take care of water! */ + //todo: take care of water! } - stable = isTrivial || S_loc <= 1 + 1e-5; + throw std::runtime_error(" Stability test did not converge"); }//end michelsen + + template + static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& globalComposition, int verbosity) + { + // Declarations + bool isTrivialL, isTrivialV; + ComponentVector x, y; + typename FlashFluidState::Scalar S_l, S_v; + ComponentVector K0 = K; + ComponentVector K1 = K; + + // Check for vapour instable phase + if (verbosity == 3 || verbosity == 4) { + std::cout << "Stability test for vapor phase:" << std::endl; + } + checkStability_(fluidState, isTrivialV, K0, y, S_v, globalComposition, /*isGas=*/true, verbosity); + bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV; + + // Check for liquids stable phase + if (verbosity == 3 || verbosity == 4) { + std::cout << "Stability test for liquid phase:" << std::endl; + } + checkStability_(fluidState, isTrivialL, K1, x, S_l, globalComposition, /*isGas=*/false, verbosity); + bool L_stable = (S_l < (1.0 + 1e-5)) || isTrivialL; + + // L-stable means success in making liquid, V-unstable means no success in making vapour + isStable = L_stable && V_unstable; + if (isStable) { + // Single phase, i.e. phase composition is equivalent to the global composition + // Update fluidstate with mole fraction + for (int compIdx=0; compIdx 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) @@ -1266,22 +1294,22 @@ template (primary_fluid_state, primary_z, pri_jac, pri_res); //corresponds to julias J_p (we miss d/dt, and have d/dL instead of d/dV) - // for (unsigned i =0; i < num_equations; ++i) { - // for (unsigned j = 0; j < primary_num_pv; ++j) { - // std::cout << " " << pri_jac[i][j]; - // } - // std::cout << std::endl; - // } - // std::cout << std::endl; + for (unsigned i =0; i < num_equations; ++i) { + for (unsigned j = 0; j < primary_num_pv; ++j) { + std::cout << " " << pri_jac[i][j]; + } + std::cout << std::endl; + } + std::cout << std::endl; //corresponds to julias J_s - // for (unsigned i = 0; i < num_equations; ++i) { - // for (unsigned j = 0; j < secondary_num_pv; ++j) { - // std::cout << " " << sec_jac[i][j] ; - // } - // std::cout << std::endl; - // } - // std::cout << std::endl; + for (unsigned i = 0; i < num_equations; ++i) { + for (unsigned j = 0; j < secondary_num_pv; ++j) { + std::cout << " " << sec_jac[i][j] ; + } + std::cout << std::endl; + } + std::cout << std::endl; SecondaryNewtonMatrix xx; pri_jac.solve(xx,sec_jac); @@ -1616,6 +1644,114 @@ template } // throw std::runtime_error("Successive substitution composition update did not converge within maxIterations"); } + + template + static void successiveSubstitutionCompositionNew_(ComponentVector& K, typename FlashFluidState::Scalar& L, FlashFluidState& fluidState, const ComponentVector& z, + const bool newton_afterwards, const int verbosity) + { + // std::cout << " Evaluation in successiveSubstitutionComposition_ is " << Dune::className(L) << std::endl; + // 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; + + // Store cout format before manipulation + std::ios_base::fmtflags f(std::cout.flags()); + + // Print initial guess + if (verbosity >= 1) + std::cout << "Initial guess: K = [" << K << "] and L = " << L << std::endl; + + if (verbosity == 2 || verbosity == 4) { + // Print header + int fugWidth = (numComponents * 12)/2; + int convWidth = fugWidth + 7; + std::cout << std::setw(10) << "Iteration" << std::setw(fugWidth) << "fL/fV" << std::setw(convWidth) << "norm2(fL/fv-1)" << std::endl; + } + // + // Successive substitution loop + // + for (int i=0; i < maxIterations; ++i){ + // Compute (normalized) liquid and vapor mole fractions + computeLiquidVapor_(fluidState, L, K, z); + + // Calculate fugacity coefficient + using ParamCache = typename FluidSystem::template ParameterCache; + ParamCache paramCache; + for (int phaseIdx=0; phaseIdx= 1) { + std::cout << "Solution converged to the following result :" << std::endl; + std::cout << "x = ["; + for (int compIdx=0; compIdx