OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

wip co2brine

Browse Source
This commit is contained in:
Trine Mykkeltvedt 2022-06-13 14:11:08 +02:00
parent d26ea55ef8
commit 4fd4fc7029
3 changed files with 282 additions and 580 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

657
opm/material/constraintsolvers/ChiFlash.hpp
View File

@ -146,17 +146,12 @@ public:
fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0))); fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0)));
// Rachford Rice equation to get initial L for composition solver
//L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity);
// Do a stability test to check if cell is single-phase (do for all cells the first time). // Do a stability test to check if cell is single-phase (do for all cells the first time).
bool isStable = false; bool isStable = false;
if ( L <= 0 || L == 1 ) { if ( L <= 0 || L == 1 ) {
if (verbosity >= 1) { if (verbosity >= 1) {
std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl; std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl;
} }
//phaseStabilityTestMichelsen_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
phaseStabilityTest_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity); phaseStabilityTest_(isStable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
} }
@ -210,123 +205,14 @@ public:
// fluid_state.setLvalue(L_scalar); //print summary after flash
if (verbosity >= 1) {
// std::cout << " ------ SUMMARY AFTER DERIVATIVES ------ " << std::endl; std::cout << " ------ SUMMARY AFTER FLASH ------ " << std::endl;
// std::cout << " L " << fluid_state.L() << std::endl; std::cout << " L " << fluid_state.L() << std::endl;
// std::cout << " K " << fluid_state.K(0) << ", " << fluid_state.K(1) << ", " << fluid_state.K(2) << std::endl; std::cout << " K " << fluid_state.K(0) << ", " << fluid_state.K(1) << ", " << fluid_state.K(2) << std::endl;
// std::cout << " x " << fluid_state.moleFraction(oilPhaseIdx, 0) << ", " << fluid_state.moleFraction(oilPhaseIdx, 1) << ", " << fluid_state.moleFraction(oilPhaseIdx, 2) << std::endl; std::cout << " x " << fluid_state.moleFraction(oilPhaseIdx, 0) << ", " << fluid_state.moleFraction(oilPhaseIdx, 1) << ", " << fluid_state.moleFraction(oilPhaseIdx, 2) << std::endl;
// std::cout << " y " << fluid_state.moleFraction(gasPhaseIdx, 0) << ", " << fluid_state.moleFraction(gasPhaseIdx, 1) << ", " << fluid_state.moleFraction(gasPhaseIdx, 2) << std::endl; std::cout << " y " << fluid_state.moleFraction(gasPhaseIdx, 0) << ", " << fluid_state.moleFraction(gasPhaseIdx, 1) << ", " << fluid_state.moleFraction(gasPhaseIdx, 2) << std::endl;
// Update phases
/* typename FluidSystem::template ParameterCache<InputEval> paramCache;
paramCache.updatePhase(fluid_state, oilPhaseIdx);
paramCache.updatePhase(fluid_state, gasPhaseIdx); */
/* // Calculate compressibility factor
const Scalar R = Opm::Constants<Scalar>::R;
InputEval Z_L = (paramCache.molarVolume(oilPhaseIdx) * fluid_state.pressure(oilPhaseIdx) ) /
(R * fluid_state.temperature(oilPhaseIdx));
InputEval Z_V = (paramCache.molarVolume(gasPhaseIdx) * fluid_state.pressure(gasPhaseIdx) ) /
(R * fluid_state.temperature(gasPhaseIdx));
std::cout << " the type of InputEval here is " << Dune::className<InputEval>() << std::endl;
// Update saturation
InputEval So = L*Z_L/(L*Z_L+(1-L)*Z_V);
InputEval Sg = 1-So;
fluid_state.setSaturation(oilPhaseIdx, So);
fluid_state.setSaturation(gasPhaseIdx, Sg);
//Update L and K to the problem for the next flash
for (int compIdx = 0; compIdx < numComponents; ++compIdx){
fluid_state.setKvalue(compIdx, K[compIdx]);
} }
fluid_state.setCompressFactor(oilPhaseIdx, Z_L);
fluid_state.setCompressFactor(gasPhaseIdx, Z_V);
fluid_state.setLvalue(L); */
// Print saturation
/* std::cout << " output the molefraction derivatives" << std::endl;
std::cout << " for vapor comp 1 " << std::endl;
fluid_state.moleFraction(gasPhaseIdx, 0).print();
std::cout << std::endl << " for vapor comp 2 " << std::endl;
fluid_state.moleFraction(gasPhaseIdx, 1).print();
std::cout << std::endl << " for vapor comp 3 " << std::endl;
fluid_state.moleFraction(gasPhaseIdx, 2).print();
std::cout << std::endl;
std::cout << " for oil comp 1 " << std::endl;
fluid_state.moleFraction(oilPhaseIdx, 0).print();
std::cout << std::endl << " for oil comp 2 " << std::endl;
fluid_state.moleFraction(oilPhaseIdx, 1).print();
std::cout << std::endl << " for oil comp 3 " << std::endl;
fluid_state.moleFraction(oilPhaseIdx, 2).print();
std::cout << std::endl;
std::cout << " for pressure " << std::endl;
fluid_state.pressure(0).print();
std::cout<< std::endl;
fluid_state.pressure(1).print();
std::cout<< std::endl; */
// Update densities
// fluid_state.setDensity(oilPhaseIdx, FluidSystem::density(fluid_state, paramCache, oilPhaseIdx));
// fluid_state.setDensity(gasPhaseIdx, FluidSystem::density(fluid_state, paramCache, gasPhaseIdx));
// check the residuals of the equations
/* using NewtonVector = Dune::FieldVector<InputEval, numMisciblePhases*numMiscibleComponents+1>;
NewtonVector newtonX;
NewtonVector newtonB;
for (int compIdx=0; compIdx<numComponents; ++compIdx){
newtonX[compIdx] = Opm::getValue(fluid_state.moleFraction(oilPhaseIdx, compIdx));
newtonX[compIdx + numMiscibleComponents] = Opm::getValue(fluid_state.moleFraction(gasPhaseIdx, compIdx));
}
newtonX[numMisciblePhases*numMiscibleComponents] = Opm::getValue(L);
evalDefect_(newtonB, newtonX, fluid_state, z); */
/* std::cout << " the residuals of the equations is " << std::endl;
for (unsigned i = 0; i < newtonB.N(); ++i) {
std::cout << newtonB[i] << std::endl;
}
std::cout << std::endl;
if (verbosity >= 5) {
std::cout << " mole fractions for oil " << std::endl;
for (int i = 0; i < numComponents; ++i) {
std::cout << " i " << i << " " << fluid_state.moleFraction(oilPhaseIdx, i) << std::endl;
}
std::cout << " mole fractions for gas " << std::endl;
for (int i = 0; i < numComponents; ++i) {
std::cout << " i " << i << " " << fluid_state.moleFraction(gasPhaseIdx, i) << std::endl;
}
std::cout << " K " << std::endl;
for (int i = 0; i < numComponents; ++i) {
std::cout << " i " << K[i] << std::endl;
}
std::cout << "Z_L = " << Z_L <<std::endl;
std::cout << "Z_V = " << Z_V <<std::endl;
std::cout << " fugacity for oil phase " << std::endl;
for (int i = 0; i < numComponents; ++i) {
auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, oilPhaseIdx, i);
std::cout << " i " << i << " " << phi << std::endl;
}
std::cout << " fugacity for gas phase " << std::endl;
for (int i = 0; i < numComponents; ++i) {
auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, gasPhaseIdx, i);
std::cout << " i " << i << " " << phi << std::endl;
}
std::cout << " density for oil phase " << std::endl;
std::cout << FluidSystem::density(fluid_state, paramCache, oilPhaseIdx) << std::endl;
std::cout << " density for gas phase " << std::endl;
std::cout << FluidSystem::density(fluid_state, paramCache, gasPhaseIdx) << std::endl;
std::cout << " viscosities for oil " << std::endl;
std::cout << FluidSystem::viscosity(fluid_state, paramCache, oilPhaseIdx) << std::endl;
std::cout << " viscosities for gas " << std::endl;
std::cout << FluidSystem::viscosity(fluid_state, paramCache, gasPhaseIdx) << std::endl;
std::cout << "So = " << So <<std::endl;
std::cout << "Sg = " << Sg <<std::endl;
} */
}//end solve }//end solve
/*! /*!
@ -564,185 +450,6 @@ protected:
throw std::runtime_error(" Rachford-Rice with bisection failed!"); throw std::runtime_error(" Rachford-Rice with bisection failed!");
} }
template <class FlashFluidState, class ComponentVector>
static void phaseStabilityTestMichelsen_(bool& stable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& z, int verbosity)
{
// Declarations
bool isTrivialL, isTrivialV;
ComponentVector x, y;
typename FlashFluidState::Scalar S_l, S_v;
ComponentVector K0 = K;
ComponentVector K1 = K;
// Check for vapour instable phase
if (verbosity == 3 || verbosity == 4) {
std::cout << "Stability test for vapor phase:" << std::endl;
}
bool stable_vapour = false;
michelsenTest_(fluidState, z, y, K0,stable_vapour,/*isGas*/true, verbosity);
bool stable_liquid = false;
michelsenTest_(fluidState, z, x, K1,stable_liquid,/*isGas*/false, verbosity);
//bool stable = false;
stable = stable_liquid && stable_vapour;
if (!stable) {
for (int compIdx = 0; compIdx<numComponents; ++compIdx) {
K[compIdx] = y[compIdx] / x[compIdx];
}
} else {
// Single phase, i.e. phase composition is equivalent to the global composition
// Update fluidstate with mole fraction
for (int compIdx=0; compIdx<numComponents; ++compIdx){
fluidState.setMoleFraction(gasPhaseIdx, compIdx, z[compIdx]);
fluidState.setMoleFraction(oilPhaseIdx, compIdx, z[compIdx]);
}
}
// printing
if (verbosity >= 1) {
std::cout << "Stability test done for - vapour - liquid - sum:" << stable_vapour << " - " << stable_liquid << " - " << stable <<std::endl;
}
}
template <class FlashFluidState, class ComponentVector>
static void michelsenTest_(const FlashFluidState& fluidState, const ComponentVector z, ComponentVector& xy_out,
ComponentVector& K, bool& stable, bool isGas, int verbosity)
{
using FlashEval = typename FlashFluidState::Scalar;
using PengRobinsonMixture = typename Opm::PengRobinsonMixture<Scalar, FluidSystem>;
// Declarations
FlashFluidState fluidState_xy = fluidState;
FlashFluidState fluidState_zi = fluidState;
ComponentVector xy_loc;
ComponentVector R;
FlashEval S_loc = 0.0;
FlashEval xy_c = 0.0;
bool isTrivial;
bool isConverged;
int phaseIdx = (isGas ? static_cast<int>(gasPhaseIdx) : static_cast<int>(oilPhaseIdx));
int phaseIdx2 = (isGas ? static_cast<int>(oilPhaseIdx) : static_cast<int>(gasPhaseIdx));
// Setup output
if (verbosity >= 3 || verbosity >= 4) {
std::cout << std::setw(10) << "Iteration" << std::setw(16) << "K-Norm" << std::setw(16) << "R-Norm" << std::endl;
}
//mixture fugacity
for (int compIdx=0; compIdx<numComponents; ++compIdx){
fluidState_zi.setMoleFraction(oilPhaseIdx, compIdx, z[compIdx]);
}
typename FluidSystem::template ParameterCache<FlashEval> paramCache_zi;
paramCache_zi.updatePhase(fluidState_zi, oilPhaseIdx);
for (int compIdx=0; compIdx<numComponents; ++compIdx){
auto phi_z = PengRobinsonMixture::computeFugacityCoefficient(fluidState_zi, paramCache_zi, oilPhaseIdx, compIdx);
fluidState_zi.setFugacityCoefficient(oilPhaseIdx, compIdx, phi_z);
auto f_zi = fluidState_zi.fugacity(oilPhaseIdx, compIdx);
//std::cout << "comp" << compIdx <<" , f_zi " << f_zi << std::endl;
}
// Michelsens stability test.
// Make two fake phases "inside" one phase and check for positive volume
int maxIter = 20000;
for (int i = 0; i < maxIter; ++i) {
S_loc = 0.0;
for (int compIdx=0; compIdx<numComponents; ++compIdx){
if (isGas) {
xy_c = K[compIdx] * z[compIdx];
} else {
xy_c = z[compIdx]/K[compIdx];
}
xy_loc[compIdx] = xy_c;
S_loc += xy_c;
}
xy_loc /= S_loc;
if (isGas)
xy_out = z;
else
xy_out = xy_loc;
//to get out fugacities each phase
for (int compIdx=0; compIdx<numComponents; ++compIdx){
if (isGas) {
fluidState_xy.setMoleFraction(gasPhaseIdx, compIdx, xy_loc[compIdx]);
} else {
fluidState_xy.setMoleFraction(oilPhaseIdx, compIdx, xy_loc[compIdx]);
}
}
typename FluidSystem::template ParameterCache<FlashEval> paramCache_xy;
paramCache_xy.updatePhase(fluidState_xy, phaseIdx);
for (int compIdx=0; compIdx<numComponents; ++compIdx){
auto phi_xy = PengRobinsonMixture::computeFugacityCoefficient(fluidState_xy, paramCache_xy, phaseIdx, compIdx);
fluidState_xy.setFugacityCoefficient(phaseIdx, compIdx, phi_xy);
auto f_xy = fluidState_xy.fugacity(phaseIdx, compIdx);
//std::cout << "comp" << compIdx <<" , f_xy " << f_xy << std::endl;
}
//RATIOS
for (int compIdx=0; compIdx<numComponents; ++compIdx){
if (isGas){
auto f_xy = fluidState_xy.fugacity(phaseIdx, compIdx);
auto f_zi = fluidState_zi.fugacity(oilPhaseIdx, compIdx);
auto fug_ratio = f_zi / f_xy;
R[compIdx] = fug_ratio / S_loc;
}
else{
auto fug_xy = fluidState_xy.fugacity(phaseIdx, compIdx);
auto fug_zi = fluidState_zi.fugacity(oilPhaseIdx, compIdx);
auto fug_ratio = fug_xy / fug_zi;
R[compIdx] = fug_ratio * S_loc;
}
}
Scalar R_norm = 0.0;
Scalar K_norm = 0.0;
for (int compIdx=0; compIdx<numComponents; ++compIdx){
K[compIdx] *= R[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) {
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);
isConverged = (R_norm < 1e-10);
bool ok = isTrivial || isConverged;
bool done = ok || i == maxIter;
if (done && !ok) {
isTrivial = true;
throw std::runtime_error(" Stability test did not converge");
//@warn "Stability test failed to converge in $maxiter iterations. Assuming stability." cond xy K_norm R_norm K
}
if (ok) {
stable = isTrivial || S_loc <= 1 + 1e-5;
return;
}
//todo: make sure that no mole fraction is smaller than 1e-8 ?
//todo: take care of water!
}
throw std::runtime_error(" Stability test did not converge");
}//end michelsen
template <class FlashFluidState, class ComponentVector> template <class FlashFluidState, class ComponentVector>
static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& globalComposition, int verbosity) static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluidState, const ComponentVector& globalComposition, int verbosity)
@ -848,7 +555,6 @@ template <class FlashFluidState, class ComponentVector>
fluidState_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal); fluidState_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
} }
ComponentVector R; ComponentVector R;
for (int compIdx=0; compIdx<numComponents; ++compIdx){ for (int compIdx=0; compIdx<numComponents; ++compIdx){
if (isGas){ if (isGas){
@ -1082,13 +788,7 @@ template <class FlashFluidState, class ComponentVector>
} }
} }
} }
// for (unsigned i = 0; i < num_equations; ++i) {
// for (unsigned j = 0; j < num_primary_variables; ++j) {
// std::cout << " " << jac[i][j] ;
// }
// std::cout << std::endl;
// }
std::cout << std::endl;
if (!converged) { if (!converged) {
throw std::runtime_error(" Newton composition update did not converge within maxIterations"); throw std::runtime_error(" Newton composition update did not converge within maxIterations");
} }
@ -1325,7 +1025,6 @@ template <class FlashFluidState, class ComponentVector>
(secondary_fluid_state, secondary_z, sec_jac, sec_res); (secondary_fluid_state, secondary_z, sec_jac, sec_res);
} }
// assembly the major matrix here // assembly the major matrix here
// primary variables are x, y and L // primary variables are x, y and L
constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1; constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1;
@ -1339,26 +1038,20 @@ template <class FlashFluidState, class ComponentVector>
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
primary_z[comp_idx] = Opm::getValue(z[comp_idx]); primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
} }
// TODO: x and y are not needed here
{
std::vector<PrimaryEval> x(numComponents), y(numComponents);
PrimaryEval l;
for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) { for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx); const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx);
x[comp_idx] = x_ii;//PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx); primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii);
primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x[comp_idx]);
const unsigned idx = comp_idx + numComponents; const unsigned idx = comp_idx + numComponents;
const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx); const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx);
y[comp_idx] = y_ii;//;PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx);
primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii); primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii);
primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii); primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii);
} }
PrimaryEval l;
l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1); l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1);
primary_fluid_state.setLvalue(l); primary_fluid_state.setLvalue(l);
}
primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx)); primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx));
primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx)); primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx));
primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0)); primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0));
// TODO: is PrimaryFlashFluidState::Scalar> PrimaryEval here? // TODO: is PrimaryFlashFluidState::Scalar> PrimaryEval here?
@ -1387,117 +1080,88 @@ template <class FlashFluidState, class ComponentVector>
(primary_fluid_state, primary_z, pri_jac, pri_res); (primary_fluid_state, primary_z, pri_jac, pri_res);
} }
//corresponds to julias J_p (we miss d/dt, and have d/dL instead of d/dV)
// for (unsigned i =0; i < num_equations; ++i) {
// for (unsigned j = 0; j < primary_num_pv; ++j) {
// std::cout << " " << pri_jac[i][j];
// }
// std::cout << std::endl;
// }
// std::cout << std::endl;
//corresponds to julias J_s
// for (unsigned i = 0; i < num_equations; ++i) {
// for (unsigned j = 0; j < secondary_num_pv; ++j) {
// std::cout << " " << sec_jac[i][j] ;
// }
// std::cout << std::endl;
// }
// std::cout << std::endl;
SecondaryNewtonMatrix xx; SecondaryNewtonMatrix xx;
pri_jac.solve(xx,sec_jac); pri_jac.solve(xx,sec_jac);
// for (unsigned i = 0; i < num_equations; ++i) {
// for (unsigned j = 0; j < secondary_num_pv; ++j) {
// std::cout << " " << xx[i][j] ;
// }
// std::cout << std::endl;
// }
using InputEval = typename FluidState::Scalar; using InputEval = typename FluidState::Scalar;
using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>; using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>;
//std::vector<InputEval> x(numComponents), y(numComponents);
ComponentVectorMoleFraction x(numComponents), y(numComponents); ComponentVectorMoleFraction x(numComponents), y(numComponents);
InputEval L_eval = L; InputEval L_eval = L;
// TODO: then beginning from that point
// use the chainrule (and using partial instead of total
// derivatives, DF / Dp = dF / dp + dF / ds * ds/dp.
// where p is the primary variables and s the secondary variables. We then obtain
// ds / dp = -inv(dF / ds)*(DF / Dp)
{ const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx);
const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx); const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx);
const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx); std::vector<double> K(numComponents);
std::vector<double> K(numComponents);
// const double L = fluid_state_scalar.L(); for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { K[compIdx] = fluid_state_scalar.K(compIdx);
K[compIdx] = fluid_state_scalar.K(compIdx); x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);//;z[compIdx] * 1. / (L + (1 - L) * K[compIdx]);
x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);//;z[compIdx] * 1. / (L + (1 - L) * K[compIdx]); y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);//;x[compIdx] * K[compIdx];
y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);//;x[compIdx] * K[compIdx];
}
// then we try to set the derivatives for x, y and K against P and x.
// p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary
// pressure joins
{
constexpr size_t num_deri = numComponents;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
std::vector<double> deri(num_deri, 0.);
// derivatives from P
// for (unsigned idx = 0; idx < num_deri; ++idx) {
// probably can use some DUNE operator for vectors or matrics
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[compIdx][0] * p_l.derivative(idx);
}
// }
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[compIdx][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
x[compIdx].setDerivative(idx, deri[idx]);
}
// handling y
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[compIdx + numComponents][0]* p_v.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[compIdx + numComponents][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
y[compIdx].setDerivative(idx, deri[idx]);
}
}
// handling derivatives of L
std::vector<double> deri(num_deri, 0.);
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[2*numComponents][0] * p_v.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[2*numComponents][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
L_eval.setDerivative(idx, deri[idx]);
}
}
} }
// x, y og L_eval
// then we try to set the derivatives for x, y and K against P and x.
// p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary
// pressure joins
constexpr size_t num_deri = numComponents;
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
std::vector<double> deri(num_deri, 0.);
// derivatives from P
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[compIdx][0] * p_l.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[compIdx][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
x[compIdx].setDerivative(idx, deri[idx]);
}
// handling y
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] = - xx[compIdx + numComponents][0]* p_v.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[compIdx + numComponents][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deri[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
y[compIdx].setDerivative(idx, deri[idx]);
}
// handling derivatives of L
std::vector<double> deriL(num_deri, 0.);
for (unsigned idx = 0; idx < num_deri; ++idx) {
deriL[idx] = - xx[2*numComponents][0] * p_v.derivative(idx);
}
for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
const double pz = -xx[2*numComponents][cIdx + 1];
const auto& zi = z[cIdx];
for (unsigned idx = 0; idx < num_deri; ++idx) {
deriL[idx] += pz * zi.derivative(idx);
}
}
for (unsigned idx = 0; idx < num_deri; ++idx) {
L_eval.setDerivative(idx, deriL[idx]);
}
// set up the mole fractions // set up the mole fractions
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
@ -1505,58 +1169,8 @@ template <class FlashFluidState, class ComponentVector>
fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]); fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
} }
fluid_state.setLvalue(L_eval); fluid_state.setLvalue(L_eval);
}
/* template <class Vector, class Matrix, class Eval, class ComponentVector>
static void evalJacobian(const ComponentVector& globalComposition,
const Vector& x,
const Vector& y,
const Eval& l,
Vector& b,
Matrix& m)
{
// TODO: all the things are going through the FluidState, which makes it difficult to get the AD correct.
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]);
} }
}//end updateDerivatives
// Compute fugacities
using ValueType = typename FluidState::Scalar;
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){
ValueType 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
ValueType 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 ComponentVector> template <class FluidState, class DefectVector, class ComponentVector>
static void evalDefect_(DefectVector& b, static void evalDefect_(DefectVector& b,
@ -1572,7 +1186,6 @@ template <class FlashFluidState, class ComponentVector>
fluidState.setMoleFraction(gasPhaseIdx, compIdx, x[compIdx + numMiscibleComponents]); fluidState.setMoleFraction(gasPhaseIdx, compIdx, x[compIdx + numMiscibleComponents]);
} }
// Compute fugacities // Compute fugacities
using ValueType = typename FluidState::Scalar; using ValueType = typename FluidState::Scalar;
using ParamCache = typename FluidSystem::template ParameterCache<typename FluidState::Scalar>; using ParamCache = typename FluidSystem::template ParameterCache<typename FluidState::Scalar>;
@ -1602,7 +1215,7 @@ template <class FlashFluidState, class ComponentVector>
// sum(x) - sum(y) = 0 // sum(x) - sum(y) = 0
b[numMiscibleComponents*numMisciblePhases] += -x[compIdx] + x[compIdx + numMiscibleComponents]; b[numMiscibleComponents*numMisciblePhases] += -x[compIdx] + x[compIdx + numMiscibleComponents];
} }
}//end valDefect }//end evalDefect
template <class FluidState, class DefectVector, class DefectMatrix, class ComponentVector> template <class FluidState, class DefectVector, class DefectMatrix, class ComponentVector>
static void evalJacobian_(DefectMatrix& A, static void evalJacobian_(DefectMatrix& A,
@ -1749,114 +1362,6 @@ template <class FlashFluidState, class ComponentVector>
// throw std::runtime_error("Successive substitution composition update did not converge within maxIterations"); // throw std::runtime_error("Successive substitution composition update did not converge within maxIterations");
} }
template <class FlashFluidState, class ComponentVector>
static void successiveSubstitutionCompositionNew_(ComponentVector& K, typename FlashFluidState::Scalar& L, FlashFluidState& fluidState, const ComponentVector& z,
const bool newton_afterwards, const int verbosity)
{
// std::cout << " Evaluation in successiveSubstitutionComposition_ is " << Dune::className(L) << std::endl;
// Determine max. iterations based on if it will be used as a standalone flash or as a pre-process to Newton (or other) method.
const int maxIterations = newton_afterwards ? 3 : 10;
// Store cout format before manipulation
std::ios_base::fmtflags f(std::cout.flags());
// Print initial guess
if (verbosity >= 1)
std::cout << "Initial guess: K = [" << K << "] and L = " << L << std::endl;
if (verbosity == 2 || verbosity == 4) {
// Print header
int fugWidth = (numComponents * 12)/2;
int convWidth = fugWidth + 7;
std::cout << std::setw(10) << "Iteration" << std::setw(fugWidth) << "fL/fV" << std::setw(convWidth) << "norm2(fL/fv-1)" << std::endl;
}
//
// Successive substitution loop
//
for (int i=0; i < maxIterations; ++i){
// Compute (normalized) liquid and vapor mole fractions
computeLiquidVapor_(fluidState, L, K, z);
// Calculate fugacity coefficient
using ParamCache = typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>;
ParamCache paramCache;
for (int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
paramCache.updatePhase(fluidState, phaseIdx);
for (int compIdx=0; compIdx<numComponents; ++compIdx){
auto 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, z, 0);
}
}
// throw std::runtime_error("Successive substitution composition update did not converge within maxIterations");
}
};//end ChiFlash };//end ChiFlash
} // namespace Opm } // namespace Opm

196
opm/material/fluidsystems/chifluid/co2brinefluidsystem.hh Normal file
View File

@ -0,0 +1,196 @@
#ifndef OPM_CO2BRINEFLUIDSYSTEM_HH
#define OPM_CO2BRINEFLUIDSYSTEM_HH
#include <opm/material/fluidsystems/BaseFluidSystem.hpp>
#include <opm/material/fluidsystems/chifluid/components.hh>
#include "ChiParameterCache.hpp"
#include "LBCviscosity.hpp"
namespace Opm {
/*!
* \ingroup FluidSystem
*
* \brief A two phase two component system, co2 brine
*/
template<class Scalar>
class Co2BrineFluidSystem
: public Opm::BaseFluidSystem<Scalar, Co2BrineFluidSystem<Scalar> > {
public:
// TODO: I do not think these should be constant in fluidsystem, will try to make it non-constant later
static constexpr int numPhases=2;
static constexpr int numComponents = 3;
// TODO: phase location should be more general
static constexpr int oilPhaseIdx = 0;
static constexpr int gasPhaseIdx = 1;
static constexpr int Comp0Idx = 0;
static constexpr int Comp1Idx = 1;
static constexpr int Comp2Idx = 2;
// TODO: needs to be more general
using Comp0 = Opm::JuliaCO2<Scalar>;
using Comp1 = Opm::ChiwomsBrine<Scalar>;
using Comp2 = Opm::JuliaC10<Scalar>;
template <class ValueType>
using ParameterCache = Opm::ChiParameterCache<ValueType, Co2BrineFluidSystem<Scalar>>;
using LBCviscosity = typename Opm::LBCviscosity<Scalar, Co2BrineFluidSystem<Scalar>>;
using PengRobinsonMixture = typename Opm::PengRobinsonMixture<Scalar, Co2BrineFluidSystem<Scalar>>;
/*!
* \brief The acentric factor of a component [].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar acentricFactor(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::acentricFactor();
case Comp1Idx: return Comp1::acentricFactor();
case Comp2Idx: return Comp2::acentricFactor();
default: throw std::runtime_error("Illegal component index for acentricFactor");
}
}
/*!
* \brief Critical temperature of a component [K].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar criticalTemperature(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::criticalTemperature();
case Comp1Idx: return Comp1::criticalTemperature();
case Comp2Idx: return Comp2::criticalTemperature();
default: throw std::runtime_error("Illegal component index for criticalTemperature");
}
}
/*!
* \brief Critical pressure of a component [Pa].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar criticalPressure(unsigned compIdx) {
switch (compIdx) {
case Comp0Idx: return Comp0::criticalPressure();
case Comp1Idx: return Comp1::criticalPressure();
case Comp2Idx: return Comp2::criticalPressure();
default: throw std::runtime_error("Illegal component index for criticalPressure");
}
}
/*!
* \brief Critical volume of a component [m3].
*
* \copydetails Doxygen::compIdxParam
*/
static Scalar criticalVolume(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::criticalVolume();
case Comp1Idx: return Comp1::criticalVolume();
case Comp2Idx: return Comp2::criticalVolume();
default: throw std::runtime_error("Illegal component index for criticalVolume");
}
}
//! \copydoc BaseFluidSystem::molarMass
static Scalar molarMass(unsigned compIdx)
{
switch (compIdx) {
case Comp0Idx: return Comp0::molarMass();
case Comp1Idx: return Comp1::molarMass();
case Comp2Idx: return Comp2::molarMass();
default: throw std::runtime_error("Illegal component index for molarMass");
}
}
/*!
* \brief Returns the interaction coefficient for two components.
*.
*/
static Scalar interactionCoefficient(unsigned /*comp1Idx*/, unsigned /*comp2Idx*/)
{
return 0.0;
}
//! \copydoc BaseFluidSystem::phaseName
static const char* phaseName(unsigned phaseIdx)
{
static const char* name[] = {"o", // oleic phase
"g"}; // gas phase
assert(0 <= phaseIdx && phaseIdx < 2);
return name[phaseIdx];
}
//! \copydoc BaseFluidSystem::componentName
static const char* componentName(unsigned compIdx)
{
static const char* name[] = {
Comp0::name(),
Comp1::name(),
Comp2::name(),
};
assert(0 <= compIdx && compIdx < 3);
return name[compIdx];
}
/*!
* \copydoc BaseFluidSystem::density
*/
template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
static LhsEval density(const FluidState& fluidState,
const ParameterCache<ParamCacheEval>& paramCache,
unsigned phaseIdx)
{
LhsEval dens;
if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) {
// paramCache.updatePhase(fluidState, phaseIdx);
dens = fluidState.averageMolarMass(phaseIdx) / paramCache.molarVolume(phaseIdx);
}
return dens;
}
//! \copydoc BaseFluidSystem::viscosity
template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
static LhsEval viscosity(const FluidState& fluidState,
const ParameterCache<ParamCacheEval>& paramCache,
unsigned phaseIdx)
{
// Use LBC method to calculate viscosity
// LhsEval mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx);
// LhsEval mu = LBCviscosity::LBC(fluidState, paramCache, phaseIdx);
LhsEval mu;
mu = LBCviscosity::LBCmod(fluidState, paramCache, phaseIdx);
// LhsEval mu = LBCviscosity::LBCJulia(fluidState, paramCache, phaseIdx);
return mu;
}
//! \copydoc BaseFluidSystem::fugacityCoefficient
template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
static LhsEval fugacityCoefficient(const FluidState& fluidState,
const ParameterCache<ParamCacheEval>& paramCache,
unsigned phaseIdx,
unsigned compIdx)
{
assert(0 <= phaseIdx && phaseIdx < numPhases);
assert(0 <= compIdx && compIdx < numComponents);
// TODO: here the derivatives for the phi are dropped. Should we keep the derivatives against the pressure
// and temperature?
LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
//Scalar phi = Opm::getValue(
// PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
return phi;
}
};
}
#endif //OPM_CO2BRINEFLUIDSYSTEM_HH

1
opm/material/fluidsystems/chifluid/components.hh
View File

@ -225,6 +225,7 @@ public:
// Critical volume [m3/kmol] // Critical volume [m3/kmol]
static Scalar criticalVolume() {return 9.412e-5; } static Scalar criticalVolume() {return 9.412e-5; }
// OLD :static Scalar criticalVolume() {return 9.4118e-2; }
}; };
template <class Scalar> template <class Scalar>