diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C index e5d3026828..4222fa4cc4 100644 --- a/applications/test/ODE/Test-ODE.C +++ b/applications/test/ODE/Test-ODE.C @@ -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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -29,7 +29,6 @@ Description #include "IOmanip.H" #include "ODESystem.H" #include "ODESolver.H" -#include "RK.H" using namespace Foam; @@ -110,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; @@ -125,27 +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); - scalarField yScale(ode.nEqns(), 1.0); - scalar hEst = 0.6; - scalar hDid, hNext; + scalar dxEst = 0.6; + scalar dxNext = dxEst; - odeSolver->solve(ode, x, y, dydx, eps, yScale, 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; @@ -161,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 0e6c49a915..a22ec1d3be 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -1,19 +1,18 @@ -ODE = ODE -ODESolvers = ODESolvers -ODESolversODESolver = ODESolvers/ODESolver -ODESolversKRR4 = ODESolvers/KRR4 -ODESolversRK = ODESolvers/RK -ODESolversSIBS = ODESolvers/SIBS +ODESolvers/ODESolver/ODESolver.C +ODESolvers/ODESolver/ODESolverNew.C -$(ODESolversODESolver)/ODESolver.C -$(ODESolversODESolver)/ODESolverNew.C - -$(ODESolversRK)/RK.C - -$(ODESolversKRR4)/KRR4.C - -$(ODESolversSIBS)/SIBS.C -$(ODESolversSIBS)/SIMPR.C -$(ODESolversSIBS)/polyExtrapolate.C +ODESolvers/adaptiveSolver/adaptiveSolver.C +ODESolvers/Euler/Euler.C +ODESolvers/EulerSI/EulerSI.C +ODESolvers/Trapezoid/Trapezoid.C +ODESolvers/RKF45/RKF45.C +ODESolvers/RKCK45/RKCK45.C +ODESolvers/RKDP45/RKDP45.C +ODESolvers/Rosenbrock21/Rosenbrock21.C +ODESolvers/Rosenbrock43/Rosenbrock43.C +ODESolvers/rodas43/rodas43.C +ODESolvers/SIBS/SIBS.C +ODESolvers/SIBS/SIMPR.C +ODESolvers/SIBS/polyExtrapolate.C LIB = $(FOAM_LIBBIN)/libODE diff --git a/src/ODE/ODESolvers/Euler/Euler.C b/src/ODE/ODESolvers/Euler/Euler.C new file mode 100644 index 0000000000..1a6edf576c --- /dev/null +++ b/src/ODE/ODESolvers/Euler/Euler.C @@ -0,0 +1,88 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "Euler.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(Euler, 0); + addToRunTimeSelectionTable(ODESolver, Euler, dictionary); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + err_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::Euler::solve +( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + // Calculate error estimate from the change in state: + forAll(err_, i) + { + err_[i] = dx*dydx0[i]; + } + + // Update the state + forAll(y, i) + { + y[i] = y0[i] + err_[i]; + } + + return normalizeError(y0, y, err_); +} + + +void Foam::Euler::solve +( + const ODESystem& odes, + scalar& x, + scalarField& y, + scalar& dxTry +) const +{ + adaptiveSolver::solve(odes, x, y, dxTry); +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/Euler/Euler.H b/src/ODE/ODESolvers/Euler/Euler.H new file mode 100644 index 0000000000..90f67f4b4c --- /dev/null +++ b/src/ODE/ODESolvers/Euler/Euler.H @@ -0,0 +1,114 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::Euler + +Description + Euler ODE solver of order (0)1. + + The method calculates the new state from: + \f[ + y_{n+1} = y_n + \delta_x f + \f] + The error is estimated directly from the change in the solution, + i.e. the difference between the 0th and 1st order solutions: + \f[ + err_{n+1} = y_{n+1} - y_n + \f] + +SourceFiles + Euler.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Euler_H +#define Euler_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Euler Declaration +\*---------------------------------------------------------------------------*/ + +class Euler +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + mutable scalarField err_; + + +public: + + //- Runtime type information + TypeName("Euler"); + + + // Constructors + + //- Construct from ODE + Euler(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/EulerSI/EulerSI.C b/src/ODE/ODESolvers/EulerSI/EulerSI.C new file mode 100644 index 0000000000..1785101758 --- /dev/null +++ b/src/ODE/ODESolvers/EulerSI/EulerSI.C @@ -0,0 +1,108 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "EulerSI.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(EulerSI, 0); + addToRunTimeSelectionTable(ODESolver, EulerSI, dictionary); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + err_(n_), + dydx_(n_), + dfdx_(n_), + dfdy_(n_, n_), + a_(n_, n_), + pivotIndices_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::EulerSI::solve +( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + ode.jacobian(x0, y0, dfdx_, dfdy_); + + for (register label i=0; i. Class - Foam::KRR4 + Foam::EulerSI Description - Foam::KRR4 + Semi-implicit Euler ODE solver of order (0)1. + + The method calculates the new state from: + \f[ + y_{n+1} = y_n + + \delta_x\left[I - \delta_x\frac{\partial f}{\partial y}\right]^{-1} + \cdot \left[f(y_n) + \delta_x\frac{\partial f}{\partial x}\right] + \f] + The error is estimated directly from the change in the solution, + i.e. the difference between the 0th and 1st order solutions: + \f[ + err_{n+1} = y_{n+1} - y_n + \f] SourceFiles - KRR4CK.C - KRR4QS.C + EulerSI.C \*---------------------------------------------------------------------------*/ -#ifndef KRR4_H -#define KRR4_H +#ifndef EulerSI_H +#define EulerSI_H #include "ODESolver.H" +#include "adaptiveSolver.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -44,66 +56,56 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class KRR4 Declaration + Class EulerSI Declaration \*---------------------------------------------------------------------------*/ -class KRR4 +class EulerSI : - public ODESolver + public ODESolver, + public adaptiveSolver { // Private data - mutable scalarField yTemp_; - mutable scalarField dydxTemp_; - mutable scalarField g1_; - mutable scalarField g2_; - mutable scalarField g3_; - mutable scalarField g4_; - mutable scalarField yErr_; + mutable scalarField err_; + mutable scalarField dydx_; mutable scalarField dfdx_; mutable scalarSquareMatrix dfdy_; mutable scalarSquareMatrix a_; mutable labelList pivotIndices_; - static const int maxtry = 40; - - static const scalar safety, grow, pgrow, shrink, pshrink, errcon; - - static const scalar - gamma, - a21, a31, a32, - c21, c31, c32, c41, c42, c43, - b1, b2, b3, b4, - e1, e2, e3, e4, - c1X, c2X, c3X, c4X, - a2X, a3X; - public: //- Runtime type information - TypeName("KRR4"); + TypeName("EulerSI"); // Constructors //- Construct from ODE - KRR4(const ODESystem& ode); + EulerSI(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, - scalarField& dydx, - const scalar eps, - const scalarField& yScale, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/KRR4/KRR4.C b/src/ODE/ODESolvers/KRR4/KRR4.C deleted file mode 100644 index e9218641a7..0000000000 --- a/src/ODE/ODESolvers/KRR4/KRR4.C +++ /dev/null @@ -1,224 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 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 "KRR4.H" -#include "addToRunTimeSelectionTable.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(KRR4, 0); - - addToRunTimeSelectionTable(ODESolver, KRR4, ODE); - -const scalar - KRR4::safety = 0.9, KRR4::grow = 1.5, KRR4::pgrow = -0.25, - KRR4::shrink = 0.5, KRR4::pshrink = (-1.0/3.0), KRR4::errcon = 0.1296; - -const scalar - KRR4::gamma = 1.0/2.0, - KRR4::a21 = 2.0, KRR4::a31 = 48.0/25.0, KRR4::a32 = 6.0/25.0, - KRR4::c21 = -8.0, KRR4::c31 = 372.0/25.0, KRR4::c32 = 12.0/5.0, - KRR4::c41 = -112.0/125.0, KRR4::c42 = -54.0/125.0, KRR4::c43 = -2.0/5.0, - KRR4::b1 = 19.0/9.0, KRR4::b2 = 1.0/2.0, KRR4::b3 = 25.0/108.0, - KRR4::b4 = 125.0/108.0, - KRR4::e1 = 17.0/54.0, KRR4::e2 = 7.0/36.0, KRR4::e3 = 0.0, - KRR4::e4 = 125.0/108.0, - KRR4::c1X = 1.0/2.0, KRR4::c2X = -3.0/2.0, KRR4::c3X = 121.0/50.0, - KRR4::c4X = 29.0/250.0, - KRR4::a2X = 1.0, KRR4::a3X = 3.0/5.0; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::KRR4::KRR4(const ODESystem& ode) -: - ODESolver(ode), - yTemp_(n_, 0.0), - dydxTemp_(n_, 0.0), - g1_(n_, 0.0), - g2_(n_, 0.0), - g3_(n_, 0.0), - g4_(n_, 0.0), - yErr_(n_, 0.0), - dfdx_(n_, 0.0), - dfdy_(n_, n_, 0.0), - a_(n_, n_, 0.0), - pivotIndices_(n_, 0.0) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::KRR4::solve -( - const ODESystem& ode, - scalar& x, - scalarField& y, - scalarField& dydx, - const scalar eps, - const scalarField& yScale, - const scalar hTry, - scalar& hDid, - scalar& hNext -) const -{ - scalar xTemp = x; - yTemp_ = y; - dydxTemp_ = dydx; - - ode.jacobian(xTemp, yTemp_, dfdx_, dfdy_); - - scalar h = hTry; - - for (register label jtry=0; jtry errcon ? safety*h*pow(maxErr, pgrow) : grow*h); - return; - } - else - { - hNext = safety*h*pow(maxErr, pshrink); - h = (h >= 0.0 ? max(hNext, shrink*h) : min(hNext, shrink*h)); - } - } - - FatalErrorIn - ( - "void Foam::KRR4::solve" - "(" - "const ODESystem&, " - "scalar&, " - "scalarField&, " - "scalarField&, " - "const scalar, " - "const scalarField&, " - "const scalar, " - "scalar&, " - "scalar&" - ") const" - ) << "Maximum number of solver iterations exceeded" - << exit(FatalError); -} - - -// ************************************************************************* // diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 5f7188a079..f4dd3588eb 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -30,82 +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()), - yScale_(n_), - dydx_(n_) + absTol_(n_, dict.lookupOrDefault("absTol", SMALL)), + relTol_(n_, dict.lookupOrDefault("relTol", 1e-4)), + 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) { - yScale_[i] = mag(y[i]) + mag(dydx_[i]*h) + SMALL; + truncated = true; + dxEst = xEnd - x; } - if ((x + h - xEnd)*(x + h - xStart) > 0.0) - { - h = xEnd - x; - hPrev = hNext; - } + // Integrate as far as possible up to dxEst + solve(ode, x, y, dxEst); - hNext = 0; - scalar hDid; - solve(ode, x, y, dydx_, eps, yScale_, h, hDid, hNext); - - 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 04315ea7dc..bd767a9063 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H @@ -53,14 +53,30 @@ class ODESolver protected: - // Private data + // Protected data + //- Size of the ODESystem label n_; - mutable scalarField yScale_; - 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&); @@ -81,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 @@ -98,8 +122,8 @@ public: //- Select null constructed static autoPtr New ( - const word& ODESolverTypeName, - const ODESystem& ode + const ODESystem& ode, + const dictionary& dict ); @@ -110,28 +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 scalarField& yScale, - 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/RK/RK.C b/src/ODE/ODESolvers/RK/RK.C deleted file mode 100644 index 06c5e6ae81..0000000000 --- a/src/ODE/ODESolvers/RK/RK.C +++ /dev/null @@ -1,204 +0,0 @@ -/*---------------------------------------------------------------------------*\ - ========= | - \\ / F ield | OpenFOAM: The Open Source CFD Toolbox - \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 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 "RK.H" -#include "addToRunTimeSelectionTable.H" - -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - defineTypeNameAndDebug(RK, 0); - addToRunTimeSelectionTable(ODESolver, RK, ODE); - -const scalar - RK::safety=0.9, RK::pGrow=-0.2, RK::pShrink=-0.25, RK::errCon=1.89e-4; - -const scalar - RK::a2 = 0.2, RK::a3 = 0.3, RK::a4 = 0.6, RK::a5 = 1.0, RK::a6 = 0.875, - RK::b21 = 0.2, RK::b31 = 3.0/40.0, RK::b32 = 9.0/40.0, - RK::b41 = 0.3, RK::b42 = -0.9, RK::b43 = 1.2, - RK::b51 = -11.0/54.0, RK::b52 = 2.5, RK::b53 = -70.0/27.0, - RK::b54 = 35.0/27.0, - RK::b61 = 1631.0/55296.0, RK::b62 = 175.0/512.0, RK::b63 = 575.0/13824.0, - RK::b64 = 44275.0/110592.0, RK::b65 = 253.0/4096.0, - RK::c1 = 37.0/378.0, RK::c3 = 250.0/621.0, - RK::c4 = 125.0/594.0, RK::c6 = 512.0/1771.0, - RK::dc1 = RK::c1 - 2825.0/27648.0, RK::dc3 = RK::c3 - 18575.0/48384.0, - RK::dc4 = RK::c4 - 13525.0/55296.0, RK::dc5 = -277.00/14336.0, - RK::dc6 = RK::c6 - 0.25; -} - - -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::RK::RK(const ODESystem& ode) -: - ODESolver(ode), - yTemp_(n_, 0.0), - ak2_(n_, 0.0), - ak3_(n_, 0.0), - ak4_(n_, 0.0), - ak5_(n_, 0.0), - ak6_(n_, 0.0), - yErr_(n_, 0.0), - yTemp2_(n_, 0.0) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::RK::solve -( - const ODESystem& ode, - const scalar x, - const scalarField& y, - const scalarField& dydx, - const scalar h, - scalarField& yout, - scalarField& yerr -) const -{ - forAll(yTemp_, i) - { - yTemp_[i] = y[i] + b21*h*dydx[i]; - } - - ode.derivatives(x + a2*h, yTemp_, ak2_); - - forAll(yTemp_, i) - { - yTemp_[i] = y[i] + h*(b31*dydx[i] + b32*ak2_[i]); - } - - ode.derivatives(x + a3*h, yTemp_, ak3_); - - forAll(yTemp_, i) - { - yTemp_[i] = y[i] + h*(b41*dydx[i] + b42*ak2_[i] + b43*ak3_[i]); - } - - ode.derivatives(x + a4*h, yTemp_, ak4_); - - forAll(yTemp_, i) - { - yTemp_[i] = y[i] - + h*(b51*dydx[i] + b52*ak2_[i] + b53*ak3_[i] + b54*ak4_[i]); - } - - ode.derivatives(x + a5*h, yTemp_, ak5_); - - forAll(yTemp_, i) - { - yTemp_[i] = y[i] - + h* - ( - b61*dydx[i] + b62*ak2_[i] + b63*ak3_[i] - + b64*ak4_[i] + b65*ak5_[i] - ); - } - - ode.derivatives(x + a6*h, yTemp_, ak6_); - - forAll(yout, i) - { - yout[i] = y[i] - + h*(c1*dydx[i] + c3*ak3_[i] + c4*ak4_[i] + c6*ak6_[i]); - } - - forAll(yerr, i) - { - yerr[i] = - h* - ( - dc1*dydx[i] + dc3*ak3_[i] + dc4*ak4_[i] - + dc5*ak5_[i] + dc6*ak6_[i] - ); - } -} - - -void Foam::RK::solve -( - const ODESystem& ode, - scalar& x, - scalarField& y, - scalarField& dydx, - const scalar eps, - const scalarField& yScale, - const scalar hTry, - scalar& hDid, - scalar& hNext -) const -{ - scalar h = hTry; - scalar maxErr = 0.0; - - for (;;) - { - solve(ode, x, y, dydx, h, yTemp2_, yErr_); - - maxErr = 0.0; - for (register label i=0; i= 0.0 ? max(hTemp, 0.1*h) : min(hTemp, 0.1*h)); - } - - if (h < VSMALL) - { - FatalErrorIn("RK::solve") - << "stepsize underflow" - << exit(FatalError); - } - } - - hDid = h; - - x += h; - y = yTemp2_; - - if (maxErr > errCon) - { - hNext = safety*h*pow(maxErr, pGrow); - } - else - { - hNext = 5.0*h; - } -} - - -// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C new file mode 100644 index 0000000000..2de3cdf99e --- /dev/null +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.C @@ -0,0 +1,167 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "RKCK45.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(RKCK45, 0); + addToRunTimeSelectionTable(ODESolver, RKCK45, dictionary); + +const scalar + 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::RKCK45::RKCK45(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::RKCK45::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_); + + forAll(y, i) + { + y[i] = y0[i] + + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b6*k6_[i]); + } + + forAll(err_, i) + { + err_[i] = + dx + *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i]); + } + + return normalizeError(y0, y, err_); +} + + +void Foam::RKCK45::solve +( + const ODESystem& odes, + scalar& x, + scalarField& y, + scalar& dxTry +) const +{ + adaptiveSolver::solve(odes, x, y, dxTry); +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H new file mode 100644 index 0000000000..a6b551a475 --- /dev/null +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.H @@ -0,0 +1,140 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::RKCK45 + +Description + 4/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 + + Based on code from: + \verbatim + "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 + +SourceFiles + RKCK45.C + +\*---------------------------------------------------------------------------*/ + +#ifndef RKCK45_H +#define RKCK45_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class RKCK45 Declaration +\*---------------------------------------------------------------------------*/ + +class RKCK45 +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + //- RKCK Constants + static const scalar + 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_; + 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("RKCK45"); + + + // Constructors + + //- Construct from ODE + RKCK45(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/RKDP45/RKDP45.C b/src/ODE/ODESolvers/RKDP45/RKDP45.C new file mode 100644 index 0000000000..cf319dca8b --- /dev/null +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.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 "RKDP45.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(RKDP45, 0); + addToRunTimeSelectionTable(ODESolver, RKDP45, dictionary); + +const scalar + RKDP45::c2 = 1.0/5.0, + RKDP45::c3 = 3.0/10.0, + RKDP45::c4 = 4.0/5.0, + RKDP45::c5 = 8.0/9.0, + + RKDP45::a21 = 1.0/5.0, + + RKDP45::a31 = 3.0/40.0, + RKDP45::a32 = 9.0/40.0, + + RKDP45::a41 = 44.0/45.0, + RKDP45::a42 = -56.0/15.0, + RKDP45::a43 = 32.0/9.0, + + RKDP45::a51 = 19372.0/6561.0, + RKDP45::a52 = -25360.0/2187.0, + RKDP45::a53 = 64448.0/6561.0, + RKDP45::a54 = -212.0/729.0, + + RKDP45::a61 = 9017.0/3168.0, + RKDP45::a62 = -355.0/33.0, + RKDP45::a63 = 46732.0/5247.0, + RKDP45::a64 = 49.0/176.0, + RKDP45::a65 = -5103.0/18656.0, + + RKDP45::b1 = 35.0/384.0, + RKDP45::b3 = 500.0/1113.0, + RKDP45::b4 = 125.0/192.0, + RKDP45::b5 = -2187.0/6784.0, + RKDP45::b6 = 11.0/84.0, + + RKDP45::e1 = 5179.0/57600.0 - RKDP45::b1, + RKDP45::e3 = 7571.0/16695.0 - RKDP45::b3, + RKDP45::e4 = 393.0/640.0 - RKDP45::b4, + RKDP45::e5 = -92097.0/339200.0 - RKDP45::b5, + RKDP45::e6 = 187.0/2100.0 - RKDP45::b6, + RKDP45::e7 = 1.0/40.0; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::RKDP45::RKDP45(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::RKDP45::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 + dx, yTemp_, k6_); + + forAll(y, i) + { + y[i] = y0[i] + + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b5*k5_[i] + b6*k6_[i]); + } + + // Reuse k2_ for the derivative of the new state + odes.derivatives(x0 + dx, y, k2_); + + forAll(err_, i) + { + err_[i] = + dx + *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i] + + e7*k2_[i]); + } + + return normalizeError(y0, y, err_); +} + + +void Foam::RKDP45::solve +( + const ODESystem& odes, + scalar& x, + scalarField& y, + scalar& dxTry +) const +{ + adaptiveSolver::solve(odes, x, y, dxTry); +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.H b/src/ODE/ODESolvers/RKDP45/RKDP45.H new file mode 100644 index 0000000000..95ae212e66 --- /dev/null +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.H @@ -0,0 +1,144 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::RKDP45 + +Description + 4/5th Order Dormand–Prince Runge-Kutta ODE solver. + + References: + \verbatim + "A family of embedded Runge-Kutta formulae" + Dormand, J. R., + Prince, P. J., + Journal of Computational and Applied Mathematics, + 6 (1), 1980: pp. 19-26. + \endverbatim + + Based on code from: + \verbatim + "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 + +SeeAlso + Foam::RKF45 + Foam::RKCK45 + +SourceFiles + RKDP45.C + +\*---------------------------------------------------------------------------*/ + +#ifndef RKDP45_H +#define RKDP45_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class RKDP45 Declaration +\*---------------------------------------------------------------------------*/ + +class RKDP45 +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + //- RKDP Constants + static const scalar + c2, c3, c4, c5, + a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, + a61, a62, a63, a64, a65, + b1, b3, b4, b5, b6, + e1, e3, e4, e5, e6, e7; + + // 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("RKDP45"); + + + // Constructors + + //- Construct from ODE + RKDP45(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/RKF45/RKF45.C b/src/ODE/ODESolvers/RKF45/RKF45.C new file mode 100644 index 0000000000..661c618255 --- /dev/null +++ b/src/ODE/ODESolvers/RKF45/RKF45.C @@ -0,0 +1,172 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::b1 = 16.0/135.0, + RKF45::b3 = 6656.0/12825.0, + RKF45::b4 = 28561.0/56430.0, + RKF45::b5 = -9.0/50.0, + RKF45::b6 = 2.0/55.0, + + RKF45::e1 = 25.0/216.0 - RKF45::b1, + RKF45::e3 = 1408.0/2565.0 - RKF45::b3, + RKF45::e4 = 2197.0/4104.0 - RKF45::b4, + RKF45::e5 = -1.0/5.0 - RKF45::b5, + RKF45::e6 = -RKF45::b6; +} + + +// * * * * * * * * * * * * * * * * 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 + *(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b5*k5_[i] + b6*k6_[i]); + } + + // Calculate the error estimate from the difference between the + // 4th-order and 5th-order solutions + forAll(err_, i) + { + err_[i] = + dx + *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[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..df91ed7808 --- /dev/null +++ b/src/ODE/ODESolvers/RKF45/RKF45.H @@ -0,0 +1,168 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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. + + Based on code from: + \verbatim + "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. + +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, + b1, b3, b4, b5, b6, + e1, e3, e4, e5, e6; + + + // 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/Rosenbrock21/Rosenbrock21.C b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C new file mode 100644 index 0000000000..de6b968348 --- /dev/null +++ b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C @@ -0,0 +1,139 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "Rosenbrock21.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(Rosenbrock21, 0); + addToRunTimeSelectionTable(ODESolver, Rosenbrock21, dictionary); + +const scalar + Rosenbrock21::gamma = 1 + 1.0/std::sqrt(2.0), + Rosenbrock21::a21 = 1.0/gamma, + Rosenbrock21::c2 = 1.0, + Rosenbrock21::c21 = -2.0/gamma, + Rosenbrock21::b1 = (3.0/2.0)/gamma, + Rosenbrock21::b2 = (1.0/2.0)/gamma, + Rosenbrock21::e1 = b1 - 1.0/gamma, + Rosenbrock21::e2 = b2, + Rosenbrock21::d1 = 1, + Rosenbrock21::d2 = -1; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + k1_(n_), + k2_(n_), + err_(n_), + dydx_(n_), + dfdx_(n_), + dfdy_(n_, n_), + a_(n_, n_), + pivotIndices_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::Rosenbrock21::solve +( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + ode.jacobian(x0, y0, dfdx_, dfdy_); + + for (register label i=0; i. + +Class + Foam::Rosenbrock21 + +Description + L-stable embedded Rosenbrock ODE solver of order (1)2. + + References: + \verbatim + "A second-order Rosenbrock method applied to + photochemical dispersion problems", + J. G. Verwer, + E. J. Spee, + J. G. Blom, + W. Hundsdorfer, + Siam Journal on Scientific Computing 01/1999; 20(4):1456-1480. + \endverbatim + +SourceFiles + Rosenbrock21.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Rosenbrock21_H +#define Rosenbrock21_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Rosenbrock21 Declaration +\*---------------------------------------------------------------------------*/ + +class Rosenbrock21 +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + mutable scalarField k1_; + mutable scalarField k2_; + mutable scalarField err_; + mutable scalarField dydx_; + mutable scalarField dfdx_; + mutable scalarSquareMatrix dfdy_; + mutable scalarSquareMatrix a_; + mutable labelList pivotIndices_; + + static const scalar + a21, + c21, + b1, b2, + gamma, + c2, + e1, e2, + d1, d2; + + +public: + + //- Runtime type information + TypeName("Rosenbrock21"); + + + // Constructors + + //- Construct from ODE + Rosenbrock21(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/Rosenbrock43/Rosenbrock43.C b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C new file mode 100644 index 0000000000..04b55dbbe6 --- /dev/null +++ b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C @@ -0,0 +1,221 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "Rosenbrock43.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(Rosenbrock43, 0); + addToRunTimeSelectionTable(ODESolver, Rosenbrock43, dictionary); + +const scalar + // L-Stable constants from Hairer et. al. + Rosenbrock43::a21 = 2, + Rosenbrock43::a31 = 1.867943637803922, + Rosenbrock43::a32 = 0.2344449711399156, + + Rosenbrock43::c21 = -7.137615036412310, + Rosenbrock43::c31 = 2.580708087951457, + Rosenbrock43::c32 = 0.6515950076447975, + Rosenbrock43::c41 = -2.137148994382534, + Rosenbrock43::c42 = -0.3214669691237626, + Rosenbrock43::c43 = -0.6949742501781779, + + Rosenbrock43::b1 = 2.255570073418735, + Rosenbrock43::b2 = 0.2870493262186792, + Rosenbrock43::b3 = 0.435317943184018, + Rosenbrock43::b4 = 1.093502252409163, + + Rosenbrock43::e1 = -0.2815431932141155, + Rosenbrock43::e2 = -0.0727619912493892, + Rosenbrock43::e3 = -0.1082196201495311, + Rosenbrock43::e4 = -1.093502252409163, + + Rosenbrock43::gamma = 0.57282, + Rosenbrock43::c2 = 1.14564, + Rosenbrock43::c3 = 0.65521686381559, + + Rosenbrock43::d1 = 0.57282, + Rosenbrock43::d2 = -1.769193891319233, + Rosenbrock43::d3 = 0.7592633437920482, + Rosenbrock43::d4 = -0.1049021087100450; + + // Constants by Shampine + // More accurate than the L-Stable coefficients for small step-size + // but less stable for large step-size + /* + Rosenbrock43::a21 = 2, + Rosenbrock43::a31 = 48.0/25.0, + Rosenbrock43::a32 = 6.0/25.0, + + Rosenbrock43::c21 = -8, + Rosenbrock43::c31 = 372.0/25.0, + Rosenbrock43::c32 = 12.0/5.0, + + Rosenbrock43::c41 = -112.0/125.0, + Rosenbrock43::c42 = -54.0/125.0, + Rosenbrock43::c43 = -2.0/5.0, + + Rosenbrock43::b1 = 19.0/9.0, + Rosenbrock43::b2 = 1.0/2.0, + Rosenbrock43::b3 = 25.0/108.0, + Rosenbrock43::b4 = 125.0/108.0, + + Rosenbrock43::e1 = 34.0/108.0, + Rosenbrock43::e2 = 7.0/36.0, + Rosenbrock43::e3 = 0, + Rosenbrock43::e4 = 125.0/108.0, + + Rosenbrock43::gamma = 0.5, + Rosenbrock43::c2 = 1, + Rosenbrock43::c3 = 3.0/5.0, + + Rosenbrock43::d1 = 1.0/2.0, + Rosenbrock43::d2 = -3.0/2.0, + Rosenbrock43::d3 = 605.0/250.0, + Rosenbrock43::d4 = 29.0/250.0; + */ +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + k1_(n_), + k2_(n_), + k3_(n_), + k4_(n_), + err_(n_), + dydx_(n_), + dfdx_(n_), + dfdy_(n_, n_), + a_(n_, n_), + pivotIndices_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::Rosenbrock43::solve +( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + ode.jacobian(x0, y0, dfdx_, dfdy_); + + for (register label i=0; i. + +Class + Foam::Rosenbrock43 + +Description + L-stable embedded Rosenbrock ODE solver of order (3)4. + + Based on code from: + \verbatim + "Solving Ordinary Differential Equations II: Stiff + and Differential-Algebraic Problems, second edition", + Hairer, E., + Nørsett, S., + Wanner, G., + Springer-Verlag, Berlin. 1996. + \endverbatim + + Also provided are the constants from: + \verbatim + "Implementation of Rosenbrock Methods" + Shampine, L. F., + ACM Transactions on Mathematical Software, vol. 8, 1982, pp. 93–113. + \endverbatim + with which the scheme is more accurate than with the L-Stable coefficients + for small step-size but less stable for large step-size. + +SourceFiles + Rosenbrock43.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Rosenbrock43_H +#define Rosenbrock43_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Rosenbrock43 Declaration +\*---------------------------------------------------------------------------*/ + +class Rosenbrock43 +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + mutable scalarField k1_; + mutable scalarField k2_; + mutable scalarField k3_; + mutable scalarField k4_; + mutable scalarField err_; + mutable scalarField dydx_; + mutable scalarField dfdx_; + mutable scalarSquareMatrix dfdy_; + mutable scalarSquareMatrix a_; + mutable labelList pivotIndices_; + + static const scalar + a21, a31, a32, + c21, c31, c32, + c41, c42, c43, + b1, b2, b3, b4, + e1, e2, e3, e4, + gamma, + c2, c3, + d1, d2, d3, d4; + + +public: + + //- Runtime type information + TypeName("Rosenbrock43"); + + + // Constructors + + //- Construct from ODE + Rosenbrock43(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 95ab44776c..eae8552a47 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.C +++ b/src/ODE/ODESolvers/SIBS/SIBS.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -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,20 +73,18 @@ void Foam::SIBS::solve const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalarField& yScale, - 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 a031b1c83b..03664fd53a 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.H +++ b/src/ODE/ODESolvers/SIBS/SIBS.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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -27,6 +27,11 @@ Class Description Foam::SIBS + Bader, G. and Deuflhard, P. + "A Semi-Implicit Mid-Point Rule for + Stiff Systems of Ordinary Differential Equations." + Numer. Math. 41, 373-398, 1983. + SourceFiles SIMPR.C polyExtrapolate.C @@ -67,6 +72,7 @@ class SIBS mutable scalarField yTemp_; mutable scalarField ySeq_; mutable scalarField yErr_; + mutable scalarField dydx0_; mutable scalarField dfdx_; mutable scalarSquareMatrix dfdy_; @@ -110,7 +116,7 @@ public: // Constructors //- Construct from ODE - SIBS(const ODESystem& ode); + SIBS(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -120,12 +126,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalarField& yScale, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/Trapezoid/Trapezoid.C b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C new file mode 100644 index 0000000000..8a2c281433 --- /dev/null +++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C @@ -0,0 +1,93 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "Trapezoid.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(Trapezoid, 0); + addToRunTimeSelectionTable(ODESolver, Trapezoid, dictionary); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::Trapezoid::Trapezoid(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + err_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::Trapezoid::solve +( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + // Predict the state using 1st-order Trapezoid method + forAll(y, i) + { + y[i] = y0[i] + dx*dydx0[i]; + } + + // Evaluate the system for the predicted state + ode.derivatives(x0 + dx, y, err_); + + // Update the state as the average between the prediction and the correction + // and estimate the error from the difference + forAll(y, i) + { + y[i] = y0[i] + 0.5*dx*(dydx0[i] + err_[i]); + err_[i] = 0.5*dx*(err_[i] - dydx0[i]); + } + + return normalizeError(y0, y, err_); +} + + +void Foam::Trapezoid::solve +( + const ODESystem& odes, + scalar& x, + scalarField& y, + scalar& dxTry +) const +{ + adaptiveSolver::solve(odes, x, y, dxTry); +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RK/RK.H b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H similarity index 62% rename from src/ODE/ODESolvers/RK/RK.H rename to src/ODE/ODESolvers/Trapezoid/Trapezoid.H index 20fa3aa789..327dbe44bc 100644 --- a/src/ODE/ODESolvers/RK/RK.H +++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.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 @@ -22,21 +22,21 @@ License along with OpenFOAM. If not, see . Class - Foam::RK + Foam::Trapezoid Description - Runge-Kutta ODE solver + Trapezoidal ODE solver of order (1)2. SourceFiles - RKCK.C - RKQS.C + Trapezoid.C \*---------------------------------------------------------------------------*/ -#ifndef RK_H -#define RK_H +#ifndef Trapezoid_H +#define Trapezoid_H #include "ODESolver.H" +#include "adaptiveSolver.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -44,71 +44,51 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class RK Declaration + Class Trapezoid Declaration \*---------------------------------------------------------------------------*/ -class RK +class Trapezoid : - public ODESolver + public ODESolver, + public adaptiveSolver { // Private data - static const scalar safety, pGrow, pShrink, errCon; - 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; + mutable scalarField err_; - mutable scalarField yTemp_; - mutable scalarField ak2_; - mutable scalarField ak3_; - mutable scalarField ak4_; - mutable scalarField ak5_; - mutable scalarField ak6_; - - mutable scalarField yErr_; - mutable scalarField yTemp2_; - public: //- Runtime type information - TypeName("RK"); + TypeName("Trapezoid"); // Constructors //- Construct from ODE - RK(const ODESystem& ode); + Trapezoid(const ODESystem& ode, const dictionary& dict); // Member Functions - void solve + //- Solve a single step dx and return the error + scalar solve ( const ODESystem& ode, - const scalar x, - const scalarField& y, - const scalarField& dydx, - const scalar h, - scalarField& yout, - scalarField& yerr + 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, - scalarField& dydx, - const scalar eps, - const scalarField& yScale, - 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 new file mode 100644 index 0000000000..13951ae7b3 --- /dev/null +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C @@ -0,0 +1,101 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-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 "adaptiveSolver.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +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 * * * * * * * * * * * * * // + +void Foam::adaptiveSolver::solve +( + const ODESystem& odes, + scalar& x, + scalarField& y, + 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, dydx0_, dx, yTemp_); + + // If error is large reduce dx + if (err > 1) + { + scalar scale = max(safeScale_*pow(err, -alphaDec_), minScale_); + dx *= scale; + + if (dx < VSMALL) + { + FatalErrorIn("adaptiveSolver::solve") + << "stepsize underflow" + << exit(FatalError); + } + } + } while (err > 1); + + // Update the state + x += dx; + y = yTemp_; + + // If the error is small increase the step-size + if (err > pow(maxScale_/safeScale_, -1.0/alphaInc_)) + { + dxTry = + min(max(safeScale_*pow(err, -alphaInc_), minScale_), maxScale_)*dx; + } + else + { + dxTry = safeScale_*maxScale_*dx; + } +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H new file mode 100644 index 0000000000..d86bb88ce1 --- /dev/null +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::adaptiveSolver + +Description + +SourceFiles + adaptiveSolver.C + +\*---------------------------------------------------------------------------*/ + +#ifndef adaptiveSolver_H +#define adaptiveSolver_H + +#include "ODESolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class adaptiveSolver Declaration +\*---------------------------------------------------------------------------*/ + +class adaptiveSolver +{ + // Private data + + //- Step-size adjustment controls + 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_; + + +public: + + // Constructors + + //- Construct from ODESystem + adaptiveSolver(const ODESystem& ode, const dictionary& dict); + + + //- Destructor + virtual ~adaptiveSolver() + {} + + + // Member Functions + + //- Solve a single step dx and return the error + virtual scalar solve + ( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y + ) const = 0; + + //- 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/rodas43/rodas43.C b/src/ODE/ODESolvers/rodas43/rodas43.C new file mode 100644 index 0000000000..74a3da1d80 --- /dev/null +++ b/src/ODE/ODESolvers/rodas43/rodas43.C @@ -0,0 +1,232 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 "rodas43.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(rodas43, 0); + addToRunTimeSelectionTable(ODESolver, rodas43, dictionary); + +const scalar + rodas43::c2 = 0.386, + rodas43::c3 = 0.21, + rodas43::c4 = 0.63, + rodas43::bet2p = 0.0317, + rodas43::bet3p = 0.0635, + rodas43::bet4p = 0.3438, + rodas43::d1 = 0.25, + rodas43::d2 = -0.1043, + rodas43::d3 = 0.1035, + rodas43::d4 = -0.3620000000000023e-01, + rodas43::a21 = 0.1544e1, + rodas43::a31 = 0.9466785280815826, + rodas43::a32 = 0.2557011698983284, + rodas43::a41 = 0.3314825187068521e1, + rodas43::a42 = 0.2896124015972201e1, + rodas43::a43 = 0.9986419139977817, + rodas43::a51 = 0.1221224509226641e1, + rodas43::a52 = 0.6019134481288629e1, + rodas43::a53 = 0.1253708332932087e2, + rodas43::a54 = -0.6878860361058950, + rodas43::c21 = -0.56688e1, + rodas43::c31 = -0.2430093356833875e1, + rodas43::c32 = -0.2063599157091915, + rodas43::c41 = -0.1073529058151375, + rodas43::c42 = -0.9594562251023355e1, + rodas43::c43 = -0.2047028614809616e2, + rodas43::c51 = 0.7496443313967647e1, + rodas43::c52 = -0.1024680431464352e2, + rodas43::c53 = -0.3399990352819905e2, + rodas43::c54 = 0.1170890893206160e2, + rodas43::c61 = 0.8083246795921522e1, + rodas43::c62 = -0.7981132988064893e1, + rodas43::c63 = -0.3152159432874371e2, + rodas43::c64 = 0.1631930543123136e2, + rodas43::c65 = -0.6058818238834054e1, + rodas43::gamma = 0.25; +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict) +: + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + k1_(n_), + k2_(n_), + k3_(n_), + k4_(n_), + k5_(n_), + dy_(n_), + err_(n_), + dydx_(n_), + dfdx_(n_), + dfdy_(n_, n_), + a_(n_, n_), + pivotIndices_(n_) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::scalar Foam::rodas43::solve +( + const ODESystem& ode, + const scalar x0, + const scalarField& y0, + const scalarField& dydx0, + const scalar dx, + scalarField& y +) const +{ + ode.jacobian(x0, y0, dfdx_, dfdy_); + + for (register label i=0; i. + +Class + Foam::rodas43 + +Description + Stiffly-stable embedded Rosenbrock ODE solver of order (3)4. + + References: + \verbatim + "Solving Ordinary Differential Equations II: Stiff + and Differential-Algebraic Problems, second edition", + Hairer, E., + Nørsett, S., + Wanner, G., + Springer-Verlag, Berlin. 1996. + \endverbatim + +SourceFiles + rodas43.C + +\*---------------------------------------------------------------------------*/ + +#ifndef rodas43_H +#define rodas43_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class rodas43 Declaration +\*---------------------------------------------------------------------------*/ + +class rodas43 +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + mutable scalarField k1_; + mutable scalarField k2_; + mutable scalarField k3_; + mutable scalarField k4_; + mutable scalarField k5_; + mutable scalarField dy_; + mutable scalarField err_; + mutable scalarField dydx_; + mutable scalarField dfdx_; + mutable scalarSquareMatrix dfdy_; + mutable scalarSquareMatrix a_; + mutable labelList pivotIndices_; + + static const scalar + c2, c3, c4, + bet2p, bet3p, bet4p, + d1, d2, d3, d4, + a21, a31, a32, + a41, a42, a43, + a51, a52, a53, a54, + c21, c31, c32, + c41, c42, c43, + c51, c52, c53, c54, + c61, c62, c63, c64, c65, + gamma; + +public: + + //- Runtime type information + TypeName("rodas43"); + + + // Constructors + + //- Construct from ODE + rodas43(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/ODESystem/ODESystem.H b/src/ODE/ODESystem/ODESystem.H index 90da37c427..50b3fbffed 100644 --- a/src/ODE/ODESystem/ODESystem.H +++ b/src/ODE/ODESystem/ODESystem.H @@ -51,7 +51,7 @@ public: // Constructors - //- Construct null8 + //- Construct null ODESystem() {} diff --git a/src/OSspecific/POSIX/signals/sigFpe.C b/src/OSspecific/POSIX/signals/sigFpe.C index 0f74e43a69..19861c86df 100644 --- a/src/OSspecific/POSIX/signals/sigFpe.C +++ b/src/OSspecific/POSIX/signals/sigFpe.C @@ -45,6 +45,7 @@ License #endif #include +#include // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -53,26 +54,7 @@ struct sigaction Foam::sigFpe::oldAction_; void Foam::sigFpe::fillSignallingNan(UList& lst) { -#ifdef LINUX - - // initialize to signalling NaN -# ifdef WM_SP - - const uint32_t sNAN = 0x7ff7fffflu; - uint32_t* dPtr = reinterpret_cast(lst.begin()); - -# else - - const uint64_t sNAN = 0x7ff7ffffffffffffllu; - uint64_t* dPtr = reinterpret_cast(lst.begin()); - -# endif - - forAll(lst, i) - { - *dPtr++ = sNAN; - } -#endif + lst = std::numeric_limits::signaling_NaN(); } diff --git a/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C b/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C index 477d724e6b..42ac74b405 100644 --- a/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C +++ b/src/OpenFOAM/containers/Lists/FixedList/FixedListIO.C @@ -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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -198,7 +198,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const FixedList& L) // Write end delimiter os << token::END_BLOCK; } - else if (Size < 11 && contiguous()) + else if (Size <= 1 ||(Size < 11 && contiguous())) { // Write start delimiter os << token::BEGIN_LIST; diff --git a/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectListIO.C b/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectListIO.C index ec7efc1f1f..c79fd18150 100644 --- a/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectListIO.C +++ b/src/OpenFOAM/containers/Lists/UIndirectList/UIndirectListIO.C @@ -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) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -67,7 +67,7 @@ Foam::Ostream& Foam::operator<< // Write end delimiter os << token::END_BLOCK; } - else if (L.size() < 11 && contiguous()) + else if(L.size() <= 1 || (L.size() < 11 && contiguous())) { // Write size and start delimiter os << L.size() << token::BEGIN_LIST; diff --git a/src/OpenFOAM/containers/Lists/UList/UListIO.C b/src/OpenFOAM/containers/Lists/UList/UListIO.C index a8bd6bf1c3..38874ca849 100644 --- a/src/OpenFOAM/containers/Lists/UList/UListIO.C +++ b/src/OpenFOAM/containers/Lists/UList/UListIO.C @@ -92,7 +92,7 @@ Foam::Ostream& Foam::operator<<(Foam::Ostream& os, const Foam::UList& L) // Write end delimiter os << token::END_BLOCK; } - else if (L.size() == 1 || (L.size() < 11 && contiguous())) + else if (L.size() <= 1 || (L.size() < 11 && contiguous())) { // Write size and start delimiter os << L.size() << token::BEGIN_LIST; diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index d2f3fee22e..33c445f416 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -1328,24 +1328,9 @@ void Foam::polyMesh::findTetFacePt { const polyMesh& mesh = *this; - tetFaceI = -1; - tetPtI = -1; - - List cellTets = - polyMeshTetDecomposition::cellTetIndices(mesh, cellI); - - forAll(cellTets, tetI) - { - const tetIndices& cellTetIs = cellTets[tetI]; - - if (cellTetIs.tet(mesh).inside(pt)) - { - tetFaceI = cellTetIs.face(); - tetPtI = cellTetIs.tetPt(); - - return; - } - } + tetIndices tet(polyMeshTetDecomposition::findTet(mesh, cellI, pt)); + tetFaceI = tet.face(); + tetPtI = tet.tetPt(); } diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C index 7bd5bd0d4b..6c1d76d4ae 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C @@ -311,15 +311,14 @@ Foam::labelList Foam::polyMeshTetDecomposition::findFaceBasePts { FatalErrorIn ( - "Foam::labelList" - "Foam::polyMeshTetDecomposition::findFaceBasePts" + "labelList" + "polyMeshTetDecomposition::findFaceBasePts" "(" "const polyMesh& mesh, " "scalar tol, " "bool report" ")" - ) - << "Coupled face base point exchange failure for face " + ) << "Coupled face base point exchange failure for face " << fI << abort(FatalError); } @@ -561,8 +560,8 @@ Foam::List Foam::polyMeshTetDecomposition::faceTetIndices { WarningIn ( - "Foam::List " - "Foam::polyMeshTetDecomposition::faceTetIndices" + "List " + "polyMeshTetDecomposition::faceTetIndices" "(" "const polyMesh&, " "label, " @@ -613,6 +612,76 @@ Foam::List Foam::polyMeshTetDecomposition::faceTetIndices } +Foam::tetIndices Foam::polyMeshTetDecomposition::triangleTetIndices +( + const polyMesh& mesh, + const label fI, + const label cI, + const label tetPtI +) +{ + static label nWarnings = 0; + static const label maxWarnings = 100; + + const face& f = mesh.faces()[fI]; + bool own = (mesh.faceOwner()[fI] == cI); + label tetBasePtI = mesh.tetBasePtIs()[fI]; + if (tetBasePtI == -1) + { + if (nWarnings < maxWarnings) + { + WarningIn + ( + "tetIndices " + "polyMeshTetDecomposition::triangleTetIndices" + "(" + "const polyMesh&, " + "label, " + "label, " + "label" + ")" + ) << "No base point for face " << fI << ", " << f + << ", produces a valid tet decomposition." + << endl; + nWarnings++; + } + if (nWarnings == maxWarnings) + { + Warning<< "Suppressing any further warnings." << endl; + nWarnings++; + } + + tetBasePtI = 0; + } + + tetIndices faceTetIs; + + label facePtI = (tetPtI + tetBasePtI) % f.size(); + label otherFacePtI = f.fcIndex(facePtI); + + faceTetIs.cell() = cI; + + faceTetIs.face() = fI; + + faceTetIs.faceBasePt() = tetBasePtI; + + if (own) + { + faceTetIs.facePtA() = facePtI; + faceTetIs.facePtB() = otherFacePtI; + } + else + { + faceTetIs.facePtA() = otherFacePtI; + faceTetIs.facePtB() = facePtI; + } + + faceTetIs.tetPt() = tetPtI; + + return faceTetIs; +} + + Foam::List Foam::polyMeshTetDecomposition::cellTetIndices ( const polyMesh& mesh, @@ -644,4 +713,56 @@ Foam::List Foam::polyMeshTetDecomposition::cellTetIndices } +Foam::tetIndices Foam::polyMeshTetDecomposition::findTet +( + const polyMesh& mesh, + label cI, + const point& pt +) +{ + const faceList& pFaces = mesh.faces(); + const cellList& pCells = mesh.cells(); + + const cell& thisCell = pCells[cI]; + + tetIndices tetContainingPt; + + + forAll(thisCell, cFI) + { + label fI = thisCell[cFI]; + const face& f = pFaces[fI]; + + for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++) + { + // Get tetIndices of face triangle + tetIndices faceTetIs + ( + triangleTetIndices + ( + mesh, + fI, + cI, + tetPtI + ) + ); + + // Check if inside + if (faceTetIs.tet(mesh).inside(pt)) + { + tetContainingPt = faceTetIs; + break; + } + } + + if (tetContainingPt.cell() != -1) + { + break; + } + } + + return tetContainingPt; +} + + // ************************************************************************* // diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H index af1defbf65..37b56a6844 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H @@ -130,6 +130,15 @@ public: label cI ); + //- Return the tet decomposition of the given triangle of the given face + static tetIndices triangleTetIndices + ( + const polyMesh& mesh, + label fI, + label cI, + const label tetPtI // offset in face + ); + //- Return the tet decomposition of the given cell, see // findFacePt for the meaning of the indices static List cellTetIndices @@ -137,6 +146,14 @@ public: const polyMesh& mesh, label cI ); + + //- Find the tet decomposition of the cell containing the given point + static tetIndices findTet + ( + const polyMesh& mesh, + label cI, + const point& pt + ); }; diff --git a/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C b/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C index 4730621e08..970e0dac47 100644 --- a/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C +++ b/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,14 +47,15 @@ extrudePatchMesh::extrudePatchMesh ( const fvMesh& mesh, const fvPatch& patch, - const dictionary& dict + const dictionary& dict, + const word regionName ) : fvMesh ( IOobject ( - dict.lookup("region"), + regionName, mesh.facesInstance(), mesh, IOobject::READ_IF_PRESENT, diff --git a/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.H b/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.H index b960988e20..0ced75eec6 100644 --- a/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.H +++ b/src/dynamicMesh/extrudePatchMesh/extrudePatchMesh.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,12 +25,11 @@ Class extrudePatchMesh Description - Mesh at a patch created on the fly. The following entried should be used + Mesh at a patch created on the fly. The following entry should be used on the field boundary dictionary: // New Shell mesh data - region "regionMesh"; extrudeModel linearNormal; linearNormalCoeffs { @@ -120,7 +119,8 @@ public: ( const fvMesh&, const fvPatch&, - const dictionary& + const dictionary&, + const word ); diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C index 00aabb3bd6..27bea7ae47 100644 --- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C +++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C @@ -1784,7 +1784,7 @@ void Foam::faceCoupleInfo::subDivisionMatch << " points:" << masterPoints[masterEdge[0]] << ' ' << masterPoints[masterEdge[1]] << " corresponding cut edge: (" << cutPoint0 - << ' ' << cutPoint1 + << ") " << cutPoint1 << abort(FatalError); } diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C index 455769d7d2..c5d24780d3 100644 --- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C +++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C @@ -38,51 +38,6 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -// Append all mapped elements of a list to a DynamicList -void Foam::polyMeshAdder::append -( - const labelList& map, - const labelList& lst, - DynamicList