From d50bce000d988d512f1f82b26fa9ebe18026e6ba Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Mon, 11 Jul 2016 17:27:04 +0100 Subject: [PATCH] ODESolvers: Add support for efficient ODE solver resizing Note: this reuses the existing storage rather than costly reallocation which requires the initial allocation to be sufficient for the largest size the ODE system might have. Attempt to set a size larger than the initial size is a fatal error. --- src/ODE/ODESolvers/Euler/Euler.C | 19 +++- src/ODE/ODESolvers/Euler/Euler.H | 16 ++- src/ODE/ODESolvers/EulerSI/EulerSI.C | 22 ++++ src/ODE/ODESolvers/EulerSI/EulerSI.H | 16 ++- src/ODE/ODESolvers/ODESolver/ODESolver.C | 101 +++++++++++++----- src/ODE/ODESolvers/ODESolver/ODESolver.H | 48 ++++++--- src/ODE/ODESolvers/ODESolver/ODESolverI.H | 67 ++++++++++++ src/ODE/ODESolvers/RKCK45/RKCK45.C | 25 ++++- src/ODE/ODESolvers/RKCK45/RKCK45.H | 16 ++- src/ODE/ODESolvers/RKDP45/RKDP45.C | 25 ++++- src/ODE/ODESolvers/RKDP45/RKDP45.H | 14 ++- src/ODE/ODESolvers/RKF45/RKF45.C | 25 ++++- src/ODE/ODESolvers/RKF45/RKF45.H | 16 ++- .../ODESolvers/Rosenbrock12/Rosenbrock12.C | 24 +++++ .../ODESolvers/Rosenbrock12/Rosenbrock12.H | 16 ++- .../ODESolvers/Rosenbrock23/Rosenbrock23.C | 25 +++++ .../ODESolvers/Rosenbrock23/Rosenbrock23.H | 16 ++- .../ODESolvers/Rosenbrock34/Rosenbrock34.C | 26 +++++ .../ODESolvers/Rosenbrock34/Rosenbrock34.H | 16 ++- src/ODE/ODESolvers/SIBS/SIBS.C | 20 ++++ src/ODE/ODESolvers/SIBS/SIBS.H | 14 ++- src/ODE/ODESolvers/Trapezoid/Trapezoid.C | 19 +++- src/ODE/ODESolvers/Trapezoid/Trapezoid.H | 16 ++- .../adaptiveSolver/adaptiveSolver.C | 11 +- .../adaptiveSolver/adaptiveSolver.H | 7 +- src/ODE/ODESolvers/rodas23/rodas23.C | 26 +++++ src/ODE/ODESolvers/rodas23/rodas23.H | 16 ++- src/ODE/ODESolvers/rodas34/rodas34.C | 28 +++++ src/ODE/ODESolvers/rodas34/rodas34.H | 16 ++- src/ODE/ODESolvers/seulex/seulex.C | 24 +++++ src/ODE/ODESolvers/seulex/seulex.H | 14 ++- src/OpenFOAM/matrices/Matrix/Matrix.H | 3 + src/OpenFOAM/matrices/Matrix/MatrixI.H | 8 ++ .../matrices/SquareMatrix/SquareMatrix.H | 3 + .../matrices/SquareMatrix/SquareMatrixI.H | 7 ++ .../chemistryModel/chemistrySolver/ode/ode.C | 11 +- .../chemistryModel/chemistrySolver/ode/ode.H | 3 +- 37 files changed, 673 insertions(+), 106 deletions(-) create mode 100644 src/ODE/ODESolvers/ODESolver/ODESolverI.H diff --git a/src/ODE/ODESolvers/Euler/Euler.C b/src/ODE/ODESolvers/Euler/Euler.C index 1496099ae8..fccecee001 100644 --- a/src/ODE/ODESolvers/Euler/Euler.C +++ b/src/ODE/ODESolvers/Euler/Euler.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,6 +47,23 @@ Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::Euler::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(err_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::Euler::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/Euler/Euler.H b/src/ODE/ODESolvers/Euler/Euler.H index fecfc57fa0..c1bf6b128d 100644 --- a/src/ODE/ODESolvers/Euler/Euler.H +++ b/src/ODE/ODESolvers/Euler/Euler.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -75,17 +75,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem Euler(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~Euler() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -95,7 +103,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.C b/src/ODE/ODESolvers/EulerSI/EulerSI.C index 0648b26649..143187d5bf 100644 --- a/src/ODE/ODESolvers/EulerSI/EulerSI.C +++ b/src/ODE/ODESolvers/EulerSI/EulerSI.C @@ -52,6 +52,28 @@ Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::EulerSI::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(err_); + resizeField(dydx_); + resizeField(dfdx_); + resizeMatrix(dfdy_); + resizeMatrix(a_); + resizeField(pivotIndices_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::EulerSI::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.H b/src/ODE/ODESolvers/EulerSI/EulerSI.H index a6703cc20c..4ad1693fe8 100644 --- a/src/ODE/ODESolvers/EulerSI/EulerSI.H +++ b/src/ODE/ODESolvers/EulerSI/EulerSI.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -82,17 +82,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem EulerSI(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~EulerSI() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -102,7 +110,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.C b/src/ODE/ODESolvers/ODESolver/ODESolver.C index 823f159097..3dce5378e3 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-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -34,34 +34,7 @@ namespace Foam } -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // - -Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict) -: - odes_(ode), - n_(ode.nEqns()), - 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 -) -: - odes_(ode), - n_(ode.nEqns()), - absTol_(absTol), - relTol_(relTol), - maxSteps_(10000) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // Foam::scalar Foam::ODESolver::normalizeError ( @@ -82,6 +55,76 @@ Foam::scalar Foam::ODESolver::normalizeError } +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict) +: + odes_(ode), + maxN_(ode.nEqns()), + n_(ode.nEqns()), + 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 +) +: + odes_(ode), + maxN_(ode.nEqns()), + n_(ode.nEqns()), + absTol_(absTol), + relTol_(relTol), + maxSteps_(10000) +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::ODESolver::resize() +{ + if (odes_.nEqns() != n_) + { + if (odes_.nEqns() > maxN_) + { + FatalErrorInFunction + << "Specified number of equations " << odes_.nEqns() + << " greater than maximum " << maxN_ + << abort(FatalError); + } + + n_ = odes_.nEqns(); + + resizeField(absTol_); + resizeField(relTol_); + + return true; + } + else + { + return false; + } +} + + +void Foam::ODESolver::solve +( + scalar& x, + scalarField& y, + scalar& dxTry +) const +{ + stepState step(dxTry); + solve(x, y, step); + dxTry = step.dxTry; +} + + void Foam::ODESolver::solve ( scalar& x, diff --git a/src/ODE/ODESolvers/ODESolver/ODESolver.H b/src/ODE/ODESolvers/ODESolver/ODESolver.H index de4ef399f8..fe443187ad 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.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) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -58,8 +58,11 @@ protected: //- Reference to ODESystem const ODESystem& odes_; - //- Size of the ODESystem - label n_; + //- Maximum size of the ODESystem + const label maxN_; + + //- Size of the ODESystem (adjustable) + mutable label n_; //- Absolute convergence tolerance per step scalarField absTol_; @@ -90,6 +93,8 @@ protected: public: + friend class ODESystem; + //- Runtime type information TypeName("ODESolver"); @@ -161,15 +166,25 @@ public: // Member Functions - scalarField& absTol() - { - return absTol_; - } + //- Return the number of equations to solve + inline label nEqns() const; - scalarField& relTol() - { - return relTol_; - } + //- Return access to the absolute tolerance field + inline scalarField& absTol(); + + //- Return access to the relative tolerance field + inline scalarField& relTol(); + + //- Resize the ODE solver + virtual bool resize() = 0; + + template + static inline void resizeField(UList& f, const label n); + + template + inline void resizeField(UList& f) const; + + inline void resizeMatrix(scalarSquareMatrix& m) const; //- Solve the ODE system as far as possible upto dxTry // adjusting the step as necessary to provide a solution within @@ -180,12 +195,7 @@ public: scalar& x, scalarField& y, scalar& dxTry - ) const //= 0; - { - stepState step(dxTry); - solve(x, y, step); - dxTry = step.dxTry; - } + ) const; //- Solve the ODE system as far as possible upto dxTry // adjusting the step as necessary to provide a solution within @@ -216,6 +226,10 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#include "ODESolverI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/ODE/ODESolvers/ODESolver/ODESolverI.H b/src/ODE/ODESolvers/ODESolver/ODESolverI.H new file mode 100644 index 0000000000..e89f7bffaf --- /dev/null +++ b/src/ODE/ODESolvers/ODESolver/ODESolverI.H @@ -0,0 +1,67 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + + +inline Foam::label Foam::ODESolver::nEqns() const +{ + return n_; +} + + +inline Foam::scalarField& Foam::ODESolver::absTol() +{ + return absTol_; +} + + +inline Foam::scalarField& Foam::ODESolver::relTol() +{ + return relTol_; +} + + +template +inline void Foam::ODESolver::resizeField(UList& f, const label n) +{ + f.shallowCopy(UList(f.begin(), n)); +} + + +template +inline void Foam::ODESolver::resizeField(UList& f) const +{ + resizeField(f, n_); +} + + +inline void Foam::ODESolver::resizeMatrix(scalarSquareMatrix& m) const +{ + m.shallowResize(n_); +} + + +// ************************************************************************* // diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C index 30ae91de1a..2e8b2685d8 100644 --- a/src/ODE/ODESolvers/RKCK45/RKCK45.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) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -87,6 +87,29 @@ Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::RKCK45::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(yTemp_); + resizeField(k2_); + resizeField(k3_); + resizeField(k4_); + resizeField(k5_); + resizeField(k6_); + resizeField(err_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::RKCK45::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H index c376a5f423..356552761e 100644 --- a/src/ODE/ODESolvers/RKCK45/RKCK45.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) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -98,17 +98,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem RKCK45(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~RKCK45() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -118,7 +126,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.C b/src/ODE/ODESolvers/RKDP45/RKDP45.C index f9a3a78836..479aecc0e9 100644 --- a/src/ODE/ODESolvers/RKDP45/RKDP45.C +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -92,6 +92,29 @@ Foam::RKDP45::RKDP45(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::RKDP45::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(yTemp_); + resizeField(k2_); + resizeField(k3_); + resizeField(k4_); + resizeField(k5_); + resizeField(k6_); + resizeField(err_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::RKDP45::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.H b/src/ODE/ODESolvers/RKDP45/RKDP45.H index 9998c15f6b..75aa756f9c 100644 --- a/src/ODE/ODESolvers/RKDP45/RKDP45.H +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.H @@ -102,17 +102,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem RKDP45(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~RKDP45() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -122,7 +130,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/RKF45/RKF45.C b/src/ODE/ODESolvers/RKF45/RKF45.C index 92804e0634..7ca9a232ea 100644 --- a/src/ODE/ODESolvers/RKF45/RKF45.C +++ b/src/ODE/ODESolvers/RKF45/RKF45.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -88,6 +88,29 @@ Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::RKF45::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(yTemp_); + resizeField(k2_); + resizeField(k3_); + resizeField(k4_); + resizeField(k5_); + resizeField(k6_); + resizeField(err_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::RKF45::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/RKF45/RKF45.H b/src/ODE/ODESolvers/RKF45/RKF45.H index 5198d5553f..cc967715e1 100644 --- a/src/ODE/ODESolvers/RKF45/RKF45.H +++ b/src/ODE/ODESolvers/RKF45/RKF45.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -102,17 +102,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem RKF45(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~RKF45() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -122,7 +130,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C index fa163aec0e..5bb3cf7513 100644 --- a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C +++ b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C @@ -66,6 +66,30 @@ Foam::Rosenbrock12::Rosenbrock12(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::Rosenbrock12::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(k1_); + resizeField(k2_); + resizeField(err_); + resizeField(dydx_); + resizeField(dfdx_); + resizeMatrix(dfdy_); + resizeMatrix(a_); + resizeField(pivotIndices_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::Rosenbrock12::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H index b0bf1c79ef..44d86a8505 100644 --- a/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H +++ b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -92,17 +92,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem Rosenbrock12(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~Rosenbrock12() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -112,7 +120,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C index 3b6c8cf463..f664726e83 100644 --- a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C +++ b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C @@ -79,6 +79,31 @@ Foam::Rosenbrock23::Rosenbrock23(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::Rosenbrock23::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(k1_); + resizeField(k2_); + resizeField(k3_); + resizeField(err_); + resizeField(dydx_); + resizeField(dfdx_); + resizeMatrix(dfdy_); + resizeMatrix(a_); + resizeField(pivotIndices_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::Rosenbrock23::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H index 38541a7c99..571a9e7ff8 100644 --- a/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H +++ b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -96,17 +96,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem Rosenbrock23(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~Rosenbrock23() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -116,7 +124,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C index 5b761b267f..477de28616 100644 --- a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C +++ b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C @@ -124,6 +124,32 @@ Foam::Rosenbrock34::Rosenbrock34(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::Rosenbrock34::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(k1_); + resizeField(k2_); + resizeField(k3_); + resizeField(k4_); + resizeField(err_); + resizeField(dydx_); + resizeField(dfdx_); + resizeMatrix(dfdy_); + resizeMatrix(a_); + resizeField(pivotIndices_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::Rosenbrock34::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H index 2af18e5d40..f72db63e7c 100644 --- a/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H +++ b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -104,17 +104,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem Rosenbrock34(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~Rosenbrock34() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -124,7 +132,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C index a095a7c507..2bd63ece59 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.C +++ b/src/ODE/ODESolvers/SIBS/SIBS.C @@ -68,6 +68,26 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +bool Foam::SIBS::resize() +{ + if (ODESolver::resize()) + { + resizeField(yTemp_); + resizeField(ySeq_); + resizeField(yErr_); + resizeField(dydx0_); + resizeField(dfdx_); + resizeMatrix(dfdy_); + + return true; + } + else + { + return false; + } +} + + void Foam::SIBS::solve ( scalar& x, diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H index 0ef0a1d40c..9f03753019 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-2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -114,13 +114,21 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODE system SIBS(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~SIBS() + {} + + // Member Functions - void solve + //- Resize the ODE solver + virtual bool resize(); + + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/Trapezoid/Trapezoid.C b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C index 027308434d..395ffc92d1 100644 --- a/src/ODE/ODESolvers/Trapezoid/Trapezoid.C +++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -47,6 +47,23 @@ Foam::Trapezoid::Trapezoid(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::Trapezoid::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(err_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::Trapezoid::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/Trapezoid/Trapezoid.H b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H index 667b6ee2ad..5dce76d5a8 100644 --- a/src/ODE/ODESolvers/Trapezoid/Trapezoid.H +++ b/src/ODE/ODESolvers/Trapezoid/Trapezoid.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -65,17 +65,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem Trapezoid(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~Trapezoid() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -85,7 +93,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C index f868e921d0..330684261c 100644 --- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -46,6 +46,15 @@ Foam::adaptiveSolver::adaptiveSolver // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::adaptiveSolver::resize(const label n) +{ + ODESolver::resizeField(dydx0_, n); + ODESolver::resizeField(yTemp_, n); + + return true; +} + + void Foam::adaptiveSolver::solve ( const ODESystem& odes, diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H index 0d7dcb035c..eef9035731 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) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -74,6 +74,9 @@ public: // Member Functions + //- Resize the ODE solver + virtual bool resize(const label n); + //- Solve a single step dx and return the error virtual scalar solve ( @@ -85,7 +88,7 @@ public: ) const = 0; //- Solve the ODE system and the update the state - void solve + virtual void solve ( const ODESystem& ode, scalar& x, diff --git a/src/ODE/ODESolvers/rodas23/rodas23.C b/src/ODE/ODESolvers/rodas23/rodas23.C index 31d5e8acee..b11b508e0e 100644 --- a/src/ODE/ODESolvers/rodas23/rodas23.C +++ b/src/ODE/ODESolvers/rodas23/rodas23.C @@ -70,6 +70,32 @@ Foam::rodas23::rodas23(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::rodas23::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(k1_); + resizeField(k2_); + resizeField(k3_); + resizeField(dy_); + resizeField(err_); + resizeField(dydx_); + resizeField(dfdx_); + resizeMatrix(dfdy_); + resizeMatrix(a_); + resizeField(pivotIndices_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::rodas23::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/rodas23/rodas23.H b/src/ODE/ODESolvers/rodas23/rodas23.H index 75f0262813..42839bfa8f 100644 --- a/src/ODE/ODESolvers/rodas23/rodas23.H +++ b/src/ODE/ODESolvers/rodas23/rodas23.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -96,17 +96,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem rodas23(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~rodas23() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -116,7 +124,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/rodas34/rodas34.C b/src/ODE/ODESolvers/rodas34/rodas34.C index 5f403238d7..2898b51c2d 100644 --- a/src/ODE/ODESolvers/rodas34/rodas34.C +++ b/src/ODE/ODESolvers/rodas34/rodas34.C @@ -93,6 +93,34 @@ Foam::rodas34::rodas34(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +bool Foam::rodas34::resize() +{ + if (ODESolver::resize()) + { + adaptiveSolver::resize(n_); + + resizeField(k1_); + resizeField(k2_); + resizeField(k3_); + resizeField(k4_); + resizeField(k5_); + resizeField(dy_); + resizeField(err_); + resizeField(dydx_); + resizeField(dfdx_); + resizeMatrix(dfdy_); + resizeMatrix(a_); + resizeField(pivotIndices_); + + return true; + } + else + { + return false; + } +} + + Foam::scalar Foam::rodas34::solve ( const scalar x0, diff --git a/src/ODE/ODESolvers/rodas34/rodas34.H b/src/ODE/ODESolvers/rodas34/rodas34.H index 35cb468f6a..9098a1e8b0 100644 --- a/src/ODE/ODESolvers/rodas34/rodas34.H +++ b/src/ODE/ODESolvers/rodas34/rodas34.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -97,17 +97,25 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem rodas34(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~rodas34() + {} + + // Member Functions //- Inherit solve from ODESolver using ODESolver::solve; + //- Resize the ODE solver + virtual bool resize(); + //- Solve a single step dx and return the error - scalar solve + virtual scalar solve ( const scalar x0, const scalarField& y0, @@ -117,7 +125,7 @@ public: ) const; //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C index b07a9f4988..df163dd22d 100644 --- a/src/ODE/ODESolvers/seulex/seulex.C +++ b/src/ODE/ODESolvers/seulex/seulex.C @@ -211,6 +211,30 @@ void Foam::seulex::extrapolate } +bool Foam::seulex::resize() +{ + if (ODESolver::resize()) + { + resizeField(dfdx_); + resizeMatrix(dfdy_); + resizeMatrix(a_); + resizeField(pivotIndices_); + resizeField(y0_); + resizeField(ySequence_); + resizeField(scale_); + resizeField(dy_); + resizeField(yTemp_); + resizeField(dydx_); + + return true; + } + else + { + return false; + } +} + + void Foam::seulex::solve ( scalar& x, diff --git a/src/ODE/ODESolvers/seulex/seulex.H b/src/ODE/ODESolvers/seulex/seulex.H index c3cd540198..8739618ee7 100644 --- a/src/ODE/ODESolvers/seulex/seulex.H +++ b/src/ODE/ODESolvers/seulex/seulex.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -133,14 +133,22 @@ public: // Constructors - //- Construct from ODE + //- Construct from ODESystem seulex(const ODESystem& ode, const dictionary& dict); + //- Destructor + virtual ~seulex() + {} + + // Member Functions + //- Resize the ODE solver + virtual bool resize(); + //- Solve the ODE system and the update the state - void solve + virtual void solve ( scalar& x, scalarField& y, diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H index 22e289a689..884979f3d8 100644 --- a/src/OpenFOAM/matrices/Matrix/Matrix.H +++ b/src/OpenFOAM/matrices/Matrix/Matrix.H @@ -249,6 +249,9 @@ public: //- Resize the matrix preserving the elements void setSize(const label m, const label n); + //- Resize the matrix without reallocating storage (unsafe) + inline void shallowResize(const label m, const label n); + //- Return the transpose of the matrix Form T() const; diff --git a/src/OpenFOAM/matrices/Matrix/MatrixI.H b/src/OpenFOAM/matrices/Matrix/MatrixI.H index 73992aac77..2ba29e355d 100644 --- a/src/OpenFOAM/matrices/Matrix/MatrixI.H +++ b/src/OpenFOAM/matrices/Matrix/MatrixI.H @@ -280,6 +280,14 @@ Foam::Matrix::col } +template +void Foam::Matrix::shallowResize(const label m, const label n) +{ + mRows_ = m; + nCols_ = n; +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H index 443ac00bd2..6baecc18e7 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -113,6 +113,9 @@ public: //- Resize the matrix preserving the elements inline void setSize(const label m); + //- Resize the matrix without reallocating storage (unsafe) + inline void shallowResize(const label m); + // Member operators diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H index 12ecd2530f..b85715ff07 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrixI.H @@ -161,6 +161,13 @@ inline void Foam::SquareMatrix::setSize(const label m) } +template +inline void Foam::SquareMatrix::shallowResize(const label m) +{ + Matrix, Type>::shallowResize(m, m); +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // template diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C index 2591f14f1a..f9acbe9606 100644 --- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C +++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -61,7 +61,14 @@ void Foam::ode::solve scalar& subDeltaT ) const { - label nSpecie = this->nSpecie(); + // Reset the size of the ODE system to the simplified size when mechanism + // reduction is active + if (odeSolver_->resize()) + { + odeSolver_->resizeField(cTp_); + } + + const label nSpecie = this->nSpecie(); // Copy the concentration, T and P to the total solve-vector for (int i=0; i odeSolver_; + + mutable autoPtr odeSolver_; // Solver data mutable scalarField cTp_;