diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C index 7848ca0d8b..f3f1658a7d 100644 --- a/applications/test/ODE/Test-ODE.C +++ b/applications/test/ODE/Test-ODE.C @@ -109,8 +109,11 @@ int main(int argc, char *argv[]) // Create the ODE system testODE ode; + dictionary dict; + dict.add("solver", args[1]); + // Create the selected ODE system solver - autoPtr odeSolver = ODESolver::New(args[1], ode); + autoPtr odeSolver = ODESolver::New(ode, dict); // Initialise the ODE system fields scalar xStart = 1.0; @@ -124,26 +127,26 @@ int main(int argc, char *argv[]) scalarField dyStart(ode.nEqns()); ode.derivatives(xStart, yStart, dyStart); - Info<< setw(10) << "eps" << setw(12) << "hEst"; - Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl; + Info<< setw(10) << "relTol" << setw(12) << "dxEst"; + Info<< setw(13) << "dxDid" << setw(14) << "dxNext" << endl; Info<< setprecision(6); for (label i=0; i<15; i++) { - scalar eps = ::Foam::exp(-scalar(i + 1)); + scalar relTol = ::Foam::exp(-scalar(i + 1)); scalar x = xStart; scalarField y(yStart); - scalarField dydx(dyStart); - scalar hEst = 0.6; - scalar hDid, hNext; + scalar dxEst = 0.6; + scalar dxNext = dxEst; - odeSolver->solve(ode, x, y, dydx, eps, hEst, hDid, hNext); + odeSolver->relTol() = relTol; + odeSolver->solve(ode, x, y, dxNext); - Info<< scientific << setw(13) << eps; - Info<< fixed << setw(11) << hEst; - Info<< setw(13) << hDid << setw(13) << hNext + Info<< scientific << setw(13) << relTol; + Info<< fixed << setw(11) << dxEst; + Info<< setw(13) << x - xStart << setw(13) << dxNext << setw(13) << y[0] << setw(13) << y[1] << setw(13) << y[2] << setw(13) << y[3] << endl; @@ -159,12 +162,13 @@ int main(int argc, char *argv[]) yEnd[2] = ::Foam::jn(2, xEnd); yEnd[3] = ::Foam::jn(3, xEnd); - scalar hEst = 0.5; + scalar dxEst = 0.5; - odeSolver->solve(ode, x, xEnd, y, 1e-4, hEst); + odeSolver->relTol() = 1e-4; + odeSolver->solve(ode, x, xEnd, y, dxEst); Info<< nl << "Analytical: y(2.0) = " << yEnd << endl; - Info << "Numerical: y(2.0) = " << y << ", hEst = " << hEst << endl; + Info << "Numerical: y(2.0) = " << y << ", dxEst = " << dxEst << endl; Info<< "\nEnd\n" << endl; diff --git a/src/ODE/Make/files b/src/ODE/Make/files index cb569e78f5..726921457c 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -2,7 +2,8 @@ ODESolvers/ODESolver/ODESolver.C ODESolvers/ODESolver/ODESolverNew.C ODESolvers/adaptiveSolver/adaptiveSolver.C -ODESolvers/RKCK5/RKCK5.C +ODESolvers/RKF45/RKF45.C +ODESolvers/RKCK45/RKCK45.C ODESolvers/KRR4/KRR4.C ODESolvers/SIBS/SIBS.C ODESolvers/SIBS/SIMPR.C diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C index 6f0b938782..db1190d1d6 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.C +++ b/src/ODE/ODESolvers/KRR4/KRR4.C @@ -31,7 +31,7 @@ License namespace Foam { defineTypeNameAndDebug(KRR4, 0); - addToRunTimeSelectionTable(ODESolver, KRR4, ODE); + addToRunTimeSelectionTable(ODESolver, KRR4, dictionary); const scalar KRR4::gamma = 1.0/2.0, @@ -50,15 +50,16 @@ const scalar // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::KRR4::KRR4(const ODESystem& ode) +Foam::KRR4::KRR4(const ODESystem& ode, const dictionary& dict) : - ODESolver(ode), - adaptiveSolver(ode), + ODESolver(ode, dict), + adaptiveSolver(ode, dict), g1_(n_), g2_(n_), g3_(n_), g4_(n_), err_(n_), + dydx_(n_), dfdx_(n_), dfdy_(n_, n_), a_(n_, n_), @@ -75,8 +76,7 @@ Foam::scalar Foam::KRR4::solve const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const { ode.jacobian(x0, y0, dfdx_, dfdy_); @@ -142,7 +142,7 @@ Foam::scalar Foam::KRR4::solve err_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i]; } - return normalizeError(eps, y0, y, err_); + return normalizeError(y0, y, err_); } @@ -151,14 +151,10 @@ void Foam::KRR4::solve const ODESystem& odes, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const { - adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext); + adaptiveSolver::solve(odes, x, y, dxTry); } diff --git a/src/ODE/ODESolvers/KRR4/KRR4.H b/src/ODE/ODESolvers/KRR4/KRR4.H index 710651fde9..5aff37095f 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.H +++ b/src/ODE/ODESolvers/KRR4/KRR4.H @@ -75,6 +75,7 @@ class KRR4 mutable scalarField g3_; mutable scalarField g4_; mutable scalarField err_; + mutable scalarField dydx_; mutable scalarField dfdx_; mutable scalarSquareMatrix dfdy_; mutable scalarSquareMatrix a_; @@ -99,7 +100,7 @@ public: // Constructors //- Construct from ODE - KRR4(const ODESystem& ode); + KRR4(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -112,8 +113,7 @@ public: const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const; //- Solve the ODE system and the update the state @@ -122,11 +122,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 60a599e8d2..00db4bf3c4 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C @@ -30,76 +30,103 @@ License namespace Foam { defineTypeNameAndDebug(ODESolver, 0); - defineRunTimeSelectionTable(ODESolver, ODE); + defineRunTimeSelectionTable(ODESolver, dictionary); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::ODESolver::ODESolver(const ODESystem& ode) +Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict) : n_(ode.nEqns()), - dydx_(n_) + absTol_(n_, dict.lookupOrDefault("absTol", SMALL)), + relTol_(n_, dict.lookupOrDefault("relTol", 1e-6)), + maxSteps_(10000) +{} + + +Foam::ODESolver::ODESolver +( + const ODESystem& ode, + const scalarField& absTol, + const scalarField& relTol +) +: + n_(ode.nEqns()), + absTol_(absTol), + relTol_(relTol), + maxSteps_(10000) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +Foam::scalar Foam::ODESolver::normalizeError +( + const scalarField& y0, + const scalarField& y, + const scalarField& err +) const +{ + // Calculate the maximum error + scalar maxErr = 0.0; + forAll(err, i) + { + scalar tol = absTol_[i] + relTol_[i]*max(mag(y0[i]), mag(y[i])); + maxErr = max(maxErr, mag(err[i])/tol); + } + + return maxErr; +} + + void Foam::ODESolver::solve ( const ODESystem& ode, const scalar xStart, const scalar xEnd, scalarField& y, - const scalar eps, - scalar& hEst + scalar& dxEst ) const { - const label MAXSTP = 10000; - scalar x = xStart; - scalar h = hEst; - scalar hNext = 0; - scalar hPrev = 0; + bool truncated = false; - for (label nStep=0; nStep 0.0) + // Check if this is a truncated step and set dxEst to integrate to xEnd + if ((x + dxEst - xEnd)*(x + dxEst - xStart) > 0) { - h = xEnd - x; - hPrev = hNext; + truncated = true; + dxEst = xEnd - x; } - hNext = 0; - scalar hDid; - solve(ode, x, y, dydx_, eps, h, hDid, hNext); + // Integrate as far as possible up to dxEst + solve(ode, x, y, dxEst); - if ((x - xEnd)*(xEnd - xStart) >= 0.0) + // Check if reached xEnd + if ((x - xEnd)*(xEnd - xStart) >= 0) { - if (hPrev != 0) + if (nStep > 0 && truncated) { - hEst = hPrev; - } - else - { - hEst = hNext; + dxEst = dxEst0; } return; } - - h = hNext; } FatalErrorIn ( "ODESolver::solve" "(const ODESystem& ode, const scalar xStart, const scalar xEnd," - "scalarField& yStart, const scalar eps, scalar& hEst) const" - ) << "Too many integration steps" + "scalarField& y, scalar& dxEst) const" + ) << "Integration steps greater than maximum " << maxSteps_ << exit(FatalError); } + // ************************************************************************* // diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H index a293088685..bd767a9063 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H @@ -53,13 +53,30 @@ class ODESolver protected: - // Private data + // Protected data + //- Size of the ODESystem label n_; - mutable scalarField dydx_; + + //- Absolute convergence tolerance per step + scalarField absTol_; + + //- Relative convergence tolerance per step + scalarField relTol_; + + //- The maximum number of sub-steps allowed for the integration step + label maxSteps_; - // Private Member Functions + // Protected Member Functions + + //- Return the nomalized scalar error + scalar normalizeError + ( + const scalarField& y0, + const scalarField& y, + const scalarField& err + ) const; //- Disallow default bitwise copy construct ODESolver(const ODESolver&); @@ -80,16 +97,24 @@ public: ( autoPtr, ODESolver, - ODE, - (const ODESystem& ode), - (ode) + dictionary, + (const ODESystem& ode, const dictionary& dict), + (ode, dict) ); // Constructors - //- Construct for given ODE - ODESolver(const ODESystem& ode); + //- Construct for given ODESystem + ODESolver(const ODESystem& ode, const dictionary& dict); + + //- Construct for given ODESystem specifying tolerances + ODESolver + ( + const ODESystem& ode, + const scalarField& absTol, + const scalarField& relTol + ); // Selectors @@ -97,8 +122,8 @@ public: //- Select null constructed static autoPtr New ( - const word& ODESolverTypeName, - const ODESystem& ode + const ODESystem& ode, + const dictionary& dict ); @@ -109,27 +134,37 @@ public: // Member Functions + scalarField& absTol() + { + return absTol_; + } + + scalarField& relTol() + { + return relTol_; + } + + //- Solve the ODE system as far as possible upto dxTry + // adjusting the step as necessary to provide a solution within + // the specified tolerance. + // Update the state and return an estimate for the next step in dxTry virtual void solve ( const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const = 0; - + //- Solve the ODE system from xStart to xEnd, update the state + // and return an estimate for the next step in dxTry virtual void solve ( const ODESystem& ode, const scalar xStart, const scalar xEnd, scalarField& y, - const scalar eps, - scalar& hEst + scalar& dxEst ) const; }; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C index aa8ed538cc..ef0d6e3454 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C @@ -29,29 +29,30 @@ License Foam::autoPtr Foam::ODESolver::New ( - const Foam::word& ODESolverTypeName, - const Foam::ODESystem& ode + const ODESystem& odes, + const dictionary& dict ) { + word ODESolverTypeName(dict.lookup("solver")); Info<< "Selecting ODE solver " << ODESolverTypeName << endl; - ODEConstructorTable::iterator cstrIter = - ODEConstructorTablePtr_->find(ODESolverTypeName); + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(ODESolverTypeName); - if (cstrIter == ODEConstructorTablePtr_->end()) + if (cstrIter == dictionaryConstructorTablePtr_->end()) { FatalErrorIn ( "ODESolver::New" - "(const word& ODESolverTypeName, const ODESystem& ode)" + "(const dictionary& dict, const ODESystem& odes)" ) << "Unknown ODESolver type " << ODESolverTypeName << nl << nl << "Valid ODESolvers are : " << endl - << ODEConstructorTablePtr_->sortedToc() + << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError); } - return autoPtr(cstrIter()(ode)); + return autoPtr(cstrIter()(odes, dict)); } diff --git a/src/ODE/ODESolvers/RKCK5/RKCK5.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C similarity index 51% rename from src/ODE/ODESolvers/RKCK5/RKCK5.C rename to src/ODE/ODESolvers/RKCK45/RKCK45.C index 5050c0e7ca..2de3cdf99e 100644 --- a/src/ODE/ODESolvers/RKCK5/RKCK5.C +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -23,41 +23,58 @@ License \*---------------------------------------------------------------------------*/ -#include "RKCK5.H" +#include "RKCK45.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(RKCK5, 0); - addToRunTimeSelectionTable(ODESolver, RKCK5, ODE); + defineTypeNameAndDebug(RKCK45, 0); + addToRunTimeSelectionTable(ODESolver, RKCK45, dictionary); const scalar - RKCK5::a2 = 0.2, RKCK5::a3 = 0.3, RKCK5::a4 = 0.6, RKCK5::a5 = 1.0, - RKCK5::a6 = 0.875, - RKCK5::b21 = 0.2, RKCK5::b31 = 3.0/40.0, RKCK5::b32 = 9.0/40.0, - RKCK5::b41 = 0.3, RKCK5::b42 = -0.9, RKCK5::b43 = 1.2, - RKCK5::b51 = -11.0/54.0, RKCK5::b52 = 2.5, RKCK5::b53 = -70.0/27.0, - RKCK5::b54 = 35.0/27.0, - RKCK5::b61 = 1631.0/55296.0, RKCK5::b62 = 175.0/512.0, - RKCK5::b63 = 575.0/13824.0, RKCK5::b64 = 44275.0/110592.0, - RKCK5::b65 = 253.0/4096.0, - RKCK5::c1 = 37.0/378.0, RKCK5::c3 = 250.0/621.0, - RKCK5::c4 = 125.0/594.0, RKCK5::c6 = 512.0/1771.0, - RKCK5::dc1 = RKCK5::c1 - 2825.0/27648.0, - RKCK5::dc3 = RKCK5::c3 - 18575.0/48384.0, - RKCK5::dc4 = RKCK5::c4 - 13525.0/55296.0, RKCK5::dc5 = -277.00/14336.0, - RKCK5::dc6 = RKCK5::c6 - 0.25; + RKCK45::c2 = 1.0/5.0, + RKCK45::c3 = 3.0/10.0, + RKCK45::c4 = 3.0/5.0, + RKCK45::c5 = 1.0, + RKCK45::c6 = 7.0/8.0, + + RKCK45::a21 = 1.0/5.0, + RKCK45::a31 = 3.0/40.0, + RKCK45::a32 = 9.0/40.0, + RKCK45::a41 = 3.0/10.0, + RKCK45::a42 = -9.0/10.0, + RKCK45::a43 = 6.0/5.0, + RKCK45::a51 = -11.0/54.0, + RKCK45::a52 = 5.0/2.0, + RKCK45::a53 = -70.0/27.0, + RKCK45::a54 = 35.0/27.0, + RKCK45::a61 = 1631.0/55296.0, + RKCK45::a62 = 175.0/512.0, + RKCK45::a63 = 575.0/13824.0, + RKCK45::a64 = 44275.0/110592.0, + RKCK45::a65 = 253.0/4096.0, + + RKCK45::b1 = 37.0/378.0, + RKCK45::b3 = 250.0/621.0, + RKCK45::b4 = 125.0/594.0, + RKCK45::b6 = 512.0/1771.0, + + RKCK45::e1 = RKCK45::b1 - 2825.0/27648.0, + RKCK45::e3 = RKCK45::b3 - 18575.0/48384.0, + RKCK45::e4 = RKCK45::b4 - 13525.0/55296.0, + RKCK45::e5 = -277.00/14336.0, + RKCK45::e6 = RKCK45::b6 - 0.25; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::RKCK5::RKCK5(const ODESystem& ode) +Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict) : - ODESolver(ode), - adaptiveSolver(ode), + ODESolver(ode, dict), + adaptiveSolver(ode, dict), yTemp_(n_), k2_(n_), k3_(n_), @@ -70,85 +87,80 @@ Foam::RKCK5::RKCK5(const ODESystem& ode) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::RKCK5::solve +Foam::scalar Foam::RKCK45::solve ( const ODESystem& odes, const scalar x0, const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const { forAll(yTemp_, i) { - yTemp_[i] = y0[i] + b21*dx*dydx0[i]; + yTemp_[i] = y0[i] + a21*dx*dydx0[i]; } - odes.derivatives(x0 + a2*dx, yTemp_, k2_); + odes.derivatives(x0 + c2*dx, yTemp_, k2_); forAll(yTemp_, i) { - yTemp_[i] = y0[i] + dx*(b31*dydx0[i] + b32*k2_[i]); + yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]); } - odes.derivatives(x0 + a3*dx, yTemp_, k3_); + odes.derivatives(x0 + c3*dx, yTemp_, k3_); forAll(yTemp_, i) { - yTemp_[i] = y0[i] + dx*(b41*dydx0[i] + b42*k2_[i] + b43*k3_[i]); + yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]); } - odes.derivatives(x0 + a4*dx, yTemp_, k4_); + odes.derivatives(x0 + c4*dx, yTemp_, k4_); forAll(yTemp_, i) { yTemp_[i] = y0[i] - + dx*(b51*dydx0[i] + b52*k2_[i] + b53*k3_[i] + b54*k4_[i]); + + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]); } - odes.derivatives(x0 + a5*dx, yTemp_, k5_); + odes.derivatives(x0 + c5*dx, yTemp_, k5_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx - *(b61*dydx0[i] + b62*k2_[i] + b63*k3_[i] + b64*k4_[i] + b65*k5_[i]); + *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]); } - odes.derivatives(x0 + a6*dx, yTemp_, k6_); + odes.derivatives(x0 + c6*dx, yTemp_, k6_); forAll(y, i) { y[i] = y0[i] - + dx*(c1*dydx0[i] + c3*k3_[i] + c4*k4_[i] + c6*k6_[i]); + + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b6*k6_[i]); } forAll(err_, i) { err_[i] = dx - *(dc1*dydx0[i] + dc3*k3_[i] + dc4*k4_[i] + dc5*k5_[i] + dc6*k6_[i]); + *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i]); } - return normalizeError(eps, y0, y, err_); + return normalizeError(y0, y, err_); } -void Foam::RKCK5::solve +void Foam::RKCK45::solve ( const ODESystem& odes, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const { - adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext); + adaptiveSolver::solve(odes, x, y, dxTry); } diff --git a/src/ODE/ODESolvers/RKCK5/RKCK5.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H similarity index 81% rename from src/ODE/ODESolvers/RKCK5/RKCK5.H rename to src/ODE/ODESolvers/RKCK45/RKCK45.H index 73a33e36bd..830e9cfc20 100644 --- a/src/ODE/ODESolvers/RKCK5/RKCK5.H +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -22,10 +22,10 @@ License along with OpenFOAM. If not, see . Class - Foam::RKCK5 + Foam::RKCK45 Description - 5th Order Cash-Karp Runge-Kutta ODE solver + 4/5th Order Cash-Karp Runge-Kutta ODE solver References: \verbatim @@ -37,12 +37,12 @@ Description \endverbatim SourceFiles - RKCK5.C + RKCK45.C \*---------------------------------------------------------------------------*/ -#ifndef RKCK5_H -#define RKCK5_H +#ifndef RKCK45_H +#define RKCK45_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -53,10 +53,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class RKCK5 Declaration + Class RKCK45 Declaration \*---------------------------------------------------------------------------*/ -class RKCK5 +class RKCK45 : public ODESolver, public adaptiveSolver @@ -65,11 +65,11 @@ class RKCK5 //- RKCK Constants static const scalar - a2, a3, a4, a5, a6, - b21, b31, b32, b41, b42, b43, - b51, b52, b53, b54, b61, b62, b63, b64, b65, - c1, c3, c4, c6, - dc1, dc3, dc4, dc5, dc6; + c2, c3, c4, c5, c6, + a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, + a61, a62, a63, a64, a65, + b1, b3, b4, b6, + e1, e3, e4, e5, e6; // Temporary fields mutable scalarField yTemp_; @@ -86,13 +86,13 @@ class RKCK5 public: //- Runtime type information - TypeName("RKCK5"); + TypeName("RKCK45"); // Constructors //- Construct from ODE - RKCK5(const ODESystem& ode); + RKCK45(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -105,8 +105,7 @@ public: const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const; //- Solve the ODE system and the update the state @@ -115,11 +114,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/RKF45/RKF45.C b/src/ODE/ODESolvers/RKF45/RKF45.C new file mode 100644 index 0000000000..b5324e049b --- /dev/null +++ b/src/ODE/ODESolvers/RKF45/RKF45.C @@ -0,0 +1,176 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "RKF45.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(RKF45, 0); + addToRunTimeSelectionTable(ODESolver, RKF45, dictionary); + +const scalar + RKF45::c2 = 1.0/4.0, + RKF45::c3 = 3.0/8.0, + RKF45::c4 = 12.0/13.0, + RKF45::c5 = 1.0, + RKF45::c6 = 1.0/2.0, + + RKF45::a21 = 1.0/4.0, + RKF45::a31 = 3.0/32.0, + RKF45::a32 = 9.0/32.0, + RKF45::a41 = 1932.0/2197.0, + RKF45::a42 = -7200.0/2197.0, + RKF45::a43 = 7296.0/2197.0, + RKF45::a51 = 439.0/216.0, + RKF45::a52 = -8.0, + RKF45::a53 = 3680.0/513.0, + RKF45::a54 = -845.0/4104.0, + RKF45::a61 = -8.0/27.0, + RKF45::a62 = 2.0, + RKF45::a63 = -3544.0/2565.0, + RKF45::a64 = 1859.0/4104.0, + RKF45::a65 = -11.0/40.0, + + RKF45::b41 = 25.0/216.0, + RKF45::b42 = 0.0, + RKF45::b43 = 1408.0/2565.0, + RKF45::b44 = 2197.0/4104.0, + RKF45::b45 = -1.0/5.0, + RKF45::b46 = 0.0, + + RKF45::b51 = 16.0/135.0, + RKF45::b52 = 0.0, + RKF45::b53 = 6656.0/12825.0, + RKF45::b54 = 28561.0/56430.0, + RKF45::b55 = -9.0/50.0, + RKF45::b56 = 2.0/55.0; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + yTemp_(n_), + k2_(n_), + k3_(n_), + k4_(n_), + k5_(n_), + k6_(n_), + err_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::RKF45::solve +( + const ODESystem& odes, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + forAll(yTemp_, i) + { + yTemp_[i] = y0[i] + a21*dx*dydx0[i]; + } + + odes.derivatives(x0 + c2*dx, yTemp_, k2_); + + forAll(yTemp_, i) + { + yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]); + } + + odes.derivatives(x0 + c3*dx, yTemp_, k3_); + + forAll(yTemp_, i) + { + yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]); + } + + odes.derivatives(x0 + c4*dx, yTemp_, k4_); + + forAll(yTemp_, i) + { + yTemp_[i] = y0[i] + + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]); + } + + odes.derivatives(x0 + c5*dx, yTemp_, k5_); + + forAll(yTemp_, i) + { + yTemp_[i] = y0[i] + + dx + *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]); + } + + odes.derivatives(x0 + c6*dx, yTemp_, k6_); + + // Calculate the 5th-order solution + forAll(y, i) + { + y[i] = y0[i] + + dx + *(b51*dydx0[i] + b53*k3_[i] + b54*k4_[i] + b55*k5_[i] + b56*k6_[i]); + } + + // Calculate the error estimate from the difference between the + // 4th-order and 5th-order solutions + forAll(err_, i) + { + err_[i] = + dx + *(b41*dydx0[i] + b43*k3_[i] + b44*k4_[i] + b45*k5_[i]) + - (y[i] - y0[i]) + ; + } + + return normalizeError(y0, y, err_); +} + + +void Foam::RKF45::solve +( + const ODESystem& odes, + scalar& x, + scalarField& y, + scalar& dxTry +) const +{ + adaptiveSolver::solve(odes, x, y, dxTry); +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RKF45/RKF45.H b/src/ODE/ODESolvers/RKF45/RKF45.H new file mode 100644 index 0000000000..34a53d0fb0 --- /dev/null +++ b/src/ODE/ODESolvers/RKF45/RKF45.H @@ -0,0 +1,175 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::RKF45 + +Description + 4/5th Order Runge-Kutta-Fehlberg ODE solver + + References: + \verbatim + "Low-order classical Runge-Kutta formulas with step size control + and their application to some heat transfer problems." + Fehlberg, E., + NASA Technical Report 315 1969 + + "Solving Ordinary Differential Equations I: Nonstiff Problems, + second edition", + Hairer, E., + Nørsett, S., + Wanner, G., + Springer-Verlag, Berlin. 1993, ISBN 3-540-56670-8. + \endverbatim + + This method embedds the 4-th order integration step into the 5-th order step + and allows to perform an adapdive step-size control using these two order + without the need of re-evaluation. + + \f{align}{ + k_1 &= h f(t_k, y_k) \\ + k_2 &= h f(t_k + \frac{1}{4}h, y_k + \frac{1}{4}k_1) \\ + k_3 &= h f(t_k + \frac{3}{8}h, y_k + \frac{3}{32}k_1 + \frac{9}{32}k_2) \\ + k_4 &= h f(t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}k_1 + \frac{7200}{2197}k_2 + \frac{7296}{2197}k_3) - \\ + k_5 &= h f(t_k + h, y_k + \frac{439}{216}k_1 - 8k_2 + \frac{3680}{513}k_3 - + \frac{845}{4104}k_4) \\ + k_6 &= h f(t_k + \frac{1}{2}h, y_k - \frac{8}{27}k_1 + 2k_2 - + \frac{3544}{2565}k_3 + \frac{1859}{4104}k_4 - \frac{11}{40}k_5) \\ + \f} + + Which yields the following update-steps for the 4th and 5th order: + + \f{align}{ + \Delta_4 &= \frac{25}{216}k_1 + \frac{1408}{2565}k_3 + + \frac{2197}{4101}k_4 - \frac{1}{5}k_5 \\ + \Delta_5 &= \frac{16}{135}k_1 + \frac{6656}{12825}k_3 + + \frac{9}{50}k_5 + \frac{2}{55}k_6 + \f} + + The difference between the 5th and 4th order (\f$\epsilon = \left|\Delta_5 + - \Delta_4\right|\f$) can be used to determine if the stepsize \f$h\f$ needs + to be adjusted. The step-size does not need to be reduced if the following + inequation holds true: + + \f[ + \left|\Delta_5 - \Delta_4\right| < \epsilon_{abs} + + \epsilon_{rel}\left|y\right| + \f] + + where \f$\left|y\right|\f$ is the element wise maximum of + \f$\left|y_k\right|\f$ and \f$\left|y_{k+1}\right|\f$. + +SourceFiles + RKF45.C + +\*---------------------------------------------------------------------------*/ + +#ifndef RKF45_H +#define RKF45_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class RKF45 Declaration +\*---------------------------------------------------------------------------*/ + +class RKF45 +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + //- RKF45 Constants + static const scalar + c2, c3, c4, c5, c6, + a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, + a61, a62, a63, a64, a65, + b41, b42, b43, b44, b45, b46, + b51, b52, b53, b54, b55, b56; + + + // Temporary fields + mutable scalarField yTemp_; + mutable scalarField k2_; + mutable scalarField k3_; + mutable scalarField k4_; + mutable scalarField k5_; + mutable scalarField k6_; + + //- Error-estimate field + mutable scalarField err_; + + +public: + + //- Runtime type information + TypeName("RKF45"); + + + // Constructors + + //- Construct from ODE + RKF45(const ODESystem& ode, const dictionary& dict); + + + // Member Functions + + //- Solve a single step dx and return the error + scalar solve + ( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y + ) const; + + //- Solve the ODE system and the update the state + void solve + ( + const ODESystem& ode, + scalar& x, + scalarField& y, + scalar& dxTry + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C index 964e288949..28d4706792 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.C +++ b/src/ODE/ODESolvers/SIBS/SIBS.C @@ -31,8 +31,7 @@ License namespace Foam { defineTypeNameAndDebug(SIBS, 0); - - addToRunTimeSelectionTable(ODESolver, SIBS, ODE); + addToRunTimeSelectionTable(ODESolver, SIBS, dictionary); const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70}; @@ -47,9 +46,9 @@ namespace Foam // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::SIBS::SIBS(const ODESystem& ode) +Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict) : - ODESolver(ode), + ODESolver(ode, dict), a_(iMaxX_, 0.0), alpha_(kMaxX_, kMaxX_, 0.0), d_p_(n_, kMaxX_, 0.0), @@ -59,6 +58,7 @@ Foam::SIBS::SIBS(const ODESystem& ode) yTemp_(n_, 0.0), ySeq_(n_, 0.0), yErr_(n_, 0.0), + dydx0_(n_), dfdx_(n_, 0.0), dfdy_(n_, n_, 0.0), first_(1), @@ -73,19 +73,18 @@ void Foam::SIBS::solve const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const { + ode.derivatives(x, y, dydx0_); + + scalar h = dxTry; bool exitflag = false; - if (eps != epsOld_) + if (relTol_[0] != epsOld_) { - hNext = xNew_ = -GREAT; - scalar eps1 = safe1*eps; + dxTry = xNew_ = -GREAT; + scalar eps1 = safe1*relTol_[0]; a_[0] = nSeq_[0] + 1; for (register label k=0; k= k && kOpt_ != kMax_ && !reduct) { scalar fact = max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX); if (a_[kOpt_ + 1]*fact <= wrkmin) { - hNext = h/fact; + dxTry = h/fact; kOpt_++; } } diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H index 1969085bad..03664fd53a 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.H +++ b/src/ODE/ODESolvers/SIBS/SIBS.H @@ -72,6 +72,7 @@ class SIBS mutable scalarField yTemp_; mutable scalarField ySeq_; mutable scalarField yErr_; + mutable scalarField dydx0_; mutable scalarField dfdx_; mutable scalarSquareMatrix dfdy_; @@ -115,7 +116,7 @@ public: // Constructors //- Construct from ODE - SIBS(const ODESystem& ode); + SIBS(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -125,11 +126,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C index 185b195158..34cc9c347b 100644 --- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C @@ -26,79 +26,50 @@ License #include "adaptiveSolver.H" #include "addToRunTimeSelectionTable.H" -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - -const scalar - adaptiveSolver::safeScale=0.9, - adaptiveSolver::alphaInc=0.2, - adaptiveSolver::alphaDec=0.25, - adaptiveSolver::minScale=0.2, - adaptiveSolver::maxScale=10; -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::adaptiveSolver::adaptiveSolver(const ODESystem& ode) +Foam::adaptiveSolver::adaptiveSolver +( + const ODESystem& ode, + const dictionary& dict +) : + safeScale_(dict.lookupOrDefault("safeScale", 0.9)), + alphaInc_(dict.lookupOrDefault("alphaIncrease", 0.2)), + alphaDec_(dict.lookupOrDefault("alphaDecrease", 0.25)), + minScale_(dict.lookupOrDefault("minScale", 0.2)), + maxScale_(dict.lookupOrDefault("maxScale", 10)), + dydx0_(ode.nEqns()), yTemp_(ode.nEqns()) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::adaptiveSolver::normalizeError -( - const scalar eps, - const scalarField& y0, - const scalarField& y, - const scalarField& err -) const -{ - scalar epsAbs = 1e-10; - scalar epsRel = eps; - - // Calculate the maximum error - scalar maxErr = 0.0; - forAll(err, i) - { - scalar eps = epsAbs + epsRel*max(mag(y0[i]), mag(y[i])); - maxErr = max(maxErr, mag(err[i])/eps); - } - - return maxErr; -} - - void Foam::adaptiveSolver::solve ( const ODESystem& odes, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const { scalar dx = dxTry; scalar err = 0.0; + odes.derivatives(x, y, dydx0_); + // Loop over solver and adjust step-size as necessary // to achieve desired error do { // Solve step and provide error estimate - err = solve(odes, x, y, dydx, dx, yTemp_, eps); + err = solve(odes, x, y, dydx0_, dx, yTemp_); // If error is large reduce dx if (err > 1) { - scalar scale = max(safeScale*pow(err, -alphaDec), minScale); + scalar scale = max(safeScale_*pow(err, -alphaDec_), minScale_); dx *= scale; if (dx < VSMALL) @@ -111,12 +82,11 @@ void Foam::adaptiveSolver::solve } while (err > 1); // Update the state - dxDid = dx; x += dx; y = yTemp_; // If the error is small increase the step-size - dxNext = min(max(safeScale*pow(err, -alphaInc), minScale), maxScale)*dx; + dxTry = min(max(safeScale_*pow(err, -alphaInc_), minScale_), maxScale_)*dx; } diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H index 40f07c2d93..d86bb88ce1 100644 --- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,16 +25,6 @@ Class Foam::adaptiveSolver Description - 5th Order Cash-Karp Runge-Kutta ODE solver - - References: - \verbatim - "A variable order Runge-Kutta method for initial value problems - with rapidly varying right-hand sides" - Cash, J.R., - Karp, A.H. - ACM Transactions on Mathematical Software, vol. 16, 1990, pp. 201–222. - \endverbatim SourceFiles adaptiveSolver.C @@ -60,8 +50,12 @@ class adaptiveSolver // Private data //- Step-size adjustment controls - static const scalar safeScale, alphaInc, alphaDec, minScale, maxScale; + scalar safeScale_, alphaInc_, alphaDec_, minScale_, maxScale_; + //- Cache for dydx at the initial time + mutable scalarField dydx0_; + + //- Temprorary for the test-step solution mutable scalarField yTemp_; @@ -70,7 +64,7 @@ public: // Constructors //- Construct from ODESystem - adaptiveSolver(const ODESystem& ode); + adaptiveSolver(const ODESystem& ode, const dictionary& dict); //- Destructor @@ -80,15 +74,6 @@ public: // Member Functions - //- Return the nomalized scalar error - scalar normalizeError - ( - const scalar eps, - const scalarField& y0, - const scalarField& y, - const scalarField& err - ) const; - //- Solve a single step dx and return the error virtual scalar solve ( @@ -97,8 +82,7 @@ public: const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const = 0; //- Solve the ODE system and the update the state @@ -107,11 +91,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const; }; diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C index b49017a4c8..afb8683086 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C @@ -36,9 +36,7 @@ Foam::ode::ode : chemistrySolver(mesh), coeffsDict_(this->subDict("odeCoeffs")), - solverName_(coeffsDict_.lookup("solver")), - odeSolver_(ODESolver::New(solverName_, *this)), - eps_(readScalar(coeffsDict_.lookup("eps"))), + odeSolver_(ODESolver::New(*this, coeffsDict_)), cTp_(this->nEqns()) {} @@ -78,7 +76,6 @@ void Foam::ode::solve 0, deltaT, cTp_, - eps_, subDeltaT ); diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H index fd1f881112..efdfd40cc2 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.H @@ -55,13 +55,8 @@ class ode // Private data dictionary coeffsDict_; - const word solverName_; autoPtr odeSolver_; - // Model constants - - scalar eps_; - // Solver data mutable scalarField cTp_;