Commit a5d73747 authored by Henry Weller's avatar Henry Weller
Browse files

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.
parent 47b9b3a8
......@@ -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,
......
......@@ -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,
......
......@@ -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,
......
......@@ -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,
......
......@@ -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,11 +34,33 @@ namespace Foam
}
// * * * * * * * * * * * * * Private 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;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict)
:
odes_(ode),
maxN_(ode.nEqns()),
n_(ode.nEqns()),
absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
......@@ -54,6 +76,7 @@ Foam::ODESolver::ODESolver
)
:
odes_(ode),
maxN_(ode.nEqns()),
n_(ode.nEqns()),
absTol_(absTol),
relTol_(relTol),
......@@ -63,22 +86,42 @@ Foam::ODESolver::ODESolver
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::ODESolver::normalizeError
(
const scalarField& y0,
const scalarField& y,
const scalarField& err
) const
bool Foam::ODESolver::resize()
{
// Calculate the maximum error
scalar maxErr = 0.0;
forAll(err, i)
if (odes_.nEqns() != n_)
{
scalar tol = absTol_[i] + relTol_[i]*max(mag(y0[i]), mag(y[i]));
maxErr = max(maxErr, mag(err[i])/tol);
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;
}
}
return maxErr;
void Foam::ODESolver::solve
(
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
stepState step(dxTry);
solve(x, y, step);
dxTry = step.dxTry;
}
......
......@@ -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;
//- 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<class Type>
static inline void resizeField(UList<Type>& f, const label n);
scalarField& relTol()
{
return relTol_;
}
template<class Type>
inline void resizeField(UList<Type>& 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
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * 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<class Type>
inline void Foam::ODESolver::resizeField(UList<Type>& f, const label n)
{
f.shallowCopy(UList<Type>(f.begin(), n));
}
template<class Type>
inline void Foam::ODESolver::resizeField(UList<Type>& f) const
{
resizeField(f, n_);
}
inline void Foam::ODESolver::resizeMatrix(scalarSquareMatrix& m) const
{
m.shallowResize(n_);
}
// ************************************************************************* //
......@@ -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,
......
......@@ -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,
......
......@@ -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,
......
......@@ -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,
......
......@@ -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,
......
......@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2013-2015 OpenFOAM Foundation