From 7fa4f4796fb0a56544e1679af17393bce65bcb14 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 5 Nov 2021 15:21:35 +0100 Subject: [PATCH] [WIP] Add ChiFlash and test. Not compiling. --- opm/material/constraintsolvers/ChiFlash.hpp | 966 ++++++++++++++++++++ tests/test_chiflash.cpp | 621 +++++++++++++ 2 files changed, 1587 insertions(+) create mode 100644 opm/material/constraintsolvers/ChiFlash.hpp create mode 100644 tests/test_chiflash.cpp diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp new file mode 100644 index 000000000..5f51aea06 --- /dev/null +++ b/opm/material/constraintsolvers/ChiFlash.hpp @@ -0,0 +1,966 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/* + 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 2 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 . + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * \copydoc Opm::ChiFlash + */ +#ifndef OPM_CHI_FLASH_HPP +#define OPM_CHI_FLASH_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include +#include +#include + +namespace Opm { + +/*! + * \brief Determines the phase compositions, pressures and saturations + * given the total mass of all components for the chiwoms problem. + * + */ +template +class ChiFlash +{ + //using Problem = GetPropType; + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + // enum { Comp0Idx = FluidSystem::Comp0Idx }; //rename for generic ? + // enum { Comp1Idx = FluidSystem::Comp1Idx }; //rename for generic ? + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx}; + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx}; + enum { numMiscibleComponents = 2}; //octane, co2 + enum { numMisciblePhases = 2}; //oil, gas + enum { + numEq = + numMisciblePhases+ + numMisciblePhases*numMiscibleComponents + };//pressure, saturation, composition + + enum { + p0PvIdx = 0, // pressure first phase primary variable index + S0PvIdx = 1, // saturation first phase primary variable index + x00PvIdx = S0PvIdx + 1, // molefraction first phase first component primary variable index + //numMiscibleComponennets*numMisciblePhases-1 molefractions/primvar follow + }; + +public: + /*! + * \brief Calculates the fluid state from the global mole fractions of the components and the phase pressures + * + */ + template + static void solve(FluidState& fluidState, + const Dune::FieldVector& globalComposition, + int spatialIdx, + int verbosity, + std::string twoPhaseMethod, + Scalar tolerance) + { + + using InputEval = typename FluidState::Scalar; + using Matrix = Dune::FieldMatrix; + using Vetor = Dune::FieldVector; + using FlashEval = Opm::DenseAd::Evaluation; + using FlashDefectVector = Dune::FieldVector; + using FlashFluidState = Opm::CompositionalFluidState; + + using ComponentVector = Dune::FieldVector; + +#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()); + + //K and L from previous timestep (wilson and -1 initially) + ComponentVector K; + for(int compIdx = 0; compIdx < numComponents; ++compIdx) { + K[compIdx] = fluidState.K(compIdx); + } + InputEval L; + L = fluidState.L(0); + + // 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 = [" << globalComposition << "], P = " << fluidState.pressure(0) << ", and T = " << fluidState.temperature(0) << std::endl; + } + + // Do a stability test to check if cell is single-phase (do for all cells the first time). + bool isStable = false; + if ( L <= 0 || L == 1 ) { + if (verbosity >= 1) { + std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl; + } + phaseStabilityTest_(isStable, K, fluidState, globalComposition, verbosity); + } + + // Update the composition if cell is two-phase + if (isStable == false) { + + // Print info + if (verbosity >= 1) { + std::cout << "Cell is two-phase! Solve Rachford-Rice with initial K = [" << K << "]" << std::endl; + } + + // Rachford Rice equation to get initial L for composition solver + L = solveRachfordRice_g_(K, globalComposition, verbosity); + + // Calculate composition using nonlinear solver + // Newton + if (twoPhaseMethod == "newton"){ + if (verbosity >= 1) { + std::cout << "Calculate composition using Newton." << std::endl; + } + newtonCompositionUpdate_(K, L, fluidState, globalComposition, verbosity); + + } + + // Successive substitution + else if (twoPhaseMethod == "ssi"){ + if (verbosity >= 1) { + std::cout << "Calculate composition using Succcessive Substitution." << std::endl; + } + successiveSubstitutionComposition_(K, L, fluidState, globalComposition, /*standAlone=*/true, verbosity); + } + } + + // Cell is one-phase. Use Li's phase labeling method to see if it's liquid or vapor + else{ + L = li_single_phase_label_(fluidState, globalComposition, verbosity); + } + + // Print footer + if (verbosity >= 1) { + std::cout << "********" << std::endl; + } + + // Update phases + typename FluidSystem::template ParameterCache paramCache; + paramCache.updatePhase(fluidState, oilPhaseIdx); + paramCache.updatePhase(fluidState, gasPhaseIdx); + + // Calculate compressibility factor + const Scalar R = Opm::Constants::R; + Evaluation Z_L = (paramCache.molarVolume(oilPhaseIdx) * fluidState.pressure(oilPhaseIdx) )/ + (R * fluidState.temperature(oilPhaseIdx)); + Evaluation Z_V = (paramCache.molarVolume(gasPhaseIdx) * fluidState.pressure(gasPhaseIdx) )/ + (R * fluidState.temperature(gasPhaseIdx)); + + // Update saturation + Evaluation So = L*Z_L/(L*Z_L+(1-L)*Z_V); + Evaluation Sg = 1-So; + + fluidState.setSaturation(oilPhaseIdx, So); + fluidState.setSaturation(gasPhaseIdx, Sg); + + //Update L and K to the problem for the next flash + for (int compIdx = 0; compIdx < numComponents; ++compIdx){ + fluidState.setKvalue(compIdx, K[compIdx]); + } + fluidState.setLvalue(L); + + + // Print saturation + if (verbosity == 5) { + std::cout << "So = " << So <= 1) { + std::cout << "Cell is single-phase, vapor (L = 0.0) due to Li's phase labeling method giving T >= Tc_est (" << T << " >= " << Tc_est << ")!" << std::endl; + } + } + + + return L; + } + + template + static typename Vector::field_type rachfordRice_g_(const Vector& K, const Evaluation L, const Vector& globalComposition) + { + typename Vector::field_type g=0; + for (int compIdx=0; compIdx + static typename Vector::field_type rachfordRice_dg_dL_(const Vector& K, const Evaluation L, const Vector& globalComposition) + { + typename Vector::field_type dg=0; + for (int compIdx=0; compIdx + static typename Vector::field_type solveRachfordRice_g_(const Vector& K, const Vector& globalComposition, int verbosity) + { + // Find min and max K. Have to do a laborious for loop to avoid water component (where K=0) + // TODO: Replace loop with Dune::min_value() and Dune::max_value() when water component is properly handled + Evaluation Kmin = K[0]; + Evaluation Kmax = K[0]; + for (int compIdx=1; compIdx= Kmax) + Kmax = K[compIdx]; + } + + // Lower and upper bound for solution + Evaluation Lmin = (Kmin / (Kmin - 1)); + Evaluation Lmax = Kmax / (Kmax - 1); + + // Check if Lmin < Lmax, and switch if not + if (Lmin > Lmax) + { + Evaluation Ltmp = Lmin; + Lmin = Lmax; + Lmax = Ltmp; + } + + // Initial guess + Evaluation L = (Lmin + Lmax)/2; + + // Print initial guess and header + if (verbosity == 3 || verbosity == 4) { + std::cout << "Initial guess: L = " << L << " and [Lmin, Lmax] = [" << Lmin << ", " << Lmax << "]" << std::endl; + std::cout << std::setw(10) << "Iteration" << std::setw(16) << "abs(step)" << std::setw(16) << "L" << std::endl; + } + + // Newton-Raphson loop + for (int iteration=1; iteration<100; ++iteration){ + // Calculate function and derivative values + Evaluation g = rachfordRice_g_(K, L, globalComposition); + Evaluation dg_dL = rachfordRice_dg_dL_(K, L, globalComposition); + + // Lnew = Lold - g/dg; + Evaluation delta = g/dg_dL; + L -= delta; + + // Check if L is within the bounds, and if not, we apply bisection method + if (L < Lmin || L > Lmax) + { + // Print info + if (verbosity == 3 || verbosity == 4) { + std::cout << "L is not within the the range [Lmin, Lmax], solve using Bisection method!" << std::endl; + } + + // Run bisection + L = bisection_g_(K, Lmin, Lmax, globalComposition, verbosity); + + // Ensure that L is in the range (0, 1) + L = Opm::min(Opm::max(L, 0.0), 1.0); + + // Print final result + if (verbosity >= 1) { + std::cout << "Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl; + } + return L; + } + + // Print iteration info + if (verbosity == 3 || verbosity == 4) { + std::cout << std::setw(10) << iteration << std::setw(16) << Opm::abs(delta) << std::setw(16) << L << std::endl; + } + // Check for convergence + if ( Opm::abs(delta) < 1e-10 ) { + // Ensure that L is in the range (0, 1) + L = Opm::min(Opm::max(L, 0.0), 1.0); + + // Print final result + if (verbosity >= 1) { + std::cout << "Rachford-Rice converged to final solution L = " << L << std::endl; + } + return L; + } + } + // Throw error if Rachford-Rice fails + throw std::runtime_error(" Rachford-Rice did not converge within maximum number of iterations" ); + } + + template + static typename Vector::field_type bisection_g_(const Vector& K, Evaluation Lmin, Evaluation Lmax, const Vector& globalComposition, int verbosity) + { + // Calculate for g(Lmin) for first comparison with gMid = g(L) + Evaluation gLmin = rachfordRice_g_(K, Lmin, globalComposition); + + // Print new header + if (verbosity == 3 || verbosity == 4) { + std::cout << std::setw(10) << "Iteration" << std::setw(16) << "g(Lmid)" << std::setw(16) << "L" << std::endl; + } + + // Bisection loop + for (int iteration=1; iteration<100; ++iteration){ + // New midpoint + Evaluation L = (Lmin + Lmax) / 2; + Evaluation gMid = rachfordRice_g_(K, L, globalComposition); + if (verbosity == 3 || verbosity == 4) { + std::cout << std::setw(10) << iteration << std::setw(16) << gMid << std::setw(16) << L << std::endl; + } + + // Check if midpoint fulfills g=0 or L - Lmin is sufficiently small + if (Opm::abs(gMid) < 1e-10 || Opm::abs((Lmax - Lmin) / 2) < 1e-10){ + return L; + } + + // Else we repeat with midpoint being either Lmin og Lmax (depending on the signs). + else if (Dune::sign(gMid) != Dune::sign(gLmin)) { + // gMid has same sign as gLmax, so we set L as the new Lmax + Lmax = L; + } + else { + // gMid and gLmin have same sign so we set L as the new Lmin + Lmin = L; + gLmin = gMid; + } + } + throw std::runtime_error(" Rachford-Rice with bisection failed!"); + } + + template + static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& globalComposition, int verbosity) + { + // Declarations + bool isTrivialL, isTrivialV; + ComponentVector x, y; + Evaluation 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 fration + for (int compIdx=0; compIdx + static void checkStability_(const FlashFluidState& fluidState, bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc, Evaluation& S_loc, const ComponentVector& globalComposition, bool isGas, int verbosity) + { + using FlashEval = typename FlashFluidState::Scalar; + using PengRobinsonMixture = typename Opm::PengRobinsonMixture; + + // Declarations + FlashFluidState fluidState_fake = fluidState; + FlashFluidState fluidState_global = fluidState; + + // 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; + } + + // Michelsens stability test. + // Make two fake phases "inside" one phase and check for positive volume + for (int i = 0; i < 20000; ++i) { + S_loc = 0.0; + if (isGas) { + for (int compIdx=0; compIdx paramCache_fake; + paramCache_fake.updatePhase(fluidState_fake, phaseIdx); + + typename FluidSystem::template ParameterCache paramCache_global; + paramCache_global.updatePhase(fluidState_global, phaseIdx2); + + //fugacity for fake phases each component + for (int compIdx=0; compIdx + static void computeLiquidVapor_(FlashFluidState& fluidState, Evaluation& L, ComponentVector& K, const ComponentVector& globalComposition) + { + // Calculate x and y, and normalize + ComponentVector x; + ComponentVector y; + Evaluation sumx=0; + Evaluation sumy=0; + for (int compIdx=0; compIdx + static void newtonCompositionUpdate_(ComponentVector& K, Evaluation& L, FlashFluidState& fluidState, const ComponentVector& globalComposition, + int verbosity) + { + // Newton declarations + using NewtonVector = Dune::FieldVector; + using NewtonMatrix = Dune::FieldMatrix; + NewtonVector newtonX; + NewtonVector newtonB; + NewtonMatrix newtonA; + NewtonVector newtonDelta; + + // Compute x and y from K, L and Z + computeLiquidVapor_(fluidState, L, K, globalComposition); + + // Print initial condition + if (verbosity >= 1) { + std::cout << "Initial guess: x = ["; + for (int compIdx=0; compIdx= 1) { + std::cout << "Solution converged to the following result :" << std::endl; + std::cout << "x = ["; + for (int compIdx=0; compIdx + static void updateCurrentSol_(DefectVector& x, DefectVector& d) + { + // Find smallest percentage update + Scalar w = 1.0; + for (int i=0; i + static bool checkFugacityEquil_(DefectVector& b) + { + // Init. fugacity vector + DefectVector fugVec; + + // Loop over b and find the fugacity equilibrium + // OBS: If the equations in b changes in evalDefect_ then you have to change here as well! + for (int compIdx=0; compIdx + static void evalDefect_(DefectVector& b, + DefectVector& x, + const FluidState& fluidStateIn, + const ComponentVector& globalComposition) + { + // Put x and y in a FluidState instance for fugacity calculations + FluidState fluidState(fluidStateIn); + ComponentVector K; + for (int compIdx=0; compIdx; + ParamCache paramCache; + for (int phaseIdx=0; phaseIdx + static void evalJacobian_(DefectMatrix& A, + const DefectVector& xIn, + const FluidState& fluidStateIn, + const ComponentVector& globalComposition) + { + // TODO: Use AD instead + // Calculate response of current state x + DefectVector x; + DefectVector b0; + for(int j=0; j + static void successiveSubstitutionComposition_(ComponentVector& K, Evaluation& L, FlashFluidState& fluidState, const ComponentVector& globalComposition, + bool standAlone, 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. + int maxIterations; + if (standAlone == true) + maxIterations = 100; + else + maxIterations = 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, globalComposition); + + // 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. + + Consult the COPYING file in the top-level source directory of this + module for the precise wording of the license and the list of + copyright holders. +*/ +/*! + * \file + * + * \brief This is test for the SPE5 fluid system (which uses the + * Peng-Robinson EOS) and the NCP flash solver. + */ +#include "config.h" + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +template +void createSurfaceGasFluidSystem(FluidState& gasFluidState) +{ + static const int gasPhaseIdx = FluidSystem::gasPhaseIdx; + + // temperature + gasFluidState.setTemperature(273.15 + 20); + + // gas pressure + gasFluidState.setPressure(gasPhaseIdx, 1e5); + + // gas saturation + gasFluidState.setSaturation(gasPhaseIdx, 1.0); + + // gas composition: mostly methane, a bit of propane + gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::H2OIdx, 0.0); + gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C1Idx, 0.94); + gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C3Idx, 0.06); + gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C6Idx, 0.00); + gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C10Idx, 0.00); + gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C15Idx, 0.00); + gasFluidState.setMoleFraction(gasPhaseIdx, FluidSystem::C20Idx, 0.00); + + // gas density + typename FluidSystem::template ParameterCache paramCache; + paramCache.updatePhase(gasFluidState, gasPhaseIdx); + gasFluidState.setDensity(gasPhaseIdx, + FluidSystem::density(gasFluidState, paramCache, gasPhaseIdx)); +} + +template +Scalar computeSumxg(FluidState& resultFluidState, + const FluidState& prestineFluidState, + const FluidState& gasFluidState, + Scalar additionalGas) +{ + static const int oilPhaseIdx = FluidSystem::oilPhaseIdx; + static const int gasPhaseIdx = FluidSystem::gasPhaseIdx; + static const int numComponents = FluidSystem::numComponents; + + typedef Dune::FieldVector ComponentVector; + typedef Opm::NcpFlash Flash; + + resultFluidState.assign(prestineFluidState); + + // add a bit of additional gas components + ComponentVector totalMolarities; + for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++ compIdx) + totalMolarities = + prestineFluidState.molarity(oilPhaseIdx, compIdx) + + additionalGas*gasFluidState.moleFraction(gasPhaseIdx, compIdx); + + // "flash" the modified fluid state + typename FluidSystem::ParameterCache paramCache; + Flash::solve(resultFluidState, totalMolarities); + + Scalar sumxg = 0; + for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) + sumxg += resultFluidState.moleFraction(gasPhaseIdx, compIdx); + + return sumxg; +} + +template +void makeOilSaturated(FluidState& fluidState, const FluidState& gasFluidState) +{ + static const int gasPhaseIdx = FluidSystem::gasPhaseIdx; + + FluidState prestineFluidState; + prestineFluidState.assign(fluidState); + + Scalar sumxg = 0; + for (unsigned compIdx = 0; compIdx < FluidSystem::numComponents; ++compIdx) + sumxg += fluidState.moleFraction(gasPhaseIdx, compIdx); + + // Newton method + Scalar tol = 1e-8; + Scalar additionalGas = 0; // [mol] + for (int i = 0; std::abs(sumxg - 1) > tol; ++i) { + if (i > 50) + throw std::runtime_error("Newton method did not converge after 50 iterations"); + + Scalar eps = std::max(1e-8, additionalGas*1e-8); + + Scalar f = 1 - computeSumxg(prestineFluidState, + fluidState, + gasFluidState, + additionalGas); + Scalar fStar = 1 - computeSumxg(prestineFluidState, + fluidState, + gasFluidState, + additionalGas + eps); + Scalar fPrime = (fStar - f)/eps; + + additionalGas -= f/fPrime; + }; +} + +template +void guessInitial(FluidState& fluidState, unsigned phaseIdx) +{ + if (phaseIdx == FluidSystem::gasPhaseIdx) { + fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.74785); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.0121364); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.00606028); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.00268136); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.000204256); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 8.78291e-06); + } + else if (phaseIdx == FluidSystem::oilPhaseIdx) { + fluidState.setMoleFraction(phaseIdx, FluidSystem::H2OIdx, 0.0); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C1Idx, 0.50); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C3Idx, 0.03); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C6Idx, 0.07); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C10Idx, 0.20); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C15Idx, 0.15); + fluidState.setMoleFraction(phaseIdx, FluidSystem::C20Idx, 0.05); + } + else { + assert(phaseIdx == FluidSystem::waterPhaseIdx); + } +} + +template +Scalar bringOilToSurface(FluidState& surfaceFluidState, Scalar alpha, const FluidState& reservoirFluidState, bool guessInitial) +{ + enum { + numPhases = FluidSystem::numPhases, + waterPhaseIdx = FluidSystem::waterPhaseIdx, + gasPhaseIdx = FluidSystem::gasPhaseIdx, + oilPhaseIdx = FluidSystem::oilPhaseIdx, + + numComponents = FluidSystem::numComponents + }; + + typedef Opm::NcpFlash Flash; + typedef Opm::ThreePhaseMaterialTraits MaterialTraits; + typedef Opm::LinearMaterial MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + typedef Dune::FieldVector ComponentVector; + + const Scalar refPressure = 1.0135e5; // [Pa] + + // set the parameters for the capillary pressure law + MaterialLawParams matParams; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + matParams.setPcMinSat(phaseIdx, 0.0); + matParams.setPcMaxSat(phaseIdx, 0.0); + } + matParams.finalize(); + + // retieve the global volumetric component molarities + surfaceFluidState.setTemperature(273.15 + 20); + + ComponentVector molarities; + for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) + molarities[compIdx] = reservoirFluidState.molarity(oilPhaseIdx, compIdx); + + if (guessInitial) { + // we start at a fluid state with reservoir oil. + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { + for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) { + surfaceFluidState.setMoleFraction(phaseIdx, + compIdx, + reservoirFluidState.moleFraction(phaseIdx, compIdx)); + } + surfaceFluidState.setDensity(phaseIdx, reservoirFluidState.density(phaseIdx)); + surfaceFluidState.setPressure(phaseIdx, reservoirFluidState.pressure(phaseIdx)); + surfaceFluidState.setSaturation(phaseIdx, 0.0); + } + surfaceFluidState.setSaturation(oilPhaseIdx, 1.0); + surfaceFluidState.setSaturation(gasPhaseIdx, 1.0 - surfaceFluidState.saturation(oilPhaseIdx)); + } + + typename FluidSystem::template ParameterCache paramCache; + paramCache.updateAll(surfaceFluidState); + + // increase volume until we are at surface pressure. use the + // newton method for this + ComponentVector tmpMolarities; + for (int i = 0;; ++i) { + if (i >= 20) + throw Opm::NumericalIssue("Newton method did not converge after 20 iterations"); + + // calculate the deviation from the standard pressure + tmpMolarities = molarities; + tmpMolarities /= alpha; + Flash::template solve(surfaceFluidState, matParams, paramCache, tmpMolarities); + Scalar f = surfaceFluidState.pressure(gasPhaseIdx) - refPressure; + + // calculate the derivative of the deviation from the standard + // pressure + Scalar eps = alpha*1e-10; + tmpMolarities = molarities; + tmpMolarities /= alpha + eps; + Flash::template solve(surfaceFluidState, matParams, paramCache, tmpMolarities); + Scalar fStar = surfaceFluidState.pressure(gasPhaseIdx) - refPressure; + Scalar fPrime = (fStar - f)/eps; + + // newton update + Scalar delta = f/fPrime; + alpha -= delta; + if (std::abs(delta) < std::abs(alpha)*1e-9) { + break; + } + } + + // calculate the final result + tmpMolarities = molarities; + tmpMolarities /= alpha; + Flash::template solve(surfaceFluidState, matParams, paramCache, tmpMolarities); + return alpha; +} + +template +void printResult(const RawTable& rawTable, + const std::string& fieldName, + size_t firstIdx, + size_t secondIdx, + double hiresThres) +{ + std::cout << "std::vector > "< hiresThres; ++numRawHires) + {} + + for (; sampleIdx < numSamples; ++sampleIdx) { + size_t rawIdx = sampleIdx*numRawHires/numSamples; + std::cout << "{ " << rawTable[rawIdx][firstIdx] << ", " + << rawTable[rawIdx][secondIdx] << " }" + << ",\n"; + } + + numSamples = 15; + for (sampleIdx = 0; sampleIdx < numSamples; ++sampleIdx) { + size_t rawIdx = sampleIdx*(rawTable.size() - numRawHires)/numSamples + numRawHires; + std::cout << "{ " << rawTable[rawIdx][firstIdx] << ", " + << rawTable[rawIdx][secondIdx] << " }"; + if (sampleIdx < numSamples - 1) + std::cout << ",\n"; + else + std::cout << "\n"; + } + + std::cout << "};\n"; +} + +template +inline void testAll() +{ + typedef Opm::Spe5FluidSystem FluidSystem; + + enum { + numPhases = FluidSystem::numPhases, + waterPhaseIdx = FluidSystem::waterPhaseIdx, + gasPhaseIdx = FluidSystem::gasPhaseIdx, + oilPhaseIdx = FluidSystem::oilPhaseIdx, + + numComponents = FluidSystem::numComponents, + H2OIdx = FluidSystem::H2OIdx, + C1Idx = FluidSystem::C1Idx, + C3Idx = FluidSystem::C3Idx, + C6Idx = FluidSystem::C6Idx, + C10Idx = FluidSystem::C10Idx, + C15Idx = FluidSystem::C15Idx, + C20Idx = FluidSystem::C20Idx + }; + + typedef Opm::NcpFlash Flash; + typedef Dune::FieldVector ComponentVector; + typedef Opm::CompositionalFluidState FluidState; + + typedef Opm::ThreePhaseMaterialTraits MaterialTraits; + typedef Opm::LinearMaterial MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + + typedef typename FluidSystem::template ParameterCache ParameterCache; + + //////////// + // Initialize the fluid system and create the capillary pressure + // parameters + //////////// + Scalar T = 273.15 + 20; // 20 deg Celsius + FluidSystem::init(/*minTemperature=*/T - 1, + /*maxTemperature=*/T + 1, + /*minPressure=*/1.0e4, + /*maxTemperature=*/40.0e6); + + // set the parameters for the capillary pressure law + MaterialLawParams matParams; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + matParams.setPcMinSat(phaseIdx, 0.0); + matParams.setPcMaxSat(phaseIdx, 0.0); + } + matParams.finalize(); + + //////////// + // Create a fluid state + //////////// + FluidState gasFluidState; + createSurfaceGasFluidSystem(gasFluidState); + + FluidState fluidState; + ParameterCache paramCache; + + // temperature + fluidState.setTemperature(T); + + // oil pressure + fluidState.setPressure(oilPhaseIdx, 4000 * 6894.7573); // 4000 PSI + + // oil saturation + fluidState.setSaturation(oilPhaseIdx, 1.0); + fluidState.setSaturation(gasPhaseIdx, 1.0 - fluidState.saturation(oilPhaseIdx)); + + // oil composition: SPE-5 reservoir oil + fluidState.setMoleFraction(oilPhaseIdx, H2OIdx, 0.0); + fluidState.setMoleFraction(oilPhaseIdx, C1Idx, 0.50); + fluidState.setMoleFraction(oilPhaseIdx, C3Idx, 0.03); + fluidState.setMoleFraction(oilPhaseIdx, C6Idx, 0.07); + fluidState.setMoleFraction(oilPhaseIdx, C10Idx, 0.20); + fluidState.setMoleFraction(oilPhaseIdx, C15Idx, 0.15); + fluidState.setMoleFraction(oilPhaseIdx, C20Idx, 0.05); + + //makeOilSaturated(fluidState, gasFluidState); + + // set the saturations and pressures of the other phases + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (phaseIdx != oilPhaseIdx) { + fluidState.setSaturation(phaseIdx, 0.0); + fluidState.setPressure(phaseIdx, fluidState.pressure(oilPhaseIdx)); + } + + // initial guess for the composition (needed by the ComputeFromReferencePhase + // constraint solver. TODO: bug in ComputeFromReferencePhase?) + guessInitial(fluidState, phaseIdx); + } + + typedef Opm::ComputeFromReferencePhase CFRP; + CFRP::solve(fluidState, + paramCache, + /*refPhaseIdx=*/oilPhaseIdx, + /*setViscosity=*/false, + /*setEnthalpy=*/false); + + //////////// + // Calculate the total molarities of the components + //////////// + ComponentVector totalMolarities; + for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) + totalMolarities[compIdx] = fluidState.saturation(oilPhaseIdx)*fluidState.molarity(oilPhaseIdx, compIdx); + + //////////// + // Gradually increase the volume for and calculate the gas + // formation factor, oil formation volume factor and gas formation + // volume factor. + //////////// + + FluidState flashFluidState, surfaceFluidState; + flashFluidState.assign(fluidState); + //Flash::guessInitial(flashFluidState, totalMolarities); + Flash::template solve(flashFluidState, matParams, paramCache, totalMolarities); + + Scalar surfaceAlpha = 1; + surfaceAlpha = bringOilToSurface(surfaceFluidState, surfaceAlpha, flashFluidState, /*guessInitial=*/true); + Scalar rho_gRef = surfaceFluidState.density(gasPhaseIdx); + Scalar rho_oRef = surfaceFluidState.density(oilPhaseIdx); + + std::vector > resultTable; + + Scalar minAlpha = 0.98; + Scalar maxAlpha = surfaceAlpha; + + std::cout << "alpha[-] p[Pa] S_g[-] rho_o[kg/m^3] rho_g[kg/m^3] [kg/mol] [kg/mol] R_s[m^3/m^3] B_g[-] B_o[-]\n"; + int n = 300; + for (int i = 0; i < n; ++i) { + // ratio between the original and the current volume + Scalar alpha = minAlpha + (maxAlpha - minAlpha)*i/(n - 1); + + // increasing the volume means decreasing the molartity + ComponentVector curTotalMolarities = totalMolarities; + curTotalMolarities /= alpha; + + // "flash" the modified reservoir oil + Flash::template solve(flashFluidState, matParams, paramCache, curTotalMolarities); + + surfaceAlpha = bringOilToSurface(surfaceFluidState, + surfaceAlpha, + flashFluidState, + /*guessInitial=*/false); + Scalar Rs = + surfaceFluidState.saturation(gasPhaseIdx) + / surfaceFluidState.saturation(oilPhaseIdx); + std::cout << alpha << " " + << flashFluidState.pressure(oilPhaseIdx) << " " + << flashFluidState.saturation(gasPhaseIdx) << " " + << flashFluidState.density(oilPhaseIdx) << " " + << flashFluidState.density(gasPhaseIdx) << " " + << flashFluidState.averageMolarMass(oilPhaseIdx) << " " + << flashFluidState.averageMolarMass(gasPhaseIdx) << " " + << Rs << " " + << rho_gRef/flashFluidState.density(gasPhaseIdx) << " " + << rho_oRef/flashFluidState.density(oilPhaseIdx) << " " + << "\n"; + + std::array tmp; + tmp[0] = alpha; + tmp[1] = flashFluidState.pressure(oilPhaseIdx); + tmp[2] = flashFluidState.saturation(gasPhaseIdx); + tmp[3] = flashFluidState.density(oilPhaseIdx); + tmp[4] = flashFluidState.density(gasPhaseIdx); + tmp[5] = flashFluidState.averageMolarMass(oilPhaseIdx); + tmp[6] = flashFluidState.averageMolarMass(gasPhaseIdx); + tmp[7] = Rs; + tmp[8] = rho_gRef/flashFluidState.density(gasPhaseIdx); + tmp[9] = rho_oRef/flashFluidState.density(oilPhaseIdx); + + resultTable.push_back(tmp); + } + + std::cout << "reference density oil [kg/m^3]: " << rho_oRef << "\n"; + std::cout << "reference density gas [kg/m^3]: " << rho_gRef << "\n"; + + Scalar hiresThresholdPressure = resultTable[20][1]; + printResult(resultTable, + "Bg", /*firstIdx=*/1, /*secondIdx=*/8, + /*hiresThreshold=*/hiresThresholdPressure); + printResult(resultTable, + "Bo", /*firstIdx=*/1, /*secondIdx=*/9, + /*hiresThreshold=*/hiresThresholdPressure); + printResult(resultTable, + "Rs", /*firstIdx=*/1, /*secondIdx=*/7, + /*hiresThreshold=*/hiresThresholdPressure); +} + +void testChiFlash() +{ + using Scalar = double; + typedef Opm::Spe5FluidSystem FluidSystem; + + enum { + numPhases = FluidSystem::numPhases, + waterPhaseIdx = FluidSystem::waterPhaseIdx, + gasPhaseIdx = FluidSystem::gasPhaseIdx, + oilPhaseIdx = FluidSystem::oilPhaseIdx, + + numComponents = FluidSystem::numComponents, + H2OIdx = FluidSystem::H2OIdx, + C1Idx = FluidSystem::C1Idx, + C3Idx = FluidSystem::C3Idx, + C6Idx = FluidSystem::C6Idx, + C10Idx = FluidSystem::C10Idx, + C15Idx = FluidSystem::C15Idx, + C20Idx = FluidSystem::C20Idx + }; + + //typedef Opm::NcpFlash Flash; + typedef Dune::FieldVector ComponentVector; + typedef Opm::CompositionalFluidState FluidState; + + typedef Opm::ThreePhaseMaterialTraits MaterialTraits; + typedef Opm::LinearMaterial MaterialLaw; + typedef typename MaterialLaw::Params MaterialLawParams; + + typedef typename FluidSystem::template ParameterCache ParameterCache; + + //////////// + // Initialize the fluid system and create the capillary pressure + // parameters + //////////// + Scalar T = 273.15 + 20; // 20 deg Celsius + FluidSystem::init(/*minTemperature=*/T - 1, + /*maxTemperature=*/T + 1, + /*minPressure=*/1.0e4, + /*maxTemperature=*/40.0e6); + + // set the parameters for the capillary pressure law + MaterialLawParams matParams; + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + matParams.setPcMinSat(phaseIdx, 0.0); + matParams.setPcMaxSat(phaseIdx, 0.0); + } + matParams.finalize(); + + //////////// + // Create a fluid state + //////////// + FluidState gasFluidState; + createSurfaceGasFluidSystem(gasFluidState); + + FluidState fluidState; + ParameterCache paramCache; + + // temperature + fluidState.setTemperature(T); + + // oil pressure + fluidState.setPressure(oilPhaseIdx, 4000 * 6894.7573); // 4000 PSI + + // oil saturation + fluidState.setSaturation(oilPhaseIdx, 1.0); + fluidState.setSaturation(gasPhaseIdx, 1.0 - fluidState.saturation(oilPhaseIdx)); + + // oil composition: SPE-5 reservoir oil + fluidState.setMoleFraction(oilPhaseIdx, H2OIdx, 0.0); + fluidState.setMoleFraction(oilPhaseIdx, C1Idx, 0.50); + fluidState.setMoleFraction(oilPhaseIdx, C3Idx, 0.03); + fluidState.setMoleFraction(oilPhaseIdx, C6Idx, 0.07); + fluidState.setMoleFraction(oilPhaseIdx, C10Idx, 0.20); + fluidState.setMoleFraction(oilPhaseIdx, C15Idx, 0.15); + fluidState.setMoleFraction(oilPhaseIdx, C20Idx, 0.05); + + //makeOilSaturated(fluidState, gasFluidState); + + // set the saturations and pressures of the other phases + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (phaseIdx != oilPhaseIdx) { + fluidState.setSaturation(phaseIdx, 0.0); + fluidState.setPressure(phaseIdx, fluidState.pressure(oilPhaseIdx)); + } + + // initial guess for the composition (needed by the ComputeFromReferencePhase + // constraint solver. TODO: bug in ComputeFromReferencePhase?) + guessInitial(fluidState, phaseIdx); + } + + typedef Opm::ComputeFromReferencePhase CFRP; + CFRP::solve(fluidState, + paramCache, + /*refPhaseIdx=*/oilPhaseIdx, + /*setViscosity=*/false, + /*setEnthalpy=*/false); + + //////////// + // Calculate the total molarities of the components + //////////// + ComponentVector totalMolarities; + for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) + totalMolarities[compIdx] = fluidState.saturation(oilPhaseIdx)*fluidState.molarity(oilPhaseIdx, compIdx); + + //////////// + // Gradually increase the volume for and calculate the gas + // formation factor, oil formation volume factor and gas formation + // volume factor. + //////////// + + FluidState flashFluidState, surfaceFluidState; + flashFluidState.assign(fluidState); + //Flash::guessInitial(flashFluidState, totalMolarities); + + using E = Opm::DenseAd::Evaluation; + using Flash = Opm::ChiFlash; + Flash::solve(flashFluidState, {0.9, 0.1}, 123456, 5, "SSI", 1e-8); + +} + +int main(int argc, char **argv) +{ + Dune::MPIHelper::instance(argc, argv); + + testAll(); + + // the Peng-Robinson test currently does not work with single-precision floating + // point scalars because of precision issues. (these are caused by the fact that the + // test uses finite differences to calculate derivatives.) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wunreachable-code" + while (0) testAll(); +#pragma GCC diagnostic pop + + + testChiFlash(); + + return 0; +}