diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C index 7848ca0d8b..f3f1658a7d 100644 --- a/applications/test/ODE/Test-ODE.C +++ b/applications/test/ODE/Test-ODE.C @@ -109,8 +109,11 @@ int main(int argc, char *argv[]) // Create the ODE system testODE ode; + dictionary dict; + dict.add("solver", args[1]); + // Create the selected ODE system solver - autoPtr odeSolver = ODESolver::New(args[1], ode); + autoPtr odeSolver = ODESolver::New(ode, dict); // Initialise the ODE system fields scalar xStart = 1.0; @@ -124,26 +127,26 @@ int main(int argc, char *argv[]) scalarField dyStart(ode.nEqns()); ode.derivatives(xStart, yStart, dyStart); - Info<< setw(10) << "eps" << setw(12) << "hEst"; - Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl; + Info<< setw(10) << "relTol" << setw(12) << "dxEst"; + Info<< setw(13) << "dxDid" << setw(14) << "dxNext" << endl; Info<< setprecision(6); for (label i=0; i<15; i++) { - scalar eps = ::Foam::exp(-scalar(i + 1)); + scalar relTol = ::Foam::exp(-scalar(i + 1)); scalar x = xStart; scalarField y(yStart); - scalarField dydx(dyStart); - scalar hEst = 0.6; - scalar hDid, hNext; + scalar dxEst = 0.6; + scalar dxNext = dxEst; - odeSolver->solve(ode, x, y, dydx, eps, hEst, hDid, hNext); + odeSolver->relTol() = relTol; + odeSolver->solve(ode, x, y, dxNext); - Info<< scientific << setw(13) << eps; - Info<< fixed << setw(11) << hEst; - Info<< setw(13) << hDid << setw(13) << hNext + Info<< scientific << setw(13) << relTol; + Info<< fixed << setw(11) << dxEst; + Info<< setw(13) << x - xStart << setw(13) << dxNext << setw(13) << y[0] << setw(13) << y[1] << setw(13) << y[2] << setw(13) << y[3] << endl; @@ -159,12 +162,13 @@ int main(int argc, char *argv[]) yEnd[2] = ::Foam::jn(2, xEnd); yEnd[3] = ::Foam::jn(3, xEnd); - scalar hEst = 0.5; + scalar dxEst = 0.5; - odeSolver->solve(ode, x, xEnd, y, 1e-4, hEst); + odeSolver->relTol() = 1e-4; + odeSolver->solve(ode, x, xEnd, y, dxEst); Info<< nl << "Analytical: y(2.0) = " << yEnd << endl; - Info << "Numerical: y(2.0) = " << y << ", hEst = " << hEst << endl; + Info << "Numerical: y(2.0) = " << y << ", dxEst = " << dxEst << endl; Info<< "\nEnd\n" << endl; diff --git a/src/ODE/Make/files b/src/ODE/Make/files index cb569e78f5..a22ec1d3be 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -2,8 +2,15 @@ ODESolvers/ODESolver/ODESolver.C ODESolvers/ODESolver/ODESolverNew.C ODESolvers/adaptiveSolver/adaptiveSolver.C -ODESolvers/RKCK5/RKCK5.C -ODESolvers/KRR4/KRR4.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 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::EulerSI + +Description + 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 + EulerSI.C + +\*---------------------------------------------------------------------------*/ + +#ifndef EulerSI_H +#define EulerSI_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class EulerSI Declaration +\*---------------------------------------------------------------------------*/ + +class EulerSI +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + mutable scalarField err_; + mutable scalarField dydx_; + mutable scalarField dfdx_; + mutable scalarSquareMatrix dfdy_; + mutable scalarSquareMatrix a_; + mutable labelList pivotIndices_; + + +public: + + //- Runtime type information + TypeName("EulerSI"); + + + // Constructors + + //- Construct from 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, + scalar& dxTry + ) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 60a599e8d2..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,76 +30,103 @@ License namespace Foam { defineTypeNameAndDebug(ODESolver, 0); - defineRunTimeSelectionTable(ODESolver, ODE); + defineRunTimeSelectionTable(ODESolver, dictionary); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::ODESolver::ODESolver(const ODESystem& ode) +Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict) : n_(ode.nEqns()), - dydx_(n_) + absTol_(n_, dict.lookupOrDefault("absTol", SMALL)), + relTol_(n_, dict.lookupOrDefault("relTol", 1e-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.0) + // Check if this is a truncated step and set dxEst to integrate to xEnd + if ((x + dxEst - xEnd)*(x + dxEst - xStart) > 0) { - h = xEnd - x; - hPrev = hNext; + truncated = true; + dxEst = xEnd - x; } - hNext = 0; - scalar hDid; - solve(ode, x, y, dydx_, eps, h, hDid, hNext); + // Integrate as far as possible up to dxEst + solve(ode, x, y, dxEst); - if ((x - xEnd)*(xEnd - xStart) >= 0.0) + // Check if reached xEnd + if ((x - xEnd)*(xEnd - xStart) >= 0) { - if (hPrev != 0) + if (nStep > 0 && truncated) { - hEst = hPrev; - } - else - { - hEst = hNext; + dxEst = dxEst0; } return; } - - h = hNext; } FatalErrorIn ( "ODESolver::solve" "(const ODESystem& ode, const scalar xStart, const scalar xEnd," - "scalarField& yStart, const scalar eps, scalar& hEst) const" - ) << "Too many integration steps" + "scalarField& y, scalar& dxEst) const" + ) << "Integration steps greater than maximum " << maxSteps_ << exit(FatalError); } + // ************************************************************************* // diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H index a293088685..bd767a9063 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H @@ -53,13 +53,30 @@ class ODESolver protected: - // Private data + // Protected data + //- Size of the ODESystem label n_; - mutable scalarField dydx_; + + //- Absolute convergence tolerance per step + scalarField absTol_; + + //- Relative convergence tolerance per step + scalarField relTol_; + + //- The maximum number of sub-steps allowed for the integration step + label maxSteps_; - // Private Member Functions + // Protected Member Functions + + //- Return the nomalized scalar error + scalar normalizeError + ( + const scalarField& y0, + const scalarField& y, + const scalarField& err + ) const; //- Disallow default bitwise copy construct ODESolver(const ODESolver&); @@ -80,16 +97,24 @@ public: ( autoPtr, ODESolver, - ODE, - (const ODESystem& ode), - (ode) + dictionary, + (const ODESystem& ode, const dictionary& dict), + (ode, dict) ); // Constructors - //- Construct for given ODE - ODESolver(const ODESystem& ode); + //- Construct for given ODESystem + ODESolver(const ODESystem& ode, const dictionary& dict); + + //- Construct for given ODESystem specifying tolerances + ODESolver + ( + const ODESystem& ode, + const scalarField& absTol, + const scalarField& relTol + ); // Selectors @@ -97,8 +122,8 @@ public: //- Select null constructed static autoPtr New ( - const word& ODESolverTypeName, - const ODESystem& ode + const ODESystem& ode, + const dictionary& dict ); @@ -109,27 +134,37 @@ public: // Member Functions + scalarField& absTol() + { + return absTol_; + } + + scalarField& relTol() + { + return relTol_; + } + + //- Solve the ODE system as far as possible upto dxTry + // adjusting the step as necessary to provide a solution within + // the specified tolerance. + // Update the state and return an estimate for the next step in dxTry virtual void solve ( const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const = 0; - + //- Solve the ODE system from xStart to xEnd, update the state + // and return an estimate for the next step in dxTry virtual void solve ( const ODESystem& ode, const scalar xStart, const scalar xEnd, scalarField& y, - const scalar eps, - scalar& hEst + scalar& dxEst ) const; }; diff --git a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C index aa8ed538cc..ef0d6e3454 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolverNew.C +++ b/src/ODE/ODESolvers/ODESolver/ODESolverNew.C @@ -29,29 +29,30 @@ License Foam::autoPtr Foam::ODESolver::New ( - const Foam::word& ODESolverTypeName, - const Foam::ODESystem& ode + const ODESystem& odes, + const dictionary& dict ) { + word ODESolverTypeName(dict.lookup("solver")); Info<< "Selecting ODE solver " << ODESolverTypeName << endl; - ODEConstructorTable::iterator cstrIter = - ODEConstructorTablePtr_->find(ODESolverTypeName); + dictionaryConstructorTable::iterator cstrIter = + dictionaryConstructorTablePtr_->find(ODESolverTypeName); - if (cstrIter == ODEConstructorTablePtr_->end()) + if (cstrIter == dictionaryConstructorTablePtr_->end()) { FatalErrorIn ( "ODESolver::New" - "(const word& ODESolverTypeName, const ODESystem& ode)" + "(const dictionary& dict, const ODESystem& odes)" ) << "Unknown ODESolver type " << ODESolverTypeName << nl << nl << "Valid ODESolvers are : " << endl - << ODEConstructorTablePtr_->sortedToc() + << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError); } - return autoPtr(cstrIter()(ode)); + return autoPtr(cstrIter()(odes, dict)); } diff --git a/src/ODE/ODESolvers/RKCK5/RKCK5.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C similarity index 51% rename from src/ODE/ODESolvers/RKCK5/RKCK5.C rename to src/ODE/ODESolvers/RKCK45/RKCK45.C index 5050c0e7ca..2de3cdf99e 100644 --- a/src/ODE/ODESolvers/RKCK5/RKCK5.C +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -23,41 +23,58 @@ License \*---------------------------------------------------------------------------*/ -#include "RKCK5.H" +#include "RKCK45.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(RKCK5, 0); - addToRunTimeSelectionTable(ODESolver, RKCK5, ODE); + defineTypeNameAndDebug(RKCK45, 0); + addToRunTimeSelectionTable(ODESolver, RKCK45, dictionary); const scalar - RKCK5::a2 = 0.2, RKCK5::a3 = 0.3, RKCK5::a4 = 0.6, RKCK5::a5 = 1.0, - RKCK5::a6 = 0.875, - RKCK5::b21 = 0.2, RKCK5::b31 = 3.0/40.0, RKCK5::b32 = 9.0/40.0, - RKCK5::b41 = 0.3, RKCK5::b42 = -0.9, RKCK5::b43 = 1.2, - RKCK5::b51 = -11.0/54.0, RKCK5::b52 = 2.5, RKCK5::b53 = -70.0/27.0, - RKCK5::b54 = 35.0/27.0, - RKCK5::b61 = 1631.0/55296.0, RKCK5::b62 = 175.0/512.0, - RKCK5::b63 = 575.0/13824.0, RKCK5::b64 = 44275.0/110592.0, - RKCK5::b65 = 253.0/4096.0, - RKCK5::c1 = 37.0/378.0, RKCK5::c3 = 250.0/621.0, - RKCK5::c4 = 125.0/594.0, RKCK5::c6 = 512.0/1771.0, - RKCK5::dc1 = RKCK5::c1 - 2825.0/27648.0, - RKCK5::dc3 = RKCK5::c3 - 18575.0/48384.0, - RKCK5::dc4 = RKCK5::c4 - 13525.0/55296.0, RKCK5::dc5 = -277.00/14336.0, - RKCK5::dc6 = RKCK5::c6 - 0.25; + RKCK45::c2 = 1.0/5.0, + RKCK45::c3 = 3.0/10.0, + RKCK45::c4 = 3.0/5.0, + RKCK45::c5 = 1.0, + RKCK45::c6 = 7.0/8.0, + + RKCK45::a21 = 1.0/5.0, + RKCK45::a31 = 3.0/40.0, + RKCK45::a32 = 9.0/40.0, + RKCK45::a41 = 3.0/10.0, + RKCK45::a42 = -9.0/10.0, + RKCK45::a43 = 6.0/5.0, + RKCK45::a51 = -11.0/54.0, + RKCK45::a52 = 5.0/2.0, + RKCK45::a53 = -70.0/27.0, + RKCK45::a54 = 35.0/27.0, + RKCK45::a61 = 1631.0/55296.0, + RKCK45::a62 = 175.0/512.0, + RKCK45::a63 = 575.0/13824.0, + RKCK45::a64 = 44275.0/110592.0, + RKCK45::a65 = 253.0/4096.0, + + RKCK45::b1 = 37.0/378.0, + RKCK45::b3 = 250.0/621.0, + RKCK45::b4 = 125.0/594.0, + RKCK45::b6 = 512.0/1771.0, + + RKCK45::e1 = RKCK45::b1 - 2825.0/27648.0, + RKCK45::e3 = RKCK45::b3 - 18575.0/48384.0, + RKCK45::e4 = RKCK45::b4 - 13525.0/55296.0, + RKCK45::e5 = -277.00/14336.0, + RKCK45::e6 = RKCK45::b6 - 0.25; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::RKCK5::RKCK5(const ODESystem& ode) +Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict) : - ODESolver(ode), - adaptiveSolver(ode), + ODESolver(ode, dict), + adaptiveSolver(ode, dict), yTemp_(n_), k2_(n_), k3_(n_), @@ -70,85 +87,80 @@ Foam::RKCK5::RKCK5(const ODESystem& ode) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::RKCK5::solve +Foam::scalar Foam::RKCK45::solve ( const ODESystem& odes, const scalar x0, const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const { forAll(yTemp_, i) { - yTemp_[i] = y0[i] + b21*dx*dydx0[i]; + yTemp_[i] = y0[i] + a21*dx*dydx0[i]; } - odes.derivatives(x0 + a2*dx, yTemp_, k2_); + odes.derivatives(x0 + c2*dx, yTemp_, k2_); forAll(yTemp_, i) { - yTemp_[i] = y0[i] + dx*(b31*dydx0[i] + b32*k2_[i]); + yTemp_[i] = y0[i] + dx*(a31*dydx0[i] + a32*k2_[i]); } - odes.derivatives(x0 + a3*dx, yTemp_, k3_); + odes.derivatives(x0 + c3*dx, yTemp_, k3_); forAll(yTemp_, i) { - yTemp_[i] = y0[i] + dx*(b41*dydx0[i] + b42*k2_[i] + b43*k3_[i]); + yTemp_[i] = y0[i] + dx*(a41*dydx0[i] + a42*k2_[i] + a43*k3_[i]); } - odes.derivatives(x0 + a4*dx, yTemp_, k4_); + odes.derivatives(x0 + c4*dx, yTemp_, k4_); forAll(yTemp_, i) { yTemp_[i] = y0[i] - + dx*(b51*dydx0[i] + b52*k2_[i] + b53*k3_[i] + b54*k4_[i]); + + dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]); } - odes.derivatives(x0 + a5*dx, yTemp_, k5_); + odes.derivatives(x0 + c5*dx, yTemp_, k5_); forAll(yTemp_, i) { yTemp_[i] = y0[i] + dx - *(b61*dydx0[i] + b62*k2_[i] + b63*k3_[i] + b64*k4_[i] + b65*k5_[i]); + *(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]); } - odes.derivatives(x0 + a6*dx, yTemp_, k6_); + odes.derivatives(x0 + c6*dx, yTemp_, k6_); forAll(y, i) { y[i] = y0[i] - + dx*(c1*dydx0[i] + c3*k3_[i] + c4*k4_[i] + c6*k6_[i]); + + dx*(b1*dydx0[i] + b3*k3_[i] + b4*k4_[i] + b6*k6_[i]); } forAll(err_, i) { err_[i] = dx - *(dc1*dydx0[i] + dc3*k3_[i] + dc4*k4_[i] + dc5*k5_[i] + dc6*k6_[i]); + *(e1*dydx0[i] + e3*k3_[i] + e4*k4_[i] + e5*k5_[i] + e6*k6_[i]); } - return normalizeError(eps, y0, y, err_); + return normalizeError(y0, y, err_); } -void Foam::RKCK5::solve +void Foam::RKCK45::solve ( const ODESystem& odes, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const { - adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext); + adaptiveSolver::solve(odes, x, y, dxTry); } diff --git a/src/ODE/ODESolvers/RKCK5/RKCK5.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H similarity index 77% rename from src/ODE/ODESolvers/RKCK5/RKCK5.H rename to src/ODE/ODESolvers/RKCK45/RKCK45.H index 73a33e36bd..a6b551a475 100644 --- a/src/ODE/ODESolvers/RKCK5/RKCK5.H +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -22,10 +22,10 @@ License along with OpenFOAM. If not, see . Class - Foam::RKCK5 + Foam::RKCK45 Description - 5th Order Cash-Karp Runge-Kutta ODE solver + 4/5th Order Cash-Karp Runge-Kutta ODE solver. References: \verbatim @@ -36,13 +36,23 @@ Description 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 - RKCK5.C + RKCK45.C \*---------------------------------------------------------------------------*/ -#ifndef RKCK5_H -#define RKCK5_H +#ifndef RKCK45_H +#define RKCK45_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -53,10 +63,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class RKCK5 Declaration + Class RKCK45 Declaration \*---------------------------------------------------------------------------*/ -class RKCK5 +class RKCK45 : public ODESolver, public adaptiveSolver @@ -65,11 +75,11 @@ class RKCK5 //- RKCK Constants static const scalar - a2, a3, a4, a5, a6, - b21, b31, b32, b41, b42, b43, - b51, b52, b53, b54, b61, b62, b63, b64, b65, - c1, c3, c4, c6, - dc1, dc3, dc4, dc5, dc6; + c2, c3, c4, c5, c6, + a21, a31, a32, a41, a42, a43, a51, a52, a53, a54, + a61, a62, a63, a64, a65, + b1, b3, b4, b6, + e1, e3, e4, e5, e6; // Temporary fields mutable scalarField yTemp_; @@ -86,13 +96,13 @@ class RKCK5 public: //- Runtime type information - TypeName("RKCK5"); + TypeName("RKCK45"); // Constructors //- Construct from ODE - RKCK5(const ODESystem& ode); + RKCK45(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -105,8 +115,7 @@ public: const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const; //- Solve the ODE system and the update the state @@ -115,11 +124,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/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/KRR4/KRR4.C b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C similarity index 50% rename from src/ODE/ODESolvers/KRR4/KRR4.C rename to src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C index 6f0b938782..de6b968348 100644 --- a/src/ODE/ODESolvers/KRR4/KRR4.C +++ b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -23,42 +23,40 @@ License \*---------------------------------------------------------------------------*/ -#include "KRR4.H" +#include "Rosenbrock21.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(KRR4, 0); - addToRunTimeSelectionTable(ODESolver, KRR4, ODE); + defineTypeNameAndDebug(Rosenbrock21, 0); + addToRunTimeSelectionTable(ODESolver, Rosenbrock21, dictionary); 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; + 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::KRR4::KRR4(const ODESystem& ode) +Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict) : - ODESolver(ode), - adaptiveSolver(ode), - g1_(n_), - g2_(n_), - g3_(n_), - g4_(n_), + ODESolver(ode, dict), + adaptiveSolver(ode, dict), + k1_(n_), + k2_(n_), err_(n_), + dydx_(n_), dfdx_(n_), dfdy_(n_, n_), a_(n_, n_), @@ -68,15 +66,14 @@ Foam::KRR4::KRR4(const ODESystem& ode) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::KRR4::solve +Foam::scalar Foam::Rosenbrock21::solve ( const ODESystem& ode, const scalar x0, const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const { ode.jacobian(x0, y0, dfdx_, dfdy_); @@ -93,72 +90,49 @@ Foam::scalar Foam::KRR4::solve LUDecompose(a_, pivotIndices_); - forAll(g1_, i) + // Calculate k1: + forAll(k1_, i) { - g1_[i] = dydx0[i] + dx*c1X*dfdx_[i]; + k1_[i] = dydx0[i] + dx*d1*dfdx_[i]; } - LUBacksubstitute(a_, pivotIndices_, g1_); + LUBacksubstitute(a_, pivotIndices_, k1_); + // Calculate k2: forAll(y, i) { - y[i] = y0[i] + a21*g1_[i]; + y[i] = y0[i] + a21*k1_[i]; } - ode.derivatives(x0 + a2X*dx, y, dydx_); + ode.derivatives(x0 + c2*dx, y, dydx_); - forAll(g2_, i) + forAll(k2_, i) { - g2_[i] = dydx_[i] + dx*c2X*dfdx_[i] + c21*g1_[i]/dx; + k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx; } - LUBacksubstitute(a_, pivotIndices_, g2_); + LUBacksubstitute(a_, pivotIndices_, k2_); + // Calculate error and update state: forAll(y, i) { - y[i] = y0[i] + a31*g1_[i] + a32*g2_[i]; + y[i] = y0[i] + b1*k1_[i] + b2*k2_[i]; + err_[i] = e1*k1_[i] + e2*k2_[i]; } - ode.derivatives(x0 + a3X*dx, y, dydx_); - - forAll(g3_, i) - { - g3_[i] = dydx0[i] + dx*c3X*dfdx_[i] + (c31*g1_[i] + c32*g2_[i])/dx; - } - - LUBacksubstitute(a_, pivotIndices_, g3_); - - forAll(g4_, i) - { - g4_[i] = dydx_[i] + dx*c4X*dfdx_[i] - + (c41*g1_[i] + c42*g2_[i] + c43*g3_[i])/dx; - } - - LUBacksubstitute(a_, pivotIndices_, g4_); - - forAll(y, i) - { - y[i] = y0[i] + b1*g1_[i] + b2*g2_[i] + b3*g3_[i] + b4*g4_[i]; - err_[i] = e1*g1_[i] + e2*g2_[i] + e3*g3_[i] + e4*g4_[i]; - } - - return normalizeError(eps, y0, y, err_); + return normalizeError(y0, y, err_); } -void Foam::KRR4::solve +void Foam::Rosenbrock21::solve ( const ODESystem& odes, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const { - adaptiveSolver::solve(odes, x, y, dydx, eps, dxTry, dxDid, dxNext); + adaptiveSolver::solve(odes, x, y, dxTry); } diff --git a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H new file mode 100644 index 0000000000..8200144751 --- /dev/null +++ b/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.H @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::Rosenbrock21 + +Description + Embedded Rosenbrock ODE solver of order (1)2. + +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..adb607b689 --- /dev/null +++ b/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C @@ -0,0 +1,186 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + + Rosenbrock43::a21 = 2.0, + Rosenbrock43::a31 = 48.0/25.0, + Rosenbrock43::a32 = 6.0/25.0, + + Rosenbrock43::c21 = -8.0, + 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.0, + Rosenbrock43::e4 = 125.0/108.0, + + Rosenbrock43::gamma = 0.5, + Rosenbrock43::c2 = 1.0, + 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::KRR4 + Foam::Rosenbrock43 Description - 4th Order Kaps Rentrop Rosenbrock stiff ODE solver + Embedded Rosenbrock ODE solver of order (3)4. - References: + Based on code from: \verbatim - "Generalized Runge-Kutta methods of order four with stepsize control - for stiff ordinary differential equations" - Kaps, P., - Rentrop, P., - Numerische Thematic, vol. 33, 1979, pp. 55–68. + "Solving Ordinary Differential Equations II: Stiff + and Differential-Algebraic Problems, second edition", + Hairer, E., + Nørsett, S., + Wanner, G., + Springer-Verlag, Berlin. 1996. \endverbatim Constants from: @@ -44,12 +45,12 @@ Description \endverbatim SourceFiles - KRR4.C + Rosenbrock43.C \*---------------------------------------------------------------------------*/ -#ifndef KRR4_H -#define KRR4_H +#ifndef Rosenbrock43_H +#define Rosenbrock43_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -60,46 +61,48 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class KRR4 Declaration + Class Rosenbrock43 Declaration \*---------------------------------------------------------------------------*/ -class KRR4 +class Rosenbrock43 : public ODESolver, public adaptiveSolver { // Private data - mutable scalarField g1_; - mutable scalarField g2_; - mutable scalarField g3_; - mutable scalarField g4_; + 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 - gamma, a21, a31, a32, - c21, c31, c32, c41, c42, c43, + c21, c31, c32, + c41, c42, c43, b1, b2, b3, b4, e1, e2, e3, e4, - c1X, c2X, c3X, c4X, - a2X, a3X; + gamma, + c2, c3, + d1, d2, d3, d4; public: //- Runtime type information - TypeName("KRR4"); + TypeName("Rosenbrock43"); // Constructors //- Construct from ODE - KRR4(const ODESystem& ode); + Rosenbrock43(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -112,8 +115,7 @@ public: const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const; //- Solve the ODE system and the update the state @@ -122,11 +124,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C index 964e288949..28d4706792 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.C +++ b/src/ODE/ODESolvers/SIBS/SIBS.C @@ -31,8 +31,7 @@ License namespace Foam { defineTypeNameAndDebug(SIBS, 0); - - addToRunTimeSelectionTable(ODESolver, SIBS, ODE); + addToRunTimeSelectionTable(ODESolver, SIBS, dictionary); const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70}; @@ -47,9 +46,9 @@ namespace Foam // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::SIBS::SIBS(const ODESystem& ode) +Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict) : - ODESolver(ode), + ODESolver(ode, dict), a_(iMaxX_, 0.0), alpha_(kMaxX_, kMaxX_, 0.0), d_p_(n_, kMaxX_, 0.0), @@ -59,6 +58,7 @@ Foam::SIBS::SIBS(const ODESystem& ode) yTemp_(n_, 0.0), ySeq_(n_, 0.0), yErr_(n_, 0.0), + dydx0_(n_), dfdx_(n_, 0.0), dfdy_(n_, n_, 0.0), first_(1), @@ -73,19 +73,18 @@ void Foam::SIBS::solve const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const { + ode.derivatives(x, y, dydx0_); + + scalar h = dxTry; bool exitflag = false; - if (eps != epsOld_) + if (relTol_[0] != epsOld_) { - hNext = xNew_ = -GREAT; - scalar eps1 = safe1*eps; + dxTry = xNew_ = -GREAT; + scalar eps1 = safe1*relTol_[0]; a_[0] = nSeq_[0] + 1; for (register label k=0; k= k && kOpt_ != kMax_ && !reduct) { scalar fact = max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX); if (a_[kOpt_ + 1]*fact <= wrkmin) { - hNext = h/fact; + dxTry = h/fact; kOpt_++; } } diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H index d77e3e128d..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,11 +126,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar hTry, - scalar& hDid, - scalar& hNext + scalar& dxTry ) const; }; diff --git a/src/ODE/ODESolvers/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/Trapezoid/Trapezoid.H b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H new file mode 100644 index 0000000000..327dbe44bc --- /dev/null +++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H @@ -0,0 +1,104 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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::Trapezoid + +Description + Trapezoidal ODE solver of order (1)2. + +SourceFiles + Trapezoid.C + +\*---------------------------------------------------------------------------*/ + +#ifndef Trapezoid_H +#define Trapezoid_H + +#include "ODESolver.H" +#include "adaptiveSolver.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class Trapezoid Declaration +\*---------------------------------------------------------------------------*/ + +class Trapezoid +: + public ODESolver, + public adaptiveSolver +{ + // Private data + + mutable scalarField err_; + + +public: + + //- Runtime type information + TypeName("Trapezoid"); + + + // Constructors + + //- Construct from ODE + Trapezoid(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/adaptiveSolver/adaptiveSolver.C b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C index 185b195158..13951ae7b3 100644 --- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C @@ -26,79 +26,50 @@ License #include "adaptiveSolver.H" #include "addToRunTimeSelectionTable.H" -// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // - -namespace Foam -{ - -const scalar - adaptiveSolver::safeScale=0.9, - adaptiveSolver::alphaInc=0.2, - adaptiveSolver::alphaDec=0.25, - adaptiveSolver::minScale=0.2, - adaptiveSolver::maxScale=10; -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::adaptiveSolver::adaptiveSolver(const ODESystem& ode) +Foam::adaptiveSolver::adaptiveSolver +( + const ODESystem& ode, + const dictionary& dict +) : + safeScale_(dict.lookupOrDefault("safeScale", 0.9)), + alphaInc_(dict.lookupOrDefault("alphaIncrease", 0.2)), + alphaDec_(dict.lookupOrDefault("alphaDecrease", 0.25)), + minScale_(dict.lookupOrDefault("minScale", 0.2)), + maxScale_(dict.lookupOrDefault("maxScale", 10)), + dydx0_(ode.nEqns()), yTemp_(ode.nEqns()) {} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::adaptiveSolver::normalizeError -( - const scalar eps, - const scalarField& y0, - const scalarField& y, - const scalarField& err -) const -{ - scalar epsAbs = 1e-10; - scalar epsRel = eps; - - // Calculate the maximum error - scalar maxErr = 0.0; - forAll(err, i) - { - scalar eps = epsAbs + epsRel*max(mag(y0[i]), mag(y[i])); - maxErr = max(maxErr, mag(err[i])/eps); - } - - return maxErr; -} - - void Foam::adaptiveSolver::solve ( const ODESystem& odes, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const { scalar dx = dxTry; scalar err = 0.0; + odes.derivatives(x, y, dydx0_); + // Loop over solver and adjust step-size as necessary // to achieve desired error do { // Solve step and provide error estimate - err = solve(odes, x, y, dydx, dx, yTemp_, eps); + err = solve(odes, x, y, dydx0_, dx, yTemp_); // If error is large reduce dx if (err > 1) { - scalar scale = max(safeScale*pow(err, -alphaDec), minScale); + scalar scale = max(safeScale_*pow(err, -alphaDec_), minScale_); dx *= scale; if (dx < VSMALL) @@ -111,12 +82,19 @@ void Foam::adaptiveSolver::solve } while (err > 1); // Update the state - dxDid = dx; x += dx; y = yTemp_; // If the error is small increase the step-size - dxNext = min(max(safeScale*pow(err, -alphaInc), minScale), maxScale)*dx; + 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 index 40f07c2d93..d86bb88ce1 100644 --- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -25,16 +25,6 @@ Class Foam::adaptiveSolver Description - 5th Order Cash-Karp Runge-Kutta ODE solver - - References: - \verbatim - "A variable order Runge-Kutta method for initial value problems - with rapidly varying right-hand sides" - Cash, J.R., - Karp, A.H. - ACM Transactions on Mathematical Software, vol. 16, 1990, pp. 201–222. - \endverbatim SourceFiles adaptiveSolver.C @@ -60,8 +50,12 @@ class adaptiveSolver // Private data //- Step-size adjustment controls - static const scalar safeScale, alphaInc, alphaDec, minScale, maxScale; + scalar safeScale_, alphaInc_, alphaDec_, minScale_, maxScale_; + //- Cache for dydx at the initial time + mutable scalarField dydx0_; + + //- Temprorary for the test-step solution mutable scalarField yTemp_; @@ -70,7 +64,7 @@ public: // Constructors //- Construct from ODESystem - adaptiveSolver(const ODESystem& ode); + adaptiveSolver(const ODESystem& ode, const dictionary& dict); //- Destructor @@ -80,15 +74,6 @@ public: // Member Functions - //- Return the nomalized scalar error - scalar normalizeError - ( - const scalar eps, - const scalarField& y0, - const scalarField& y, - const scalarField& err - ) const; - //- Solve a single step dx and return the error virtual scalar solve ( @@ -97,8 +82,7 @@ public: const scalarField& y0, const scalarField& dydx0, const scalar dx, - scalarField& y, - const scalar eps + scalarField& y ) const = 0; //- Solve the ODE system and the update the state @@ -107,11 +91,7 @@ public: const ODESystem& ode, scalar& x, scalarField& y, - scalarField& dydx, - const scalar eps, - const scalar dxTry, - scalar& dxDid, - scalar& dxNext + scalar& dxTry ) const; }; diff --git a/src/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/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/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