[WIP] Add ChiFlash and test. Not compiling.
This commit is contained in:
Atgeirr Flø Rasmussen
Trine Mykkeltvedt
parent
fc98447a97
commit
7fa4f4796f
966
opm/material/constraintsolvers/ChiFlash.hpp
Normal file
966
opm/material/constraintsolvers/ChiFlash.hpp
Normal file
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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 <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
|
||||
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
|
||||
#include <opm/material/densead/Evaluation.hpp>
|
||||
#include <opm/material/densead/Math.hpp>
|
||||
#include <opm/material/common/MathToolbox.hpp>
|
||||
#include <opm/material/common/Valgrind.hpp>
|
||||
#include <opm/material/Constants.hpp>
|
||||
#include <opm/material/eos/PengRobinsonMixture.hpp>
|
||||
|
||||
#include <opm/material/common/Exceptions.hpp>
|
||||
|
||||
#include <dune/common/fvector.hh>
|
||||
#include <dune/common/fmatrix.hh>
|
||||
#include <dune/common/version.hh>
|
||||
|
||||
#include <limits>
|
||||
#include <iostream>
|
||||
#include <iomanip>
|
||||
|
||||
namespace Opm {
|
||||
|
||||
/*!
|
||||
* \brief Determines the phase compositions, pressures and saturations
|
||||
* given the total mass of all components for the chiwoms problem.
|
||||
*
|
||||
*/
|
||||
template <class Scalar, class Evaluation, class FluidSystem>
|
||||
class ChiFlash
|
||||
{
|
||||
//using Problem = GetPropType<TypeTag, Properties::Problem>;
|
||||
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 <class FluidState>
|
||||
static void solve(FluidState& fluidState,
|
||||
const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalComposition,
|
||||
int spatialIdx,
|
||||
int verbosity,
|
||||
std::string twoPhaseMethod,
|
||||
Scalar tolerance)
|
||||
{
|
||||
|
||||
using InputEval = typename FluidState::Scalar;
|
||||
using Matrix = Dune::FieldMatrix<InputEval, numEq, numEq>;
|
||||
using Vetor = Dune::FieldVector<InputEval, numEq>;
|
||||
using FlashEval = Opm::DenseAd::Evaluation</*Scalar=*/InputEval, /*numDerivs=*/numEq>;
|
||||
using FlashDefectVector = Dune::FieldVector<FlashEval, numEq>;
|
||||
using FlashFluidState = Opm::CompositionalFluidState<FlashEval, FluidSystem, /*energy=*/false>;
|
||||
|
||||
using ComponentVector = Dune::FieldVector<typename FluidState::Scalar, numComponents>;
|
||||
|
||||
#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
|
||||
Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
|
||||
#endif
|
||||
|
||||
if (tolerance <= 0)
|
||||
tolerance = std::min<Scalar>(1e-3, 1e8*std::numeric_limits<Scalar>::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<Evaluation> paramCache;
|
||||
paramCache.updatePhase(fluidState, oilPhaseIdx);
|
||||
paramCache.updatePhase(fluidState, gasPhaseIdx);
|
||||
|
||||
// Calculate compressibility factor
|
||||
const Scalar R = Opm::Constants<Scalar>::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 <<std::endl;
|
||||
std::cout << "Sg = " << Sg <<std::endl;
|
||||
std::cout << "Z_L = " << Z_L <<std::endl;
|
||||
std::cout << "Z_V = " << Z_V <<std::endl;
|
||||
}
|
||||
|
||||
// Update densities
|
||||
fluidState.setDensity(oilPhaseIdx, FluidSystem::density(fluidState, paramCache, oilPhaseIdx));
|
||||
fluidState.setDensity(gasPhaseIdx, FluidSystem::density(fluidState, paramCache, gasPhaseIdx));
|
||||
}//end solve
|
||||
|
||||
/*!
|
||||
* \brief Calculates the chemical equilibrium from the component
|
||||
* fugacities in a phase.
|
||||
*
|
||||
* This is a convenience method which assumes that the capillary pressure is
|
||||
* zero...
|
||||
*/
|
||||
template <class FluidState, class ComponentVector>
|
||||
static void solve(FluidState& fluidState,
|
||||
const ComponentVector& globalMolarities,
|
||||
Scalar tolerance = 0.0)
|
||||
{
|
||||
using MaterialTraits = NullMaterialTraits<Scalar, numPhases>;
|
||||
using MaterialLaw = NullMaterial<MaterialTraits>;
|
||||
using MaterialLawParams = typename MaterialLaw::Params;
|
||||
|
||||
MaterialLawParams matParams;
|
||||
solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
|
||||
}
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
template <class FlashFluidState>
|
||||
static typename FlashFluidState::Scalar wilsonK_(const FlashFluidState& fluidState, int compIdx)
|
||||
{
|
||||
using FlashEval = typename FlashFluidState::Scalar;
|
||||
const auto& acf = FluidSystem::acentricFactor(compIdx);
|
||||
const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
|
||||
const auto& T = fluidState.temperature(0);
|
||||
const auto& p_crit = FluidSystem::criticalPressure(compIdx);
|
||||
const auto& p = fluidState.pressure(0); //for now assume no capillary pressure
|
||||
|
||||
const auto& tmp = Opm::exp(5.3727 * (1+acf) * (1-T_crit/T)) * (p_crit/p);
|
||||
return tmp;
|
||||
}
|
||||
|
||||
template <class Vector, class FlashFluidState>
|
||||
static typename Vector::field_type li_single_phase_label_(const FlashFluidState& fluidState, const Vector& globalComposition, int verbosity)
|
||||
{
|
||||
// Calculate intermediate sum
|
||||
typename Vector::field_type sumVz = 0.0;
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
// Get component information
|
||||
const auto& V_crit = FluidSystem::criticalVolume(compIdx);
|
||||
|
||||
// Sum calculation
|
||||
sumVz += (V_crit * globalComposition[compIdx]);
|
||||
}
|
||||
|
||||
// Calculate approximate (pseudo) critical temperature using Li's method
|
||||
typename Vector::field_type Tc_est = 0.0;
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
// Get component information
|
||||
const auto& V_crit = FluidSystem::criticalVolume(compIdx);
|
||||
const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
|
||||
|
||||
// Sum calculation
|
||||
Tc_est += (V_crit * T_crit * globalComposition[compIdx] / sumVz);
|
||||
}
|
||||
|
||||
// Get temperature
|
||||
const auto& T = fluidState.temperature(0);
|
||||
|
||||
// If temperature is below estimated critical temperature --> phase = liquid; else vapor
|
||||
Evaluation L;
|
||||
if (T < Tc_est) {
|
||||
// Liquid
|
||||
L = 1.0;
|
||||
|
||||
// Print
|
||||
if (verbosity >= 1) {
|
||||
std::cout << "Cell is single-phase, liquid (L = 1.0) due to Li's phase labeling method giving T < Tc_est (" << T << " < " << Tc_est << ")!" << std::endl;
|
||||
}
|
||||
}
|
||||
else {
|
||||
// Vapor
|
||||
L = 0.0;
|
||||
|
||||
// Print
|
||||
if (verbosity >= 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 <class Vector>
|
||||
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<numComponents; ++compIdx){
|
||||
g += (globalComposition[compIdx]*(K[compIdx]-1))/(K[compIdx]-L*(K[compIdx]-1));
|
||||
}
|
||||
return g;
|
||||
}
|
||||
|
||||
template <class Vector>
|
||||
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<numComponents; ++compIdx){
|
||||
dg += (globalComposition[compIdx]*(K[compIdx]-1)*(K[compIdx]-1))/((K[compIdx]-L*(K[compIdx]-1))*(K[compIdx]-L*(K[compIdx]-1)));
|
||||
}
|
||||
return dg;
|
||||
}
|
||||
|
||||
template <class Vector>
|
||||
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<numComponents; ++compIdx){
|
||||
if (K[compIdx] < Kmin)
|
||||
Kmin = K[compIdx];
|
||||
else if (K[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 <class Vector>
|
||||
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 <class FlashFluidState, class ComponentVector>
|
||||
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<numComponents; ++compIdx){
|
||||
fluidState.setMoleFraction(gasPhaseIdx, compIdx, globalComposition[compIdx]);
|
||||
fluidState.setMoleFraction(oilPhaseIdx, compIdx, globalComposition[compIdx]);
|
||||
}
|
||||
}
|
||||
// If not stable: use the mole fractions from Michelsen's test to update K
|
||||
else {
|
||||
for (int compIdx = 0; compIdx<numComponents; ++compIdx) {
|
||||
K[compIdx] = y[compIdx] / x[compIdx];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template <class FlashFluidState, class ComponentVector>
|
||||
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<Scalar, FluidSystem>;
|
||||
|
||||
// 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<numComponents; ++compIdx){
|
||||
xy_loc[compIdx] = K[compIdx] * globalComposition[compIdx];
|
||||
S_loc += xy_loc[compIdx];
|
||||
}
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
xy_loc[compIdx] /= S_loc;
|
||||
fluidState_fake.setMoleFraction(gasPhaseIdx, compIdx, xy_loc[compIdx]);
|
||||
}
|
||||
}
|
||||
else {
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
xy_loc[compIdx] = globalComposition[compIdx]/K[compIdx];
|
||||
S_loc += xy_loc[compIdx];
|
||||
}
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
xy_loc[compIdx] /= S_loc;
|
||||
fluidState_fake.setMoleFraction(oilPhaseIdx, compIdx, xy_loc[compIdx]);
|
||||
}
|
||||
}
|
||||
|
||||
int phaseIdx = (isGas?gasPhaseIdx:oilPhaseIdx);
|
||||
int phaseIdx2 = (isGas?oilPhaseIdx:gasPhaseIdx);
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
fluidState_global.setMoleFraction(phaseIdx2, compIdx, globalComposition[compIdx]);
|
||||
}
|
||||
|
||||
typename FluidSystem::template ParameterCache<FlashEval> paramCache_fake;
|
||||
paramCache_fake.updatePhase(fluidState_fake, phaseIdx);
|
||||
|
||||
typename FluidSystem::template ParameterCache<FlashEval> paramCache_global;
|
||||
paramCache_global.updatePhase(fluidState_global, phaseIdx2);
|
||||
|
||||
//fugacity for fake phases each component
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
auto phiFake = PengRobinsonMixture::computeFugacityCoefficient(fluidState_fake, paramCache_fake, phaseIdx, compIdx);
|
||||
auto phiGlobal = PengRobinsonMixture::computeFugacityCoefficient(fluidState_global, paramCache_global, phaseIdx2, compIdx);
|
||||
|
||||
fluidState_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake);
|
||||
fluidState_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
|
||||
}
|
||||
|
||||
|
||||
ComponentVector R;
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
if (isGas){
|
||||
auto fug_fake = fluidState_fake.fugacity(phaseIdx, compIdx);
|
||||
auto fug_global = fluidState_global.fugacity(phaseIdx2, compIdx);
|
||||
auto fug_ratio = fug_global / fug_fake;
|
||||
R[compIdx] = fug_ratio / S_loc;
|
||||
}
|
||||
else{
|
||||
auto fug_fake = fluidState_fake.fugacity(phaseIdx, compIdx);
|
||||
auto fug_global = fluidState_global.fugacity(phaseIdx2, compIdx);
|
||||
auto fug_ratio = fug_fake / fug_global;
|
||||
R[compIdx] = fug_ratio * S_loc;
|
||||
}
|
||||
}
|
||||
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
K[compIdx] *= R[compIdx];
|
||||
}
|
||||
Scalar R_norm = 0.0;
|
||||
Scalar K_norm = 0.0;
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
auto a = Opm::getValue(R[compIdx]) - 1.0;
|
||||
auto b = Opm::log(Opm::getValue(K[compIdx]));
|
||||
R_norm += a*a;
|
||||
K_norm += b*b;
|
||||
}
|
||||
|
||||
// Print iteration info
|
||||
if (verbosity == 3 || verbosity == 4) {
|
||||
std::cout << std::setw(10) << i << std::setw(16) << K_norm << std::setw(16) << R_norm << std::endl;
|
||||
}
|
||||
|
||||
// Check convergence
|
||||
isTrivial = (K_norm < 1e-5);
|
||||
if (isTrivial || R_norm < 1e-10)
|
||||
return;
|
||||
//todo: make sure that no molefraction is smaller than 1e-8 ?
|
||||
//todo: take care of water!
|
||||
}
|
||||
throw std::runtime_error(" Stability test did not converge");
|
||||
}//end checkStability
|
||||
|
||||
template <class FlashFluidState, class ComponentVector>
|
||||
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<numComponents; ++compIdx){
|
||||
x[compIdx] = globalComposition[compIdx]/(L + (1-L)*K[compIdx]);
|
||||
sumx += x[compIdx];
|
||||
y[compIdx] = (K[compIdx]*globalComposition[compIdx])/(L + (1-L)*K[compIdx]);
|
||||
sumy += y[compIdx];
|
||||
}
|
||||
x /= sumx;
|
||||
y /= sumy;
|
||||
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
fluidState.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]);
|
||||
fluidState.setMoleFraction(gasPhaseIdx, compIdx, y[compIdx]);
|
||||
//SET k ??
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
template <class FlashFluidState, class ComponentVector>
|
||||
static void newtonCompositionUpdate_(ComponentVector& K, Evaluation& L, FlashFluidState& fluidState, const ComponentVector& globalComposition,
|
||||
int verbosity)
|
||||
{
|
||||
// Newton declarations
|
||||
using NewtonVector = Dune::FieldVector<Evaluation, numMisciblePhases*numMiscibleComponents+1>;
|
||||
using NewtonMatrix = Dune::FieldMatrix<Evaluation, numMisciblePhases*numMiscibleComponents+1, numMisciblePhases*numMiscibleComponents+1>;
|
||||
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<numComponents; ++compIdx){
|
||||
if (compIdx < numComponents - 1)
|
||||
std::cout << fluidState.moleFraction(oilPhaseIdx, compIdx) << " ";
|
||||
else
|
||||
std::cout << fluidState.moleFraction(oilPhaseIdx, compIdx);
|
||||
}
|
||||
std::cout << "], y = [";
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
if (compIdx < numComponents - 1)
|
||||
std::cout << fluidState.moleFraction(gasPhaseIdx, compIdx) << " ";
|
||||
else
|
||||
std::cout << fluidState.moleFraction(gasPhaseIdx, compIdx);
|
||||
}
|
||||
std::cout << "], and " << "L = " << L << std::endl;
|
||||
}
|
||||
|
||||
// Print header
|
||||
if (verbosity == 2 || verbosity == 4) {
|
||||
std::cout << std::setw(10) << "Iteration" << std::setw(16) << "Norm2(step)" << std::setw(16) << "Norm2(Residual)" << std::endl;
|
||||
}
|
||||
|
||||
// Assign primary variables (x, y and L)
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
newtonX[compIdx] = Opm::getValue(fluidState.moleFraction(oilPhaseIdx, compIdx));
|
||||
newtonX[compIdx + numMiscibleComponents] = Opm::getValue(fluidState.moleFraction(gasPhaseIdx, compIdx));
|
||||
}
|
||||
newtonX[numMiscibleComponents*numMiscibleComponents] = Opm::getValue(L);
|
||||
|
||||
//
|
||||
// Main Newton loop
|
||||
//
|
||||
bool convFug = false; /* for convergence check */
|
||||
for (int i = 0; i< 100; ++i){
|
||||
// Evaluate residuals (newtonB)
|
||||
evalDefect_(newtonB, newtonX, fluidState, globalComposition);
|
||||
|
||||
// Print iteration info
|
||||
if (verbosity == 2 || verbosity == 4) {
|
||||
if (i == 0) {
|
||||
std::cout << std::setw(10) << i << std::setw(16) << "N/A" << std::setw(16) << newtonB.two_norm() << std::endl;
|
||||
}
|
||||
else {
|
||||
std::cout << std::setw(10) << i << std::setw(16) << newtonDelta.two_norm() << std::setw(16) << newtonB.two_norm() << std::endl;
|
||||
}
|
||||
}
|
||||
|
||||
// Check fugacity equilibrium for convergence
|
||||
convFug = checkFugacityEquil_(newtonB);
|
||||
|
||||
// If convergence have been met, we abort; else we update step and loop once more
|
||||
if (convFug == true) {
|
||||
// Extract x, y and L together with calculation of K
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
fluidState.setMoleFraction(oilPhaseIdx, compIdx, newtonX[compIdx]);
|
||||
fluidState.setMoleFraction(gasPhaseIdx, compIdx, newtonX[compIdx + numMiscibleComponents]);
|
||||
K[compIdx] = newtonX[compIdx + numMiscibleComponents] / newtonX[compIdx];
|
||||
}
|
||||
L = newtonX[numMiscibleComponents*numMiscibleComponents];
|
||||
|
||||
// Print info
|
||||
if (verbosity >= 1) {
|
||||
std::cout << "Solution converged to the following result :" << std::endl;
|
||||
std::cout << "x = [";
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
if (compIdx < numComponents - 1)
|
||||
std::cout << newtonX[compIdx] << " ";
|
||||
else
|
||||
std::cout << newtonX[compIdx];
|
||||
}
|
||||
std::cout << "]" << std::endl;
|
||||
std::cout << "y = [";
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
if (compIdx < numComponents - 1)
|
||||
std::cout << newtonX[compIdx + numMiscibleComponents] << " ";
|
||||
else
|
||||
std::cout << newtonX[compIdx + numMiscibleComponents];
|
||||
}
|
||||
std::cout << "]" << std::endl;
|
||||
std::cout << "K = [" << K << "]" << std::endl;
|
||||
std::cout << "L = " << L << std::endl;
|
||||
}
|
||||
return;
|
||||
}
|
||||
else {
|
||||
// Calculate Jacobian (newtonA)
|
||||
evalJacobian_(newtonA, newtonX, fluidState, globalComposition);
|
||||
|
||||
// Solve system J * x = -r, which in our case is newtonA*newtonX = newtonB, to get next step (newtonDelta)
|
||||
newtonA.solve(newtonDelta, newtonB);
|
||||
|
||||
// Update current solution (newtonX) with simple relaxation method (percentage of step applied)
|
||||
updateCurrentSol_(newtonX, newtonDelta);
|
||||
}
|
||||
}
|
||||
throw std::runtime_error(" Newton composition update did not converge within maxIterations");
|
||||
}
|
||||
|
||||
template <class DefectVector>
|
||||
static void updateCurrentSol_(DefectVector& x, DefectVector& d)
|
||||
{
|
||||
// Find smallest percentage update
|
||||
Scalar w = 1.0;
|
||||
for (int i=0; i<x.size(); ++i){
|
||||
Scalar w_tmp = Opm::getValue(Opm::min(Opm::max(x[i] + d[i], 0.0), 1.0) - x[i]) / Opm::getValue(d[i]);
|
||||
w = Opm::min(w, w_tmp);
|
||||
}
|
||||
|
||||
// Loop over the solution vector and apply the smallest percentage update
|
||||
for (int i=0; i<x.size(); ++i){
|
||||
x[i] += w*d[i];
|
||||
}
|
||||
}
|
||||
|
||||
template <class DefectVector>
|
||||
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<numComponents; ++compIdx){
|
||||
fugVec[compIdx] = 0.0;
|
||||
fugVec[compIdx+numMiscibleComponents] = b[compIdx+numMiscibleComponents];
|
||||
}
|
||||
fugVec[numMiscibleComponents*numMisciblePhases] = 0.0;
|
||||
|
||||
// Check if norm(fugVec) is less than tolerance
|
||||
bool conv = fugVec.two_norm() < 1e-6;
|
||||
return conv;
|
||||
}
|
||||
|
||||
template <class FluidState, class DefectVector, class ComponentVector>
|
||||
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<numComponents; ++compIdx){
|
||||
fluidState.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]);
|
||||
fluidState.setMoleFraction(gasPhaseIdx, compIdx, x[compIdx + numMiscibleComponents]);
|
||||
}
|
||||
|
||||
|
||||
// Compute fugacities
|
||||
using ParamCache = typename FluidSystem::template ParameterCache<typename FluidState::Scalar>;
|
||||
ParamCache paramCache;
|
||||
for (int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
|
||||
paramCache.updatePhase(fluidState, phaseIdx);
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
Evaluation phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
|
||||
fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
|
||||
}
|
||||
}
|
||||
|
||||
// Compute residuals for Newton update. Primary variables are: x, y, and L
|
||||
// TODO: Make this AD
|
||||
// Extract L
|
||||
Evaluation L = x[numMiscibleComponents*numMisciblePhases];
|
||||
|
||||
// Residuals
|
||||
// OBS: the residuals are negative in the newton system!
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
// z - L*x - (1-L) * y
|
||||
b[compIdx] = -globalComposition[compIdx] + L*x[compIdx] + (1-L)*x[compIdx + numMiscibleComponents];
|
||||
|
||||
// (f_liquid/f_vapor) - 1 = 0
|
||||
b[compIdx + numMiscibleComponents] = -(fluidState.fugacity(oilPhaseIdx, compIdx) / fluidState.fugacity(gasPhaseIdx, compIdx)) + 1.0;
|
||||
|
||||
// sum(x) - sum(y) = 0
|
||||
b[numMiscibleComponents*numMisciblePhases] += -x[compIdx] + x[compIdx + numMiscibleComponents];
|
||||
}
|
||||
}//end valDefect
|
||||
|
||||
template <class FluidState, class DefectVector, class DefectMatrix, class ComponentVector>
|
||||
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<xIn.size(); ++j){
|
||||
x[j] = xIn[j];
|
||||
}
|
||||
|
||||
evalDefect_(b0, x, fluidStateIn, globalComposition);
|
||||
|
||||
// Make the jacobian A in Newton system Ax=b
|
||||
Scalar epsilon = 1e-10;
|
||||
for(int i=0; i<b0.size(); ++i){
|
||||
// Permutate x and calculate response
|
||||
x[i] += epsilon;
|
||||
DefectVector bEps;
|
||||
evalDefect_(bEps, x, fluidStateIn, globalComposition);
|
||||
x[i] -= epsilon;
|
||||
|
||||
// Forward difference of all eqs wrt primary variable i
|
||||
DefectVector derivI;
|
||||
for(int j=0; j<b0.size(); ++j){
|
||||
derivI[j] = bEps[j];
|
||||
derivI[j] -= b0[j];
|
||||
derivI[j] /= epsilon;
|
||||
A[j][i] = derivI[j];
|
||||
}
|
||||
}
|
||||
}//end evalJacobian
|
||||
|
||||
template <class FlashFluidState, class ComponentVector>
|
||||
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<typename FlashFluidState::Scalar>;
|
||||
ParamCache paramCache;
|
||||
for (int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
|
||||
paramCache.updatePhase(fluidState, phaseIdx);
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
Evaluation phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
|
||||
fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
|
||||
}
|
||||
}
|
||||
|
||||
// Calculate fugacity ratio
|
||||
ComponentVector newFugRatio;
|
||||
ComponentVector convFugRatio;
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
newFugRatio[compIdx] = fluidState.fugacity(oilPhaseIdx, compIdx)/fluidState.fugacity(gasPhaseIdx, compIdx);
|
||||
convFugRatio[compIdx] = newFugRatio[compIdx] - 1.0;
|
||||
}
|
||||
|
||||
// Print iteration info
|
||||
if (verbosity == 2 || verbosity == 4) {
|
||||
int prec = 5;
|
||||
int fugWidth = (prec + 3);
|
||||
int convWidth = prec + 9;
|
||||
std::cout << std::defaultfloat;
|
||||
std::cout << std::fixed;
|
||||
std::cout << std::setw(5) << i;
|
||||
std::cout << std::setw(fugWidth);
|
||||
std::cout << std::setprecision(prec);
|
||||
std::cout << newFugRatio;
|
||||
std::cout << std::scientific;
|
||||
std::cout << std::setw(convWidth) << convFugRatio.two_norm() << std::endl;
|
||||
}
|
||||
|
||||
// Check convergence
|
||||
if (convFugRatio.two_norm() < 1e-6){
|
||||
// Restore cout format
|
||||
std::cout.flags(f);
|
||||
|
||||
// Print info
|
||||
if (verbosity >= 1) {
|
||||
std::cout << "Solution converged to the following result :" << std::endl;
|
||||
std::cout << "x = [";
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
if (compIdx < numComponents - 1)
|
||||
std::cout << fluidState.moleFraction(oilPhaseIdx, compIdx) << " ";
|
||||
else
|
||||
std::cout << fluidState.moleFraction(oilPhaseIdx, compIdx);
|
||||
}
|
||||
std::cout << "]" << std::endl;
|
||||
std::cout << "y = [";
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
if (compIdx < numComponents - 1)
|
||||
std::cout << fluidState.moleFraction(gasPhaseIdx, compIdx) << " ";
|
||||
else
|
||||
std::cout << fluidState.moleFraction(gasPhaseIdx, compIdx);
|
||||
}
|
||||
std::cout << "]" << std::endl;
|
||||
std::cout << "K = [" << K << "]" << std::endl;
|
||||
std::cout << "L = " << L << std::endl;
|
||||
}
|
||||
// Restore cout format format
|
||||
return;
|
||||
}
|
||||
|
||||
// If convergence is not met, K is updated in a successive substitution manner
|
||||
else{
|
||||
// Update K
|
||||
for (int compIdx=0; compIdx<numComponents; ++compIdx){
|
||||
K[compIdx] *= newFugRatio[compIdx];
|
||||
}
|
||||
|
||||
// Solve Rachford-Rice to get L from updated K
|
||||
L = solveRachfordRice_g_(K, globalComposition, 0);
|
||||
}
|
||||
|
||||
}
|
||||
throw std::runtime_error("Successive substitution composition update did not converge within maxIterations");
|
||||
}
|
||||
|
||||
};//end ChiFlash
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
#endif
|
||||
621
tests/test_chiflash.cpp
Normal file
621
tests/test_chiflash.cpp
Normal file
@ -0,0 +1,621 @@
|
||||
// -*- 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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
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 <opm/material/constraintsolvers/ChiFlash.hpp>
|
||||
|
||||
#include <opm/material/densead/Evaluation.hpp>
|
||||
#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
|
||||
#include <opm/material/constraintsolvers/NcpFlash.hpp>
|
||||
#include <opm/material/fluidstates/CompositionalFluidState.hpp>
|
||||
#include <opm/material/fluidsystems/Spe5FluidSystem.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
|
||||
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
|
||||
|
||||
#include <dune/common/parallel/mpihelper.hh>
|
||||
|
||||
template <class FluidSystem, class FluidState>
|
||||
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<typename FluidState::Scalar> paramCache;
|
||||
paramCache.updatePhase(gasFluidState, gasPhaseIdx);
|
||||
gasFluidState.setDensity(gasPhaseIdx,
|
||||
FluidSystem::density(gasFluidState, paramCache, gasPhaseIdx));
|
||||
}
|
||||
|
||||
template <class Scalar, class FluidSystem, class FluidState>
|
||||
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<Scalar, numComponents> ComponentVector;
|
||||
typedef Opm::NcpFlash<Scalar, FluidSystem> 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 <class Scalar, class FluidSystem, class FluidState>
|
||||
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<Scalar, FluidSystem>(prestineFluidState,
|
||||
fluidState,
|
||||
gasFluidState,
|
||||
additionalGas);
|
||||
Scalar fStar = 1 - computeSumxg<Scalar, FluidSystem>(prestineFluidState,
|
||||
fluidState,
|
||||
gasFluidState,
|
||||
additionalGas + eps);
|
||||
Scalar fPrime = (fStar - f)/eps;
|
||||
|
||||
additionalGas -= f/fPrime;
|
||||
};
|
||||
}
|
||||
|
||||
template <class FluidSystem, class FluidState>
|
||||
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 <class Scalar, class FluidSystem, class FluidState>
|
||||
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<Scalar, FluidSystem> Flash;
|
||||
typedef Opm::ThreePhaseMaterialTraits<Scalar, waterPhaseIdx, oilPhaseIdx, gasPhaseIdx> MaterialTraits;
|
||||
typedef Opm::LinearMaterial<MaterialTraits> MaterialLaw;
|
||||
typedef typename MaterialLaw::Params MaterialLawParams;
|
||||
typedef Dune::FieldVector<Scalar, numComponents> 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<Scalar> 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<MaterialLaw>(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<MaterialLaw>(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<MaterialLaw>(surfaceFluidState, matParams, paramCache, tmpMolarities);
|
||||
return alpha;
|
||||
}
|
||||
|
||||
template <class RawTable>
|
||||
void printResult(const RawTable& rawTable,
|
||||
const std::string& fieldName,
|
||||
size_t firstIdx,
|
||||
size_t secondIdx,
|
||||
double hiresThres)
|
||||
{
|
||||
std::cout << "std::vector<std::pair<Scalar, Scalar> > "<<fieldName<<" = {\n";
|
||||
|
||||
size_t sampleIdx = 0;
|
||||
size_t numSamples = 20;
|
||||
size_t numRawHires = 0;
|
||||
for (; rawTable[numRawHires][firstIdx] > 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 <class Scalar>
|
||||
inline void testAll()
|
||||
{
|
||||
typedef Opm::Spe5FluidSystem<Scalar> 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<Scalar, FluidSystem> Flash;
|
||||
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
|
||||
typedef Opm::CompositionalFluidState<Scalar, FluidSystem> FluidState;
|
||||
|
||||
typedef Opm::ThreePhaseMaterialTraits<Scalar, waterPhaseIdx, oilPhaseIdx, gasPhaseIdx> MaterialTraits;
|
||||
typedef Opm::LinearMaterial<MaterialTraits> MaterialLaw;
|
||||
typedef typename MaterialLaw::Params MaterialLawParams;
|
||||
|
||||
typedef typename FluidSystem::template ParameterCache<Scalar> 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<FluidSystem>(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<Scalar, FluidSystem>(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<FluidSystem>(fluidState, phaseIdx);
|
||||
}
|
||||
|
||||
typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> 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<MaterialLaw>(flashFluidState, matParams, paramCache, totalMolarities);
|
||||
|
||||
Scalar surfaceAlpha = 1;
|
||||
surfaceAlpha = bringOilToSurface<Scalar, FluidSystem>(surfaceFluidState, surfaceAlpha, flashFluidState, /*guessInitial=*/true);
|
||||
Scalar rho_gRef = surfaceFluidState.density(gasPhaseIdx);
|
||||
Scalar rho_oRef = surfaceFluidState.density(oilPhaseIdx);
|
||||
|
||||
std::vector<std::array<Scalar, 10> > 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] <M_o>[kg/mol] <M_g>[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<MaterialLaw>(flashFluidState, matParams, paramCache, curTotalMolarities);
|
||||
|
||||
surfaceAlpha = bringOilToSurface<Scalar, FluidSystem>(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<Scalar, 10> 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<Scalar> 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<Scalar, FluidSystem> Flash;
|
||||
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
|
||||
typedef Opm::CompositionalFluidState<Scalar, FluidSystem> FluidState;
|
||||
|
||||
typedef Opm::ThreePhaseMaterialTraits<Scalar, waterPhaseIdx, oilPhaseIdx, gasPhaseIdx> MaterialTraits;
|
||||
typedef Opm::LinearMaterial<MaterialTraits> MaterialLaw;
|
||||
typedef typename MaterialLaw::Params MaterialLawParams;
|
||||
|
||||
typedef typename FluidSystem::template ParameterCache<Scalar> 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<FluidSystem>(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<Scalar, FluidSystem>(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<FluidSystem>(fluidState, phaseIdx);
|
||||
}
|
||||
|
||||
typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> 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<double, 3>;
|
||||
using Flash = Opm::ChiFlash<double, E, FluidSystem>;
|
||||
Flash::solve(flashFluidState, {0.9, 0.1}, 123456, 5, "SSI", 1e-8);
|
||||
|
||||
}
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
Dune::MPIHelper::instance(argc, argv);
|
||||
|
||||
testAll<double>();
|
||||
|
||||
// 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<float>();
|
||||
#pragma GCC diagnostic pop
|
||||
|
||||
|
||||
testChiFlash();
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user