fixed the derivatives in PengRobinsonMixture with a rewrite of a term to get pri_jac correct, also a getValue to get sec_jac equal, did rewriting and printing in ugly way jus t compare with Julia-code values
This commit is contained in:
Trine Mykkeltvedt
parent
03685e5bf2
commit
03f5a39985
@ -955,6 +955,9 @@ protected:
|
||||
auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
|
||||
fluid_state.fugacity(gasPhaseIdx, compIdx));
|
||||
res[compIdx + numComponents] = Opm::getValue(local_res);
|
||||
//std::cout << "fugacity oil = " << local_res.derivative(i) << " gas = " << fluid_state.fugacity(gasPhaseIdx, compIdx) << " comp " << compIdx << std::endl; trine
|
||||
|
||||
|
||||
for (unsigned i = 0; i < num_primary; ++i) {
|
||||
jac[compIdx + numComponents][i] = local_res.derivative(i);
|
||||
}
|
||||
@ -1084,11 +1087,66 @@ protected:
|
||||
assembleNewton_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
|
||||
(primary_fluid_state, primary_z, pri_jac, pri_res);
|
||||
|
||||
// not totally sure the following matrix operations are correct
|
||||
pri_jac.invert();
|
||||
sec_jac.template leftmultiply(pri_jac);
|
||||
//corresponds to julias J_s
|
||||
std::cout << "sec_jac:" << std::endl;
|
||||
std::cout << "[" << sec_jac[0][0] << " " << sec_jac[0][1] << " " << sec_jac[0][2] << " " << sec_jac[0][3] << "] " << std::endl;
|
||||
std::cout << "[" << sec_jac[1][0] << " " << sec_jac[1][1] << " " << sec_jac[1][2] << " " << sec_jac[1][3] << "] " << std::endl;
|
||||
std::cout << "[" << sec_jac[2][0] << " " << sec_jac[2][1] << " " << sec_jac[2][2] << " " << sec_jac[2][3] << "] " << std::endl;
|
||||
std::cout << "[" << sec_jac[3][0] << " " << sec_jac[3][1] << " " << sec_jac[3][2] << " " << sec_jac[3][3] << "] " << std::endl;
|
||||
std::cout << "[" << sec_jac[4][0] << " " << sec_jac[4][1] << " " << sec_jac[4][2] << " " << sec_jac[4][3] << "] " << std::endl;
|
||||
std::cout << "[" << sec_jac[5][0] << " " << sec_jac[5][1] << " " << sec_jac[5][2] << " " << sec_jac[5][3] << "] " << std::endl;
|
||||
std::cout << "[" << sec_jac[6][0] << " " << sec_jac[6][1] << " " << sec_jac[6][2] << " " << sec_jac[6][3] << "] " << std::endl;
|
||||
|
||||
//corresponds to julias J_p (we miss d/dt, and have d/dL instead of d/dV)
|
||||
std::cout << "pri_jac:" << std::endl;
|
||||
std::cout << "[" << pri_jac[0][0] << " " << pri_jac[0][1] << " " << pri_jac[0][2] << " " << pri_jac[0][3] << " " << pri_jac[0][4]<< " " << pri_jac[0][5] << " " << pri_jac[0][6]<< "] " << std::endl;
|
||||
std::cout << "[" << pri_jac[1][0] << " " << pri_jac[1][1] << " " << pri_jac[1][2] << " " << pri_jac[1][3] << " " << pri_jac[1][4]<< " " << pri_jac[1][5] << " " << pri_jac[1][6]<< "] " << std::endl;
|
||||
std::cout << "[" << pri_jac[2][0] << " " << pri_jac[2][1] << " " << pri_jac[2][2] << " " << pri_jac[2][3] << " " << pri_jac[2][4]<< " " << pri_jac[2][5] << " " << pri_jac[2][6]<< "] " << std::endl;
|
||||
std::cout << "[" << pri_jac[3][0] << " " << pri_jac[3][1] << " " << pri_jac[3][2] << " " << pri_jac[3][3] << " " << pri_jac[3][4]<< " " << pri_jac[3][5] << " " << pri_jac[3][6]<< "] " << std::endl;
|
||||
std::cout << "[" << pri_jac[4][0] << " " << pri_jac[4][1] << " " << pri_jac[4][2] << " " << pri_jac[4][3] << " " << pri_jac[4][4]<< " " << pri_jac[4][5] << " " << pri_jac[4][6]<< "] " << std::endl;
|
||||
std::cout << "[" << pri_jac[5][0] << " " << pri_jac[5][1] << " " << pri_jac[5][2] << " " << pri_jac[5][3] << " " << pri_jac[5][4]<< " " << pri_jac[5][5] << " " << pri_jac[5][6]<< "] " << std::endl;
|
||||
std::cout << "[" << pri_jac[6][0] << " " << pri_jac[6][1] << " " << pri_jac[6][2] << " " << pri_jac[6][3] << " " << pri_jac[6][4]<< " " << pri_jac[6][5] << " " << pri_jac[6][6]<< "] " << std::endl;
|
||||
|
||||
|
||||
SecondaryNewtonMatrix xx;
|
||||
pri_jac.solve(xx,sec_jac);
|
||||
std::cout << " corresponding to julia-code value and updated J_s " << std::endl;
|
||||
std::cout << "x1 = [" << x[0] << " " << xx[0][0] << " " << xx[0][1] << " " << xx[0][2] << " " << xx[0][3] << "] " << std::endl;
|
||||
std::cout << "x2 = [" << x[1] << " " << xx[1][0] << " " <<xx[1][1] << " " << xx[1][2] << " " << xx[1][3] << "] " << std::endl;
|
||||
std::cout << "x3 = [" << x[2] << " " << xx[2][0] << " " << xx[2][1] << " " << xx[2][2] << " " << xx[2][3] << "] " << std::endl;
|
||||
std::cout << "y1 = [" << y[0] << " " << xx[3][0] << " " << xx[3][1] << " " << xx[3][2] << " " << xx[3][3] << "] " << std::endl;
|
||||
std::cout << "y2 = [" << y[1] << " " << xx[4][0] << " " << xx[4][1] << " " << xx[4][2] << " " << xx[4][3] << "] " << std::endl;
|
||||
std::cout << "y3 = [" << y[2] << " " << xx[5][0] << " " << xx[5][1] << " " << xx[5][2] << " " << xx[5][3] << "] " << std::endl;
|
||||
std::cout << "L = [" << L << " " << xx[6][0] << " " << xx[6][1] << " " << xx[6][2] << " " << xx[6][3] << "] " << std::endl;
|
||||
|
||||
// rewrite like Olav xx --> xxx (do this properly to clean up =)
|
||||
// z3 = 1 -z2 -z1;
|
||||
// dx1/dp dx1/dt (dx1/dz1-dx1/dz3) (dx1/dz2 - dx1/dz3)
|
||||
using TertiaryMatrix = Dune::FieldMatrix<Scalar, num_equations, secondary_num_pv-1>;
|
||||
TertiaryMatrix xxx;
|
||||
xxx[0][0] = xx[0][0];
|
||||
xxx[1][0] = xx[1][0];
|
||||
xxx[2][0] = xx[2][0];
|
||||
xxx[3][0] = xx[0][0];
|
||||
xxx[4][0] = xx[1][0];
|
||||
xxx[5][0] = xx[2][0];
|
||||
xxx[6][0] = xx[3][0];
|
||||
for (unsigned i = 0; i < primary_num_pv; ++i) { // 7 rekker
|
||||
xxx[i][1] = Opm::getValue(xx[i][1])-Opm::getValue(xx[i][3]);
|
||||
xxx[i][2] = Opm::getValue(xx[i][2])-Opm::getValue(xx[i][3]);
|
||||
}
|
||||
std::cout << " corresponding to julia-code value and derivatives listed in test_setup NB CHANGE SIGN, and we dont have d/dT " << std::endl;
|
||||
std::cout << "x1 = [" << x[0] << " " << xxx[0][0] << " " << xxx[0][1] << " " << xxx[0][2] << "] " << std::endl;
|
||||
std::cout << "x2 = [" << x[1] << " " << xxx[1][0] << " " <<xxx[1][1] << " " << xxx[1][2] << "] " << std::endl;
|
||||
std::cout << "x3 = [" << x[2] << " " << xxx[2][0] << " " << xxx[2][1] << " " << xxx[2][2] << "] " << std::endl;
|
||||
std::cout << "y1 = [" << y[0] << " " << xxx[3][0] << " " << xxx[3][1] << " " << xxx[3][2] << "] " << std::endl;
|
||||
std::cout << "y2 = [" << y[1] << " " << xxx[4][0] << " " << xxx[4][1] << " " << xxx[4][2] << "] " << std::endl;
|
||||
std::cout << "y3 = [" << y[2] << " " << xxx[5][0] << " " << xxx[5][1] << " " << xxx[5][2] << "] " << std::endl;
|
||||
std::cout << "L = [" << L << " " << xxx[6][0] << " " << xxx[6][1] << " " << xxx[6][2] << " " << xx[6][3] << "] " << std::endl;
|
||||
|
||||
// TODO: then beginning from that point
|
||||
}
|
||||
|
||||
}//end updateDerivative
|
||||
|
||||
/* template <class Vector, class Matrix, class Eval, class ComponentVector>
|
||||
static void evalJacobian(const ComponentVector& globalComposition,
|
||||
|
||||
@ -94,21 +94,21 @@ public:
|
||||
unsigned phaseIdx,
|
||||
unsigned compIdx)
|
||||
{
|
||||
#warning also a HACK, should investigate why
|
||||
auto fs = fs2;
|
||||
double sumx = 0.0;
|
||||
for (int i = 0; i < FluidState::numComponents; ++i)
|
||||
sumx += Opm::scalarValue(fs.moleFraction(phaseIdx, i));
|
||||
if (sumx < 0.95) {
|
||||
double alpha = 0.95/sumx;
|
||||
std::cerr << "normalize: " << sumx
|
||||
<< " alpha: " << alpha << "\n";
|
||||
for (int i = 0; i < FluidState::numComponents; ++i)
|
||||
fs.setMoleFraction(phaseIdx, i, alpha*fs.moleFraction(phaseIdx, i));
|
||||
}
|
||||
// #warning also a HACK, should investigate why
|
||||
auto fs = fs2;
|
||||
// double sumx = 0.0;
|
||||
// for (int i = 0; i < FluidState::numComponents; ++i)
|
||||
// sumx += Opm::scalarValue(fs.moleFraction(phaseIdx, i));
|
||||
// if (sumx < 0.95) {
|
||||
// double alpha = 0.95/sumx;
|
||||
// std::cerr << "normalize: " << sumx
|
||||
// << " alpha: " << alpha << "\n";
|
||||
// for (int i = 0; i < FluidState::numComponents; ++i)
|
||||
// fs.setMoleFraction(phaseIdx, i, alpha*fs.moleFraction(phaseIdx, i));
|
||||
// }
|
||||
|
||||
auto params = params2;
|
||||
params.updatePhase(fs, phaseIdx);
|
||||
auto params = params2;
|
||||
// params.updatePhase(fs, phaseIdx);
|
||||
// note that we normalize the component mole fractions, so
|
||||
// that their sum is 100%. This increases numerical stability
|
||||
// considerably if the fluid state is not physical.
|
||||
@ -126,29 +126,69 @@ public:
|
||||
LhsEval Astar = params.a(phaseIdx)*p/(RT*RT);
|
||||
LhsEval Bstar = params.b(phaseIdx)*p/(RT);
|
||||
|
||||
//MAYBE REMOVE THIS COMPLETELY ??
|
||||
// calculate delta_i (see: Reid, p. 145)
|
||||
LhsEval sumMoleFractions = 0.0;
|
||||
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx)
|
||||
sumMoleFractions += fs.moleFraction(phaseIdx, compJIdx);
|
||||
LhsEval deltai = 2*sqrt(params.aPure(phaseIdx, compIdx))/params.a(phaseIdx);
|
||||
LhsEval tmp = 0;
|
||||
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
|
||||
tmp +=
|
||||
fs.moleFraction(phaseIdx, compJIdx)
|
||||
/ sumMoleFractions
|
||||
* sqrt(params.aPure(phaseIdx, compJIdx))
|
||||
* (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
|
||||
};
|
||||
deltai *= tmp;
|
||||
//std::cout << "params.aPure [" << params.aPure(phaseIdx, compIdx) << "] " << std::endl;
|
||||
LhsEval tmp = 0;
|
||||
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
|
||||
tmp +=
|
||||
fs.moleFraction(phaseIdx, compJIdx)
|
||||
/ sumMoleFractions
|
||||
* sqrt(params.aPure(phaseIdx, compJIdx));
|
||||
//* (1.0 - StaticParameters::interactionCoefficient(compIdx, compJIdx));
|
||||
};
|
||||
deltai *= tmp;
|
||||
|
||||
// Calculate A_s LIKE IN JULIA CODE
|
||||
|
||||
LhsEval A_s = 0.0;
|
||||
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx)
|
||||
{
|
||||
//param
|
||||
A_s += params.aCache(phaseIdx, compIdx, compJIdx) * fs.moleFraction(phaseIdx, compJIdx) * p / (RT * RT);
|
||||
}
|
||||
//LhsEval A_s = 0.5*Astar*deltai;
|
||||
|
||||
|
||||
LhsEval julia_alpha;
|
||||
LhsEval julia_betta;
|
||||
LhsEval julia_gamma;
|
||||
LhsEval ln_phi;
|
||||
LhsEval fugCoeff;
|
||||
|
||||
Scalar m1;
|
||||
Scalar m2;
|
||||
|
||||
m1 = 0.5*(u + std::sqrt(u*u - 4*w));
|
||||
m2 = 0.5*(u - std::sqrt(u*u - 4*w));
|
||||
|
||||
|
||||
LhsEval base =
|
||||
(2*Z + Bstar*(u + std::sqrt(u*u - 4*w))) /
|
||||
(2*Z + Bstar*(u - std::sqrt(u*u - 4*w)));
|
||||
(Z + Bstar*m1) /
|
||||
(Z + Bstar*m2);
|
||||
LhsEval expo = Astar/(Bstar*std::sqrt(u*u - 4*w))*(bi_b - deltai);
|
||||
|
||||
LhsEval fugCoeff =
|
||||
exp(bi_b*(Z - 1))/max(1e-9, Z - Bstar) *
|
||||
pow(base, expo);
|
||||
//julia:
|
||||
//α = -log(Z - B) + (B_i[i]/B)*(Z - 1)
|
||||
//β = log((Z + m2*B)/(Z + m1*B))*A/(Δm*B)
|
||||
//γ = (2/A)*A_s - B_i[i]/B
|
||||
//return α + β*γ
|
||||
|
||||
julia_alpha = -log(Z-Bstar)+bi_b*(Z - 1);
|
||||
julia_betta = log((Z+m2*Bstar)/(Z+m1*Bstar))*Astar/((m1-m2)*Bstar);
|
||||
//julia_gamma = -bi_b+deltai;//expo;//(2/Astar)*A_s -bi_b/Bstar;
|
||||
julia_gamma = (2 / Astar )* A_s - bi_b;
|
||||
|
||||
ln_phi = julia_alpha + (julia_betta*julia_gamma);
|
||||
|
||||
|
||||
|
||||
fugCoeff = exp(ln_phi);
|
||||
//exp(bi_b*(Z - 1))/max(1e-9, Z - Bstar) * pow(base, expo);
|
||||
|
||||
////////
|
||||
// limit the fugacity coefficient to a reasonable range:
|
||||
|
||||
@ -66,6 +66,15 @@ class PengRobinsonParamsMixture
|
||||
typedef MathToolbox<Scalar> Toolbox;
|
||||
|
||||
public:
|
||||
|
||||
/*!
|
||||
* \brief TODO
|
||||
*/
|
||||
Scalar getaCache(unsigned compIIdx, unsigned compJIdx )
|
||||
{
|
||||
return aCache_[compIIdx][compJIdx];
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Update Peng-Robinson parameters for the pure components.
|
||||
*/
|
||||
@ -99,12 +108,12 @@ public:
|
||||
|
||||
Scalar f_omega;
|
||||
|
||||
if (useSpe5Relations) {
|
||||
// if (useSpe5Relations) {
|
||||
if (omega < 0.49) f_omega = 0.37464 + omega*(1.54226 + omega*(-0.26992));
|
||||
else f_omega = 0.379642 + omega*(1.48503 + omega*(-0.164423 + omega*0.016666));
|
||||
}
|
||||
else
|
||||
f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
|
||||
// }
|
||||
// else
|
||||
// f_omega = 0.37464 + omega*(1.54226 - omega*0.26992);
|
||||
|
||||
Valgrind::CheckDefined(f_omega);
|
||||
|
||||
@ -135,6 +144,7 @@ public:
|
||||
template <class FluidState>
|
||||
void updateMix(const FluidState& fs)
|
||||
{
|
||||
using FlashEval = typename FluidState::Scalar;
|
||||
Scalar sumx = 0.0;
|
||||
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
|
||||
sumx += fs.moleFraction(phaseIdx, compIdx);
|
||||
@ -144,16 +154,16 @@ public:
|
||||
//
|
||||
// See: R. Reid, et al.: The Properties of Gases and Liquids,
|
||||
// 4th edition, McGraw-Hill, 1987, p. 82
|
||||
Scalar newA = 0;
|
||||
Scalar newB = 0;
|
||||
FlashEval newA = 0;
|
||||
FlashEval newB = 0;
|
||||
for (unsigned compIIdx = 0; compIIdx < numComponents; ++compIIdx) {
|
||||
const Scalar moleFracI = fs.moleFraction(phaseIdx, compIIdx);
|
||||
Scalar xi = max(0.0, min(1.0, moleFracI));
|
||||
const FlashEval moleFracI = fs.moleFraction(phaseIdx, compIIdx);
|
||||
FlashEval xi = max(0.0, min(1.0, moleFracI));
|
||||
Valgrind::CheckDefined(xi);
|
||||
|
||||
for (unsigned compJIdx = 0; compJIdx < numComponents; ++compJIdx) {
|
||||
const Scalar moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
|
||||
Scalar xj = max(0.0, min(1.0, moleFracJ));
|
||||
const FlashEval moleFracJ = fs.moleFraction(phaseIdx, compJIdx );
|
||||
FlashEval xj = max(0.0, min(1.0, moleFracJ));
|
||||
Valgrind::CheckDefined(xj);
|
||||
|
||||
// mixing rule from Reid, page 82
|
||||
|
||||
@ -177,6 +177,27 @@ public:
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*!
|
||||
* \brief TODO
|
||||
*
|
||||
* \param phaseIdx The fluid phase of interest
|
||||
* \param compIdx The component phase of interest
|
||||
* \param compJIdx Additional component index
|
||||
*/
|
||||
Scalar aCache(unsigned phaseIdx, unsigned compIdx, unsigned compJIdx)
|
||||
{
|
||||
switch (phaseIdx)
|
||||
{
|
||||
case oilPhaseIdx: return oilPhaseParams_.getaCache(compIdx,compJIdx);
|
||||
case gasPhaseIdx: return gasPhaseParams_.getaCache(compIdx,compJIdx);
|
||||
default:
|
||||
throw std::logic_error("The aCache() parameter is only defined for "
|
||||
"oil phase");
|
||||
};
|
||||
}
|
||||
|
||||
/*!
|
||||
* \brief Returns the molar volume of a phase [m^3/mol]
|
||||
*
|
||||
|
||||
@ -160,8 +160,9 @@ namespace Opm {
|
||||
|
||||
// TODO: here the derivatives for the phi are dropped. Should we keep the derivatives against the pressure
|
||||
// and temperature?
|
||||
Scalar phi = Opm::getValue(
|
||||
PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
|
||||
LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
|
||||
//Scalar phi = Opm::getValue(
|
||||
// PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
|
||||
return phi;
|
||||
}
|
||||
|
||||
|
||||
@ -272,8 +272,9 @@ public:
|
||||
assert(0 <= phaseIdx && phaseIdx < numPhases);
|
||||
assert(0 <= compIdx && compIdx < numComponents);
|
||||
|
||||
Scalar phi = Opm::getValue(
|
||||
PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
|
||||
LhsEval phi = PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
|
||||
//Scalar phi = Opm::getValue(
|
||||
// PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
|
||||
return phi;
|
||||
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user