From 16f7fb8e9d560f104ce5ffad6d9c0d4b6e297ef4 Mon Sep 17 00:00:00 2001 From: Trine Mykkeltvedt Date: Wed, 25 May 2022 11:33:26 +0200 Subject: [PATCH] added Svenns cubic solver --- opm/material/common/MathToolbox.hpp | 40 +++++++++- opm/material/common/PolynomialUtils.hpp | 100 ++++++++++++++++++++++++ opm/material/densead/Math.hpp | 68 ++++++++++++++++ opm/material/eos/PengRobinson.hpp | 73 ++++++++--------- 4 files changed, 242 insertions(+), 39 deletions(-) diff --git a/opm/material/common/MathToolbox.hpp b/opm/material/common/MathToolbox.hpp index 94f4635cf..c34ba9005 100644 --- a/opm/material/common/MathToolbox.hpp +++ b/opm/material/common/MathToolbox.hpp @@ -100,7 +100,7 @@ public: * This basically boils down to creating an uninitialized object of sufficient size. * This is method only non-trivial for dynamically-sized Evaluation objects. */ - static Scalar createBlank(Scalar) + static Scalar createBlank(Scalar value OPM_UNUSED) { return Scalar(); } /*! @@ -136,7 +136,7 @@ public: * function. In general, this returns an evaluation object for which all derivatives * are zero. */ - static Scalar createConstant(Scalar, Scalar value) + static Scalar createConstant(Scalar x OPM_UNUSED, Scalar value) { return value; } /*! @@ -146,7 +146,7 @@ public: * regard to x. For scalars (which do not consider derivatives), this method does * nothing. */ - static Scalar createVariable(Scalar, unsigned) + static Scalar createVariable(Scalar value OPM_UNUSED, unsigned varIdx OPM_UNUSED) { throw std::logic_error("Plain floating point objects cannot represent variables"); } /*! @@ -157,7 +157,7 @@ public: * regard to x. For scalars (which do not consider derivatives), this method does * nothing. */ - static Scalar createVariable(Scalar, Scalar, unsigned) + static Scalar createVariable(Scalar x OPM_UNUSED, Scalar value OPM_UNUSED, unsigned varIdx OPM_UNUSED) { throw std::logic_error("Plain floating point objects cannot represent variables"); } /*! @@ -227,6 +227,14 @@ public: static Scalar asin(Scalar arg) { return std::asin(arg); } + //! The sine hyperbolicus of a value + static Scalar sinh(Scalar arg) + { return std::sinh(arg); } + + //! The arcus sine hyperbolicus of a value + static Scalar asinh(Scalar arg) + { return std::asinh(arg); } + //! The cosine of a value static Scalar cos(Scalar arg) { return std::cos(arg); } @@ -235,6 +243,14 @@ public: static Scalar acos(Scalar arg) { return std::acos(arg); } + //! The cosine hyperbolicus of a value + static Scalar cosh(Scalar arg) + { return std::cosh(arg); } + + //! The arcus cosine hyperbolicus of a value + static Scalar acosh(Scalar arg) + { return std::acosh(arg); } + //! The square root of a value static Scalar sqrt(Scalar arg) { return std::sqrt(arg); } @@ -357,6 +373,14 @@ template Evaluation asin(const Evaluation& value) { return MathToolbox::asin(value); } +template +Evaluation sinh(const Evaluation& value) +{ return MathToolbox::sinh(value); } + +template +Evaluation asinh(const Evaluation& value) +{ return MathToolbox::asinh(value); } + template Evaluation cos(const Evaluation& value) { return MathToolbox::cos(value); } @@ -365,6 +389,14 @@ template Evaluation acos(const Evaluation& value) { return MathToolbox::acos(value); } +template +Evaluation cosh(const Evaluation& value) +{ return MathToolbox::cosh(value); } + +template +Evaluation acosh(const Evaluation& value) +{ return MathToolbox::acosh(value); } + template Evaluation sqrt(const Evaluation& value) { return MathToolbox::sqrt(value); } diff --git a/opm/material/common/PolynomialUtils.hpp b/opm/material/common/PolynomialUtils.hpp index c5def1eec..451f0c4a8 100644 --- a/opm/material/common/PolynomialUtils.hpp +++ b/opm/material/common/PolynomialUtils.hpp @@ -309,6 +309,106 @@ unsigned invertCubicPolynomial(SolContainer* sol, return 3; } + +/*! + * \ingroup Math + * \brief Invert a cubic polynomial analytically + * + * The polynomial is defined as + * \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f] + * + * This method teturns the number of solutions which are in the real + * numbers. The "sol" argument contains the real roots of the cubic + * polynomial in order with the smallest root first. + * + * \param sol Container into which the solutions are written + * \param a The coefficient for the cubic term + * \param b The coefficient for the quadratic term + * \param c The coefficient for the linear term + * \param d The coefficient for the constant term + */ +template +unsigned cubicRoots(SolContainer* sol, + Scalar a, + Scalar b, + Scalar c, + Scalar d) +{ + // reduces to a quadratic polynomial + if (std::abs(scalarValue(a)) < 1e-30) + return invertQuadraticPolynomial(sol, b, c, d); + + // We need to reduce the cubic equation to its "depressed cubic" form (however strange that sounds) + // Depressed cubic form: t^3 + p*t + q, where x = t - b/3*a is the transform we use when we have + // roots for t. p and q are defined below. + // Formula for p and q: + Scalar p = (3.0 * a * c - b * b) / (3.0 * a * a); + Scalar q = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * d * a * a) / (27.0 * a * a * a); + + // Check if we have three or one real root by looking at the discriminant, and solve accordingly with + // correct formula + Scalar discr = 4.0 * p * p * p + 27.0 * q * q; + if (discr < 0.0) { + // Find three real roots of a depressed cubic, using the trigonometric method + // Help calculation + Scalar theta = (1.0 / 3.0) * acos( ((3.0 * q) / (2.0 * p)) * sqrt(-3.0 / p) ); + + // Calculate the three roots + sol[0] = 2.0 * sqrt(-p / 3.0) * cos( theta ); + sol[1] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((2.0 * M_PI) / 3.0) ); + sol[2] = 2.0 * sqrt(-p / 3.0) * cos( theta - ((4.0 * M_PI) / 3.0) ); + + // Sort in ascending order + std::sort(sol, sol + 3); + + // Return confirmation of three roots + // std::cout << "Z (discr < 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl; + return 3; + } + else if (discr > 0.0) { + // Find one real root of a depressed cubic using hyperbolic method. Different solutions depending on + // sign of p + Scalar t; + if (p < 0) { + // Help calculation + Scalar theta = (1.0 / 3.0) * acosh( ((-3.0 * abs(q)) / (2.0 * p)) * sqrt(-3.0 / p) ); + + // Root + t = ( (-2.0 * abs(q)) / q ) * sqrt(-p / 3.0) * cosh(theta); + } + else if (p > 0) { + // Help calculation + Scalar theta = (1.0 / 3.0) * asinh( ((3.0 * q) / (2.0 * p)) * sqrt(3.0 / p) ); + // Root + t = -2.0 * sqrt(p / 3.0) * sinh(theta); + + } + else { + std::runtime_error(" p = 0 in cubic root solver!"); + } + + // Transform t to output solution + sol[0] = t - b / (3.0 * a); + // std::cout << "Z (discr > 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl; + return 1; + + } + else { + // The discriminant, 4*p^3 + 27*q^2 = 0, thus we have simple (real) roots + // If p = 0 then also q = 0, and t = 0 is a triple root + if (p == 0) { + sol[0] = sol[1] = sol[2] = 0.0 - b / (3.0 * a); + } + // If p != 0, the we have a simple root and a double root + else { + sol[0] = (3.0 * q / p) - b / (3.0 * a); + sol[1] = sol[2] = (-3.0 * q) / (2.0 * p) - b / (3.0 * a); + std::sort(sol, sol + 3); + } + // std::cout << "Z (disc = 0) = " << sol[0] << " " << sol[1] << " " << sol[2] << std::endl; + return 3; + } } +} // end Opm #endif diff --git a/opm/material/densead/Math.hpp b/opm/material/densead/Math.hpp index 05c151eee..6899f3001 100644 --- a/opm/material/densead/Math.hpp +++ b/opm/material/densead/Math.hpp @@ -225,6 +225,40 @@ Evaluation asin(const Evaluation +Evaluation sinh(const Evaluation& x) +{ + typedef MathToolbox ValueTypeToolbox; + + Evaluation result(x); + + result.setValue(ValueTypeToolbox::sinh(x.value())); + + // derivatives use the chain rule + const ValueType& df_dx = ValueTypeToolbox::cosh(x.value()); + for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx) + result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); + + return result; +} + +template +Evaluation asinh(const Evaluation& x) +{ + typedef MathToolbox ValueTypeToolbox; + + Evaluation result(x); + + result.setValue(ValueTypeToolbox::asinh(x.value())); + + // derivatives use the chain rule + const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(x.value()*x.value() + 1); + for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx) + result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); + + return result; +} + template Evaluation cos(const Evaluation& x) { @@ -259,6 +293,40 @@ Evaluation acos(const Evaluation +Evaluation cosh(const Evaluation& x) +{ + typedef MathToolbox ValueTypeToolbox; + + Evaluation result(x); + + result.setValue(ValueTypeToolbox::cosh(x.value())); + + // derivatives use the chain rule + const ValueType& df_dx = ValueTypeToolbox::sinh(x.value()); + for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx) + result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); + + return result; +} + +template +Evaluation acosh(const Evaluation& x) +{ + typedef MathToolbox ValueTypeToolbox; + + Evaluation result(x); + + result.setValue(ValueTypeToolbox::acosh(x.value())); + + // derivatives use the chain rule + const ValueType& df_dx = 1.0/ValueTypeToolbox::sqrt(x.value()*x.value() - 1); + for (int curVarIdx = 0; curVarIdx < result.size(); ++curVarIdx) + result.setDerivative(curVarIdx, df_dx*x.derivative(curVarIdx)); + + return result; +} + template Evaluation sqrt(const Evaluation& x) { diff --git a/opm/material/eos/PengRobinson.hpp b/opm/material/eos/PengRobinson.hpp index 46438525d..191e9a5d3 100644 --- a/opm/material/eos/PengRobinson.hpp +++ b/opm/material/eos/PengRobinson.hpp @@ -55,6 +55,8 @@ namespace Opm { template class PengRobinson { + //! The ideal gas constant [Pa * m^3/mol/K] + static const Scalar R; PengRobinson() { } @@ -127,7 +129,7 @@ public: const Evaluation& delta = f/df_dp; pVap = pVap - delta; - if (std::abs(scalarValue(delta/pVap)) < 1e-10) + if (std::abs(Opm::scalarValue(delta/pVap)) < 1e-10) break; } @@ -160,13 +162,13 @@ public: const Evaluation& a = params.a(phaseIdx); // "attractive factor" const Evaluation& b = params.b(phaseIdx); // "co-volume" - if (!std::isfinite(scalarValue(a)) - || std::abs(scalarValue(a)) < 1e-30) + if (!std::isfinite(Opm::scalarValue(a)) + || std::abs(Opm::scalarValue(a)) < 1e-30) return std::numeric_limits::quiet_NaN(); - if (!std::isfinite(scalarValue(b)) || b <= 0) + if (!std::isfinite(Opm::scalarValue(b)) || b <= 0) return std::numeric_limits::quiet_NaN(); - const Evaluation& RT= Constants::R*T; + const Evaluation& RT= R*T; const Evaluation& Astar = a*p/(RT*RT); const Evaluation& Bstar = b*p/RT; @@ -184,27 +186,26 @@ public: Valgrind::CheckDefined(a2); Valgrind::CheckDefined(a3); Valgrind::CheckDefined(a4); - int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4); + // std::cout << "Cubic params : " << a1 << " " << a2 << " " << a3 << " " << a4 << std::endl; + // int numSol = invertCubicPolynomial(Z, a1, a2, a3, a4); + int numSol = cubicRoots(Z, a1, a2, a3, a4); + // std::cout << "Z = " << Z[0] << " " << Z[1] << " " << Z[2] << std::endl; if (numSol == 3) { // the EOS has three intersections with the pressure, // i.e. the molar volume of gas is the largest one and the // molar volume of liquid is the smallest one -#warning HACK, should investigate why if (isGasPhase) - // Vm = Z[2]*RT/p; - Vm = max(1e-7, Z[2]*RT/p); + Vm = Opm::max(1e-7, Z[2]*RT/p); else - // Vm = Z[2]*RT/p; - Vm = max(1e-7, Z[0]*RT/p); + Vm = Opm::max(1e-7, Z[0]*RT/p); } else if (numSol == 1) { // the EOS only has one intersection with the pressure, // for the other phase, we take the extremum of the EOS // with the largest distance from the intersection. - Evaluation VmCubic = Z[0]*RT/p; + Evaluation VmCubic = Opm::max(1e-7, Z[0]*RT/p); Vm = VmCubic; -#warning, should investigate why here // find the extrema (if they are present) // Evaluation Vmin, Vmax, pmin, pmax; // if (findExtrema_(Vmin, Vmax, @@ -229,7 +230,7 @@ public: } Valgrind::CheckDefined(Vm); - assert(std::isfinite(scalarValue(Vm))); + assert(Opm::isfinite(Vm)); assert(Vm > 0); return Vm; } @@ -251,7 +252,7 @@ public: const Evaluation& p = params.pressure(); const Evaluation& Vm = params.molarVolume(); - const Evaluation& RT = Constants::R*T; + const Evaluation& RT = R*T; const Evaluation& Z = p*Vm/RT; const Evaluation& Bstar = p*params.b() / RT; @@ -260,8 +261,8 @@ public: (Vm + params.b()*(1 - std::sqrt(2))); const Evaluation& expo = - params.a()/(RT * 2 * params.b() * std::sqrt(2)); const Evaluation& fugCoeff = - exp(Z - 1) / (Z - Bstar) * - pow(tmp, expo); + Opm::exp(Z - 1) / (Z - Bstar) * + Opm::pow(tmp, expo); return fugCoeff; } @@ -299,9 +300,9 @@ protected: //Evaluation Vcrit = criticalMolarVolume_.eval(params.a(phaseIdx), params.b(phaseIdx)); if (isGasPhase) - Vm = max(Vm, Vcrit); + Vm = Opm::max(Vm, Vcrit); else - Vm = min(Vm, Vcrit); + Vm = Opm::min(Vm, Vcrit); } template @@ -351,14 +352,14 @@ protected: const Scalar eps = - 1e-11; bool hasExtrema OPM_OPTIM_UNUSED = findExtrema_(minVm, maxVm, minP, maxP, a, b, T + eps); assert(hasExtrema); - assert(std::isfinite(scalarValue(maxVm))); + assert(std::isfinite(Opm::scalarValue(maxVm))); Evaluation fStar = maxVm - minVm; // derivative of the difference between the maximum's // molar volume and the minimum's molar volume regarding // temperature Evaluation fPrime = (fStar - f)/eps; - if (std::abs(scalarValue(fPrime)) < 1e-40) { + if (std::abs(Opm::scalarValue(fPrime)) < 1e-40) { Tcrit = T; pcrit = (minP + maxP)/2; Vcrit = (maxVm + minVm)/2; @@ -367,7 +368,7 @@ protected: // update value for the current iteration Evaluation delta = f/fPrime; - assert(std::isfinite(scalarValue(delta))); + assert(std::isfinite(Opm::scalarValue(delta))); if (delta > 0) delta = -10; @@ -415,7 +416,7 @@ protected: Scalar u = 2; Scalar w = -1; - const Evaluation& RT = Constants::R*T; + const Evaluation& RT = R*T; // calculate coefficients of the 4th order polynominal in // monomial basis @@ -425,11 +426,11 @@ protected: 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; - assert(std::isfinite(scalarValue(a1))); - assert(std::isfinite(scalarValue(a2))); - assert(std::isfinite(scalarValue(a3))); - assert(std::isfinite(scalarValue(a4))); - assert(std::isfinite(scalarValue(a5))); + assert(std::isfinite(Opm::scalarValue(a1))); + assert(std::isfinite(Opm::scalarValue(a2))); + assert(std::isfinite(Opm::scalarValue(a3))); + assert(std::isfinite(Opm::scalarValue(a4))); + assert(std::isfinite(Opm::scalarValue(a5))); // Newton method to find first root @@ -438,11 +439,11 @@ protected: // above the covolume Evaluation V = b*1.1; Evaluation delta = 1.0; - for (unsigned i = 0; std::abs(scalarValue(delta)) > 1e-12; ++i) { + for (unsigned i = 0; std::abs(Opm::scalarValue(delta)) > 1e-12; ++i) { 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)); - if (std::abs(scalarValue(fPrime)) < 1e-20) { + if (std::abs(Opm::scalarValue(fPrime)) < 1e-20) { // give up if the derivative is zero return false; } @@ -456,7 +457,7 @@ protected: return false; } } - assert(std::isfinite(scalarValue(V))); + assert(std::isfinite(Opm::scalarValue(V))); // polynomial division Evaluation b1 = a1; @@ -467,7 +468,7 @@ protected: // invert resulting cubic polynomial analytically Evaluation allV[4]; allV[0] = V; - int numSol = 1 + invertCubicPolynomial(allV + 1, b1, b2, b3, b4); + int numSol = 1 + Opm::invertCubicPolynomial(allV + 1, b1, b2, b3, b4); // sort all roots of the derivative std::sort(allV + 0, allV + numSol); @@ -507,9 +508,9 @@ protected: const Evaluation& tau = 1 - Tr; const Evaluation& omega = Component::acentricFactor(); - 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 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr; - const Evaluation& f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr; + const Evaluation& f0 = (tau*(-5.97616 + Opm::sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*Opm::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& f2 = (tau*(-0.64771 + Opm::sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*Opm::pow(tau, 5))/Tr; return Component::criticalPressure()*std::exp(f0 + omega * (f1 + omega*f2)); } @@ -539,6 +540,8 @@ protected: */ }; +template +const Scalar PengRobinson::R = Opm::Constants::R; /* template