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

added new cubic solver from Svenns old branch, also made a phaseStabilityTestMichelsen_ with the same wrapping as julia code - this does not work with newton right now. BUT the fix of the roots makes the stabilitytest give similar K as Olavs Julia code. Will never be completely equal due to minimization through gibbs (not implemented in opm) which will different choize of roots

Browse Source
This commit is contained in:
Trine Mykkeltvedt 2022-06-01 12:23:29 +02:00
parent 16f7fb8e9d
commit eed51f4d55
3 changed files with 41 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
opm/material/common/PolynomialUtils.hpp
View File

@ -354,18 +354,21 @@ unsigned cubicRoots(SolContainer* sol,
Scalar theta = (1.0 / 3.0) * acos( ((3.0 * q) / (2.0 * p)) * sqrt(-3.0 / p) ); Scalar theta = (1.0 / 3.0) * acos( ((3.0 * q) / (2.0 * p)) * sqrt(-3.0 / p) );
// Calculate the three roots // Calculate the three roots
sol[0] = 2.0 * sqrt(-p / 3.0) * cos( theta ); sol[0] = 2.0 * sqrt(-p / 3.0) * cos( theta ) - b / (3.0 * a);
sol[1] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((2.0 * M_PI) / 3.0) ); sol[1] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((2.0 * M_PI) / 3.0) ) - b / (3.0 * a);
sol[2] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((4.0 * M_PI) / 3.0) ); sol[2] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((4.0 * M_PI) / 3.0) ) - b / (3.0 * a);
//std::cout << "Z (discr < 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl;
// Sort in ascending order // Sort in ascending order
std::sort(sol, sol + 3); std::sort(sol, sol + 3);
// Return confirmation of three roots // Return confirmation of three roots
// std::cout << "Z (discr < 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl;
return 3; return 3;
} }
else if (discr > 0.0) { else if (discr > 0.0) {
// Find one real root of a depressed cubic using hyperbolic method. Different solutions depending on // Find one real root of a depressed cubic using hyperbolic method. Different solutions depending on
// sign of p // sign of p
Scalar t; Scalar t;

63
opm/material/eos/PengRobinson.hpp
View File

@ -129,7 +129,7 @@ public:
const Evaluation& delta = f/df_dp; const Evaluation& delta = f/df_dp;
pVap = pVap - delta; pVap = pVap - delta;
if (std::abs(Opm::scalarValue(delta/pVap)) < 1e-10) if (std::abs(scalarValue(delta/pVap)) < 1e-10)
break; break;
} }
@ -162,13 +162,13 @@ public:
const Evaluation& a = params.a(phaseIdx); // "attractive factor" const Evaluation& a = params.a(phaseIdx); // "attractive factor"
const Evaluation& b = params.b(phaseIdx); // "co-volume" const Evaluation& b = params.b(phaseIdx); // "co-volume"
if (!std::isfinite(Opm::scalarValue(a)) if (!std::isfinite(scalarValue(a))
|| std::abs(Opm::scalarValue(a)) < 1e-30) || std::abs(scalarValue(a)) < 1e-30)
return std::numeric_limits<Scalar>::quiet_NaN(); return std::numeric_limits<Scalar>::quiet_NaN();
if (!std::isfinite(Opm::scalarValue(b)) || b <= 0) if (!std::isfinite(scalarValue(b)) || b <= 0)
return std::numeric_limits<Scalar>::quiet_NaN(); return std::numeric_limits<Scalar>::quiet_NaN();
const Evaluation& RT= R*T; const Evaluation& RT= Constants<Scalar>::R*T;
const Evaluation& Astar = a*p/(RT*RT); const Evaluation& Astar = a*p/(RT*RT);
const Evaluation& Bstar = b*p/RT; const Evaluation& Bstar = b*p/RT;
@ -195,15 +195,15 @@ public:
// i.e. the molar volume of gas is the largest one and the // i.e. the molar volume of gas is the largest one and the
// molar volume of liquid is the smallest one // molar volume of liquid is the smallest one
if (isGasPhase) if (isGasPhase)
Vm = Opm::max(1e-7, Z[2]*RT/p); Vm = max(1e-7, Z[2]*RT/p);
else else
Vm = Opm::max(1e-7, Z[0]*RT/p); Vm = max(1e-7, Z[0]*RT/p);
} }
else if (numSol == 1) { else if (numSol == 1) {
// the EOS only has one intersection with the pressure, // the EOS only has one intersection with the pressure,
// for the other phase, we take the extremum of the EOS // for the other phase, we take the extremum of the EOS
// with the largest distance from the intersection. // with the largest distance from the intersection.
Evaluation VmCubic = Opm::max(1e-7, Z[0]*RT/p); Evaluation VmCubic = max(1e-7, Z[0]*RT/p);
Vm = VmCubic; Vm = VmCubic;
// find the extrema (if they are present) // find the extrema (if they are present)
@ -230,7 +230,7 @@ public:
} }
Valgrind::CheckDefined(Vm); Valgrind::CheckDefined(Vm);
assert(Opm::isfinite(Vm)); assert(std::isfinite(scalarValue(Vm)));
assert(Vm > 0); assert(Vm > 0);
return Vm; return Vm;
} }
@ -252,7 +252,7 @@ public:
const Evaluation& p = params.pressure(); const Evaluation& p = params.pressure();
const Evaluation& Vm = params.molarVolume(); const Evaluation& Vm = params.molarVolume();
const Evaluation& RT = R*T; const Evaluation& RT = Constants<Scalar>::R*T;
const Evaluation& Z = p*Vm/RT; const Evaluation& Z = p*Vm/RT;
const Evaluation& Bstar = p*params.b() / RT; const Evaluation& Bstar = p*params.b() / RT;
@ -261,8 +261,8 @@ public:
(Vm + params.b()*(1 - std::sqrt(2))); (Vm + params.b()*(1 - std::sqrt(2)));
const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2)); const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2));
const Evaluation& fugCoeff = const Evaluation& fugCoeff =
Opm::exp(Z - 1) / (Z - Bstar) * exp(Z - 1) / (Z - Bstar) *
Opm::pow(tmp, expo); pow(tmp, expo);
return fugCoeff; return fugCoeff;
} }
@ -300,9 +300,9 @@ protected:
//Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx)); //Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx));
if (isGasPhase) if (isGasPhase)
Vm = Opm::max(Vm, Vcrit); Vm = max(Vm, Vcrit);
else else
Vm = Opm::min(Vm, Vcrit); Vm = min(Vm, Vcrit);
} }
template <class Evaluation> template <class Evaluation>
@ -352,14 +352,14 @@ protected:
const Scalar eps = - 1e-11; const Scalar eps = - 1e-11;
bool hasExtrema OPM_OPTIM_UNUSED = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps); bool hasExtrema OPM_OPTIM_UNUSED = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps);
assert(hasExtrema); assert(hasExtrema);
assert(std::isfinite(Opm::scalarValue(maxVm))); assert(std::isfinite(scalarValue(maxVm)));
Evaluation fStar = maxVm - minVm; Evaluation fStar = maxVm - minVm;
// derivative of the difference between the maximum's // derivative of the difference between the maximum's
// molar volume and the minimum's molar volume regarding // molar volume and the minimum's molar volume regarding
// temperature // temperature
Evaluation fPrime = (fStar - f)/eps; Evaluation fPrime = (fStar - f)/eps;
if (std::abs(Opm::scalarValue(fPrime)) < 1e-40) { if (std::abs(scalarValue(fPrime)) < 1e-40) {
Tcrit = T; Tcrit = T;
pcrit = (minP + maxP)/2; pcrit = (minP + maxP)/2;
Vcrit = (maxVm + minVm)/2; Vcrit = (maxVm + minVm)/2;
@ -368,7 +368,7 @@ protected:
// update value for the current iteration // update value for the current iteration
Evaluation delta = f/fPrime; Evaluation delta = f/fPrime;
assert(std::isfinite(Opm::scalarValue(delta))); assert(std::isfinite(scalarValue(delta)));
if (delta > 0) if (delta > 0)
delta = -10; delta = -10;
@ -416,8 +416,7 @@ protected:
Scalar u = 2; Scalar u = 2;
Scalar w = -1; Scalar w = -1;
const Evaluation& RT = R*T; const Evaluation& RT = Constants<Scalar>::R*T;
// calculate coefficients of the 4th order polynominal in // calculate coefficients of the 4th order polynominal in
// monomial basis // monomial basis
const Evaluation& a1 = RT; const Evaluation& a1 = RT;
@ -426,11 +425,11 @@ protected:
const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b; const Evaluation& a4 = 2*RT*u*w*b*b*b + 2*u*a*b*b - 2*a*b*b;
const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b; const Evaluation& a5 = RT*w*w*b*b*b*b - u*a*b*b*b;
assert(std::isfinite(Opm::scalarValue(a1))); assert(std::isfinite(scalarValue(a1)));
assert(std::isfinite(Opm::scalarValue(a2))); assert(std::isfinite(scalarValue(a2)));
assert(std::isfinite(Opm::scalarValue(a3))); assert(std::isfinite(scalarValue(a3)));
assert(std::isfinite(Opm::scalarValue(a4))); assert(std::isfinite(scalarValue(a4)));
assert(std::isfinite(Opm::scalarValue(a5))); assert(std::isfinite(scalarValue(a5)));
// Newton method to find first root // Newton method to find first root
@ -439,11 +438,11 @@ protected:
// above the covolume // above the covolume
Evaluation V = b*1.1; Evaluation V = b*1.1;
Evaluation delta = 1.0; Evaluation delta = 1.0;
for (unsigned i = 0; std::abs(Opm::scalarValue(delta)) > 1e-12; ++i) { for (unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) {
const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1))); const Evaluation& f = a5 + V*(a4 + V*(a3 + V*(a2 + V*a1)));
const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1)); const Evaluation& fPrime = a4 + V*(2*a3 + V*(3*a2 + V*4*a1));
if (std::abs(Opm::scalarValue(fPrime)) < 1e-20) { if (std::abs(scalarValue(fPrime)) < 1e-20) {
// give up if the derivative is zero // give up if the derivative is zero
return false; return false;
} }
@ -457,7 +456,7 @@ protected:
return false; return false;
} }
} }
assert(std::isfinite(Opm::scalarValue(V))); assert(std::isfinite(scalarValue(V)));
// polynomial division // polynomial division
Evaluation b1 = a1; Evaluation b1 = a1;
@ -468,7 +467,7 @@ protected:
// invert resulting cubic polynomial analytically // invert resulting cubic polynomial analytically
Evaluation allV[4]; Evaluation allV[4];
allV[0] = V; allV[0] = V;
int numSol = 1 + Opm::invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4); int numSol = 1 + invertCubicPolynomial<Evaluation>(allV + 1, b1, b2, b3, b4);
// sort all roots of the derivative // sort all roots of the derivative
std::sort(allV + 0, allV + numSol); std::sort(allV + 0, allV + numSol);
@ -508,9 +507,9 @@ protected:
const Evaluation& tau = 1 - Tr; const Evaluation& tau = 1 - Tr;
const Evaluation& omega = Component::acentricFactor(); const Evaluation& omega = Component::acentricFactor();
const Evaluation& f0 = (tau*(-5.97616 + Opm::sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*Opm::pow(tau, 5))/Tr; const Evaluation& f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
const Evaluation& f1 = (tau*(-5.03365 + Opm::sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*Opm::pow(tau, 5))/Tr; const Evaluation& f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
const Evaluation& f2 = (tau*(-0.64771 + Opm::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*Opm::pow(tau, 5))/Tr; const Evaluation& f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;
return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2)); return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2));
} }
@ -540,10 +539,10 @@ protected:
*/ */
}; };
/*
template <class Scalar> template <class Scalar>
const Scalar PengRobinson<Scalar>::R = Opm::Constants<Scalar>::R; const Scalar PengRobinson<Scalar>::R = Opm::Constants<Scalar>::R;
/*
template <class Scalar> template <class Scalar>
UniformTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalTemperature_; UniformTabulated2DFunction<Scalar> PengRobinson<Scalar>::criticalTemperature_;

4
tests/test_chiflash.cpp
View File

@ -117,8 +117,8 @@ void testChiFlash()
// TODO: only, p, z need the derivatives. // TODO: only, p, z need the derivatives.
const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional const double flash_tolerance = 1.e-12; // just to test the setup in co2-compositional
const int flash_verbosity = 1; const int flash_verbosity = 1;
const std::string flash_twophase_method = "newton"; // "ssi" //const std::string flash_twophase_method = "ssi"; // "ssi"
// const std::string flash_twophase_method = "ssi"; const std::string flash_twophase_method = "newton";
// const std::string flash_twophase_method = "ssi+newton"; // const std::string flash_twophase_method = "ssi+newton";
// TODO: should we set these? // TODO: should we set these?