diff --git a/opm/material/common/PolynomialUtils.hpp b/opm/material/common/PolynomialUtils.hpp index 397fc7135..5a36b3645 100644 --- a/opm/material/common/PolynomialUtils.hpp +++ b/opm/material/common/PolynomialUtils.hpp @@ -26,6 +26,8 @@ #include #include +#include + namespace Opm { /*! * \ingroup Math @@ -75,6 +77,8 @@ int invertQuadraticPolynomial(SolContainer &sol, Scalar b, Scalar c) { + typedef MathToolbox Toolbox; + // check for a line if (a == 0.0) return invertLinearPolynomial(sol, b, c); @@ -84,7 +88,7 @@ int invertQuadraticPolynomial(SolContainer &sol, if (Delta < 0) return 0; // no real roots - Delta = std::sqrt(Delta); + Delta = Toolbox::sqrt(Delta); sol[0] = (- b + Delta)/(2*a); sol[1] = (- b - Delta)/(2*a); @@ -103,6 +107,8 @@ void invertCubicPolynomialPostProcess_(SolContainer &sol, Scalar c, Scalar d) { + typedef MathToolbox Toolbox; + // do one Newton iteration on the analytic solution if the // precision is increased for (int i = 0; i < numSol; ++i) { @@ -115,7 +121,7 @@ void invertCubicPolynomialPostProcess_(SolContainer &sol, x -= fOld/fPrime; Scalar fNew = d + x*(c + x*(b + x*a)); - if (std::abs(fNew) < std::abs(fOld)) + if (std::abs(Toolbox::value(fNew)) < std::abs(Toolbox::value(fOld))) sol[i] = x; } } @@ -145,6 +151,8 @@ int invertCubicPolynomial(SolContainer *sol, Scalar c, Scalar d) { + typedef MathToolbox Toolbox; + // reduces to a quadratic polynomial if (a == 0) return invertQuadraticPolynomial(sol, b, c, d); @@ -193,9 +201,9 @@ int invertCubicPolynomial(SolContainer *sol, Scalar wDisc = q*q/4 + p*p*p/27; if (wDisc >= 0) { // the positive discriminant case: // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27) - Scalar u = - q/2 + std::sqrt(wDisc); - if (u < 0) u = - std::pow(-u, 1.0/3); - else u = std::pow(u, 1.0/3); + Scalar u = - q/2 + Toolbox::sqrt(wDisc); + if (u < 0) u = - Toolbox::pow(-u, 1.0/3); + else u = Toolbox::pow(u, 1.0/3); // at this point, u != 0 since p^3 = 0 is necessary in order // for u = 0 to hold, so @@ -208,10 +216,10 @@ int invertCubicPolynomial(SolContainer *sol, } else { // the negative discriminant case: Scalar uCubedRe = - q/2; - Scalar uCubedIm = std::sqrt(-wDisc); + Scalar uCubedIm = Toolbox::sqrt(-wDisc); // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27) - Scalar uAbs = std::pow(std::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3); - Scalar phi = std::atan2(uCubedIm, uCubedRe)/3; + Scalar uAbs = Toolbox::pow(Toolbox::sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm), 1.0/3); + Scalar phi = Toolbox::atan2(uCubedIm, uCubedRe)/3; // calculate the length and the angle of the primitive root @@ -253,7 +261,7 @@ int invertCubicPolynomial(SolContainer *sol, // values for phi which differ by 2/3*pi. This allows to // calculate the three real roots of the polynomial: for (int i = 0; i < 3; ++i) { - sol[i] = std::cos(phi)*(uAbs - p/(3*uAbs)) - b/3; + sol[i] = Toolbox::cos(phi)*(uAbs - p/(3*uAbs)) - b/3; phi += 2*M_PI/3; } @@ -279,8 +287,8 @@ int invertCubicPolynomial(SolContainer *sol, // // i. e. single real root at t=curt(q) Scalar t; - if (-q > 0) t = std::pow(-q, 1./3); - else t = - std::pow(q, 1./3); + if (-q > 0) t = Toolbox::pow(-q, 1./3); + else t = - Toolbox::pow(q, 1./3); sol[0] = t - b/3; return 1; @@ -297,9 +305,9 @@ int invertCubicPolynomial(SolContainer *sol, } // two additional real roots at t = sqrt(-p) and t = -sqrt(-p) - sol[0] = -std::sqrt(-p) - b/3;; + sol[0] = -Toolbox::sqrt(-p) - b/3;; sol[1] = 0.0 - b/3; - sol[2] = std::sqrt(-p) - b/3; + sol[2] = Toolbox::sqrt(-p) - b/3; return 3; } diff --git a/opm/material/common/Spline.hpp b/opm/material/common/Spline.hpp index c90670ab1..003a9be66 100644 --- a/opm/material/common/Spline.hpp +++ b/opm/material/common/Spline.hpp @@ -812,6 +812,39 @@ public: return eval_(x, segmentIdx_(x)); } + /*! + * \brief Evaluate the spline for a given function evaluation. + * + * \param x The value on the abscissa where the spline ought to be + * evaluated + * \param extrapolate If this parameter is set to true, the spline + * will be extended beyond its range by + * straight lines, if false calling extrapolate + * for \f$ x \not [x_{min}, x_{max}]\f$ will + * cause a failed assertation. + */ + template + Evaluation eval(const Evaluation& x, bool extrapolate=false) const + { + assert(extrapolate || applies(x.value)); + + // handle extrapolation + if (extrapolate) { + if (x.value < xMin()) { + Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0); + Scalar y0 = y_(0); + return Evaluation::createConstant(y0 + m*(x.value - xMin())); + } + else if (x > xMax()) { + Scalar m = evalDerivative_(xMax(), /*segmentIdx=*/numSamples()-2); + Scalar y0 = y_(numSamples() - 1); + return Evaluation::createConstant(y0 + m*(x.value - xMax())); + } + } + + return eval_(x, segmentIdx_(x.value)); + } + /*! * \brief Evaluate the spline's derivative at a given position. * @@ -885,7 +918,11 @@ public: * Opm::MathError exception if there is more or less than * one solution. */ - Scalar intersect(Scalar a, Scalar b, Scalar c, Scalar d) const + template + Evaluation intersect(const Evaluation& a, + const Evaluation& b, + const Evaluation& c, + const Evaluation& d) const { return intersectInterval(xMin(), xMax(), a, b, c, d); } @@ -896,12 +933,16 @@ public: * Opm::MathError exception if there is more or less than * one solution. */ - Scalar intersectInterval(Scalar x0, Scalar x1, - Scalar a, Scalar b, Scalar c, Scalar d) const + template + Evaluation intersectInterval(Scalar x0, Scalar x1, + const Evaluation& a, + const Evaluation& b, + const Evaluation& c, + const Evaluation& d) const { assert(applies(x0) && applies(x1)); - Scalar tmpSol[3], sol = 0; + Evaluation tmpSol[3], sol = 0; int nSol = 0; int iFirst = segmentIdx_(x0); int iLast = segmentIdx_(x1); @@ -1485,6 +1526,29 @@ protected: + h11_(t) * slope_(i + 1)*delta; } + // evaluate the spline at a given the position and given the + // segment index + template + Evaluation eval_(const Evaluation& x, int i) const + { + // See http://en.wikipedia.org/wiki/Cubic_Hermite_spline + Scalar delta = h_(i + 1); + Scalar t = (x.value - x_(i))/delta; + + Evaluation result; + result.value = + h00_(t) * y_(i) + + h10_(t) * slope_(i)*delta + + h01_(t) * y_(i + 1) + + h11_(t) * slope_(i + 1)*delta; + + Scalar df_dg = evalDerivative_(x.value, i); + for (unsigned varIdx = 0; varIdx < result.derivatives.size(); ++ varIdx) + result.derivatives[varIdx] = df_dg*x.derivatives[varIdx]; + + return result; + } + // evaluate the derivative of a spline given the actual position // and the segment index Scalar evalDerivative_(Scalar x, int i) const @@ -1663,9 +1727,13 @@ protected: * \brief Find all the intersections of a segment of the spline * with a cubic polynomial within a specified interval. */ - int intersectSegment_(Scalar *sol, + template + int intersectSegment_(Evaluation *sol, int segIdx, - Scalar a, Scalar b, Scalar c, Scalar d, + const Evaluation& a, + const Evaluation& b, + const Evaluation& c, + const Evaluation& d, Scalar x0 = -1e100, Scalar x1 = 1e100) const { int n = Opm::invertCubicPolynomial(sol, @@ -1736,22 +1804,22 @@ protected: // returns the coefficient in front of the x^3 term. In Stoer this // is delta. Scalar a_(int i) const - { return evalDerivative3_(/*x=*/0, i)/6.0; } + { return evalDerivative3_(/*x=*/Scalar(0.0), i)/6.0; } // returns the coefficient in front of the x^2 term In Stoer this // is gamma. Scalar b_(int i) const - { return evalDerivative2_(/*x=*/0, i)/2.0; } + { return evalDerivative2_(/*x=*/Scalar(0.0), i)/2.0; } // returns the coefficient in front of the x^1 term. In Stoer this // is beta. Scalar c_(int i) const - { return evalDerivative_(/*x=*/0, i); } + { return evalDerivative_(/*x=*/Scalar(0.0), i); } // returns the coefficient in front of the x^0 term. In Stoer this // is alpha. Scalar d_(int i) const - { return eval_(/*x=*/0, i); } + { return eval_(/*x=*/Scalar(0.0), i); } Vector xPos_; Vector yPos_; diff --git a/opm/material/eos/PengRobinson.hpp b/opm/material/eos/PengRobinson.hpp index eeb914c07..9f9f87abb 100644 --- a/opm/material/eos/PengRobinson.hpp +++ b/opm/material/eos/PengRobinson.hpp @@ -452,7 +452,7 @@ protected: // invert resulting cubic polynomial analytically Scalar allV[4]; allV[0] = V; - int numSol = 1 + Opm::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); diff --git a/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp b/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp index 428a29b5d..83dcb09dd 100644 --- a/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp +++ b/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp @@ -244,7 +244,7 @@ public: const Spline& spline = params.pcnwHighSpline(); return spline.intersectInterval(/*x0=*/SwThHigh, /*x1=*/1.0, - /*a=*/0, /*b=*/0, /*c=*/0, /*d=*/pC); + /*a=*/0.0, /*b=*/0.0, /*c=*/0.0, /*d=*/pC); } return Sw; diff --git a/tests/test_spline.cpp b/tests/test_spline.cpp index 23ac2de49..27be79508 100644 --- a/tests/test_spline.cpp +++ b/tests/test_spline.cpp @@ -195,8 +195,8 @@ void testMonotonic(const Spline &sp, // test the intersection methods double d = (y[i] + y[i+1])/2; - double interX = sp.intersectInterval(x[i], x[i+1], - /*a=*/0, /*b=*/0, /*c=*/0, d); + double interX = sp.template intersectInterval(x[i], x[i+1], + /*a=*/0, /*b=*/0, /*c=*/0, d); double interY = sp.eval(interX); if (std::abs(interY - d) > 1e-5) OPM_THROW(std::runtime_error,