From d3368fcb1eec2ba222f8a96a4ca45413f86f61c6 Mon Sep 17 00:00:00 2001 From: Henry Date: Mon, 4 Nov 2013 18:22:49 +0000 Subject: [PATCH] Added Euler and Trapezoid ODE solvers --- src/ODE/Make/files | 2 + src/ODE/ODESolvers/Euler/Euler.C | 88 +++++++++++++++++ src/ODE/ODESolvers/Euler/Euler.H | 114 +++++++++++++++++++++++ src/ODE/ODESolvers/Trapezoid/Trapezoid.C | 93 ++++++++++++++++++ src/ODE/ODESolvers/Trapezoid/Trapezoid.H | 104 +++++++++++++++++++++ 5 files changed, 401 insertions(+) create mode 100644 src/ODE/ODESolvers/Euler/Euler.C create mode 100644 src/ODE/ODESolvers/Euler/Euler.H create mode 100644 src/ODE/ODESolvers/Trapezoid/Trapezoid.C create mode 100644 src/ODE/ODESolvers/Trapezoid/Trapezoid.H diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 765fb7b207..98bbffc29c 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -2,7 +2,9 @@ ODESolvers/ODESolver/ODESolver.C ODESolvers/ODESolver/ODESolverNew.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 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/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 + +// ************************************************************************* //