Commit 5b7dff49 authored by Henry's avatar Henry
Browse files

ODESolvers: Use the ODESystem protected data rather than pass redundant argument to solve

parent 63595880
......@@ -142,7 +142,7 @@ int main(int argc, char *argv[])
scalar dxNext = dxEst;
odeSolver->relTol() = relTol;
odeSolver->solve(ode, x, y, dxNext);
odeSolver->solve(x, y, dxNext);
Info<< scientific << setw(13) << relTol;
Info<< fixed << setw(11) << dxEst;
......@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
scalar dxEst = 0.5;
odeSolver->relTol() = 1e-4;
odeSolver->solve(ode, x, xEnd, y, dxEst);
odeSolver->solve(x, xEnd, y, dxEst);
Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
Info << "Numerical: y(2.0) = " << y << ", dxEst = " << dxEst << endl;
......
......@@ -49,7 +49,6 @@ Foam::Euler::Euler(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::Euler::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -75,13 +74,12 @@ Foam::scalar Foam::Euler::solve
void Foam::Euler::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -84,7 +84,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -95,7 +94,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -54,7 +54,6 @@ Foam::EulerSI::EulerSI(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::EulerSI::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -62,7 +61,7 @@ Foam::scalar Foam::EulerSI::solve
scalarField& y
) const
{
ode.jacobian(x0, y0, dfdx_, dfdy_);
odes_.jacobian(x0, y0, dfdx_, dfdy_);
for (register label i=0; i<n_; i++)
{
......@@ -95,13 +94,12 @@ Foam::scalar Foam::EulerSI::solve
void Foam::EulerSI::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -91,7 +91,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -102,7 +101,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -38,6 +38,7 @@ namespace Foam
Foam::ODESolver::ODESolver(const ODESystem& ode, const dictionary& dict)
:
odes_(ode),
n_(ode.nEqns()),
absTol_(n_, dict.lookupOrDefault<scalar>("absTol", SMALL)),
relTol_(n_, dict.lookupOrDefault<scalar>("relTol", 1e-4)),
......@@ -52,6 +53,7 @@ Foam::ODESolver::ODESolver
const scalarField& relTol
)
:
odes_(ode),
n_(ode.nEqns()),
absTol_(absTol),
relTol_(relTol),
......@@ -82,7 +84,6 @@ Foam::scalar Foam::ODESolver::normalizeError
void Foam::ODESolver::solve
(
const ODESystem& ode,
const scalar xStart,
const scalar xEnd,
scalarField& y,
......@@ -105,7 +106,7 @@ void Foam::ODESolver::solve
}
// Integrate as far as possible up to dxEst
solve(ode, x, y, dxEst);
solve(x, y, dxEst);
// Check if reached xEnd
if ((x - xEnd)*(xEnd - xStart) >= 0)
......@@ -122,7 +123,7 @@ void Foam::ODESolver::solve
FatalErrorIn
(
"ODESolver::solve"
"(const ODESystem& ode, const scalar xStart, const scalar xEnd,"
"(const scalar xStart, const scalar xEnd,"
"scalarField& y, scalar& dxEst) const"
) << "Integration steps greater than maximum " << maxSteps_
<< exit(FatalError);
......
......@@ -55,6 +55,9 @@ protected:
// Protected data
//- Reference to ODESystem
const ODESystem& odes_;
//- Size of the ODESystem
label n_;
......@@ -150,7 +153,6 @@ public:
// Update the state and return an estimate for the next step in dxTry
virtual void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......@@ -160,7 +162,6 @@ public:
// and return an estimate for the next step in dxTry
virtual void solve
(
const ODESystem& ode,
const scalar xStart,
const scalar xEnd,
scalarField& y,
......
......@@ -89,7 +89,6 @@ Foam::RKCK45::RKCK45(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::RKCK45::solve
(
const ODESystem& odes,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -102,21 +101,21 @@ Foam::scalar Foam::RKCK45::solve
yTemp_[i] = y0[i] + a21*dx*dydx0[i];
}
odes.derivatives(x0 + c2*dx, yTemp_, k2_);
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_);
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_);
odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
forAll(yTemp_, i)
{
......@@ -124,7 +123,7 @@ Foam::scalar Foam::RKCK45::solve
+ dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
}
odes.derivatives(x0 + c5*dx, yTemp_, k5_);
odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
forAll(yTemp_, i)
{
......@@ -133,7 +132,7 @@ Foam::scalar Foam::RKCK45::solve
*(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
}
odes.derivatives(x0 + c6*dx, yTemp_, k6_);
odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
forAll(y, i)
{
......@@ -154,13 +153,12 @@ Foam::scalar Foam::RKCK45::solve
void Foam::RKCK45::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -110,7 +110,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -121,7 +120,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -94,7 +94,6 @@ Foam::RKDP45::RKDP45(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::RKDP45::solve
(
const ODESystem& odes,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -107,21 +106,21 @@ Foam::scalar Foam::RKDP45::solve
yTemp_[i] = y0[i] + a21*dx*dydx0[i];
}
odes.derivatives(x0 + c2*dx, yTemp_, k2_);
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_);
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_);
odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
forAll(yTemp_, i)
{
......@@ -129,7 +128,7 @@ Foam::scalar Foam::RKDP45::solve
+ dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
}
odes.derivatives(x0 + c5*dx, yTemp_, k5_);
odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
forAll(yTemp_, i)
{
......@@ -138,7 +137,7 @@ Foam::scalar Foam::RKDP45::solve
*(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
}
odes.derivatives(x0 + dx, yTemp_, k6_);
odes_.derivatives(x0 + dx, yTemp_, k6_);
forAll(y, i)
{
......@@ -147,7 +146,7 @@ Foam::scalar Foam::RKDP45::solve
}
// Reuse k2_ for the derivative of the new state
odes.derivatives(x0 + dx, y, k2_);
odes_.derivatives(x0 + dx, y, k2_);
forAll(err_, i)
{
......@@ -163,13 +162,12 @@ Foam::scalar Foam::RKDP45::solve
void Foam::RKDP45::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -114,7 +114,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -125,7 +124,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -90,7 +90,6 @@ Foam::RKF45::RKF45(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::RKF45::solve
(
const ODESystem& odes,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -103,21 +102,21 @@ Foam::scalar Foam::RKF45::solve
yTemp_[i] = y0[i] + a21*dx*dydx0[i];
}
odes.derivatives(x0 + c2*dx, yTemp_, k2_);
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_);
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_);
odes_.derivatives(x0 + c4*dx, yTemp_, k4_);
forAll(yTemp_, i)
{
......@@ -125,7 +124,7 @@ Foam::scalar Foam::RKF45::solve
+ dx*(a51*dydx0[i] + a52*k2_[i] + a53*k3_[i] + a54*k4_[i]);
}
odes.derivatives(x0 + c5*dx, yTemp_, k5_);
odes_.derivatives(x0 + c5*dx, yTemp_, k5_);
forAll(yTemp_, i)
{
......@@ -134,7 +133,7 @@ Foam::scalar Foam::RKF45::solve
*(a61*dydx0[i] + a62*k2_[i] + a63*k3_[i] + a64*k4_[i] + a65*k5_[i]);
}
odes.derivatives(x0 + c6*dx, yTemp_, k6_);
odes_.derivatives(x0 + c6*dx, yTemp_, k6_);
// Calculate the 5th-order solution
forAll(y, i)
......@@ -159,13 +158,12 @@ Foam::scalar Foam::RKF45::solve
void Foam::RKF45::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -138,7 +138,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -149,7 +148,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -68,7 +68,6 @@ Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::Rosenbrock21::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -76,7 +75,7 @@ Foam::scalar Foam::Rosenbrock21::solve
scalarField& y
) const
{
ode.jacobian(x0, y0, dfdx_, dfdy_);
odes_.jacobian(x0, y0, dfdx_, dfdy_);
for (register label i=0; i<n_; i++)
{
......@@ -104,7 +103,7 @@ Foam::scalar Foam::Rosenbrock21::solve
y[i] = y0[i] + a21*k1_[i];
}
ode.derivatives(x0 + c2*dx, y, dydx_);
odes_.derivatives(x0 + c2*dx, y, dydx_);
forAll(k2_, i)
{
......@@ -126,13 +125,12 @@ Foam::scalar Foam::Rosenbrock21::solve
void Foam::Rosenbrock21::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -101,7 +101,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -112,7 +111,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -81,7 +81,6 @@ Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::Rosenbrock32::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -89,7 +88,7 @@ Foam::scalar Foam::Rosenbrock32::solve
scalarField& y
) const
{
ode.jacobian(x0, y0, dfdx_, dfdy_);
odes_.jacobian(x0, y0, dfdx_, dfdy_);
for (register label i=0; i<n_; i++)
{
......@@ -117,7 +116,7 @@ Foam::scalar Foam::Rosenbrock32::solve
y[i] = y0[i] + a21*k1_[i];
}
ode.derivatives(x0 + c2*dx, y, dydx_);
odes_.derivatives(x0 + c2*dx, y, dydx_);
forAll(k2_, i)
{
......@@ -148,13 +147,12 @@ Foam::scalar Foam::Rosenbrock32::solve
void Foam::Rosenbrock32::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -105,7 +105,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -116,7 +115,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -126,7 +126,6 @@ Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict)
Foam::scalar Foam::Rosenbrock43::solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -134,7 +133,7 @@ Foam::scalar Foam::Rosenbrock43::solve
scalarField& y
) const
{
ode.jacobian(x0, y0, dfdx_, dfdy_);
odes_.jacobian(x0, y0, dfdx_, dfdy_);
for (register label i=0; i<n_; i++)
{
......@@ -162,7 +161,7 @@ Foam::scalar Foam::Rosenbrock43::solve
y[i] = y0[i] + a21*k1_[i];
}
ode.derivatives(x0 + c2*dx, y, dydx_);
odes_.derivatives(x0 + c2*dx, y, dydx_);
forAll(k2_, i)
{
......@@ -177,7 +176,7 @@ Foam::scalar Foam::Rosenbrock43::solve
y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
}
ode.derivatives(x0 + c3*dx, y, dydx_);
odes_.derivatives(x0 + c3*dx, y, dydx_);
forAll(k3_, i)
{
......@@ -208,13 +207,12 @@ Foam::scalar Foam::Rosenbrock43::solve
void Foam::Rosenbrock43::solve
(
const ODESystem& odes,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
adaptiveSolver::solve(odes, x, y, dxTry);
adaptiveSolver::solve(odes_, x, y, dxTry);
}
......
......@@ -112,7 +112,6 @@ public:
//- Solve a single step dx and return the error
scalar solve
(
const ODESystem& ode,
const scalar x0,
const scalarField& y0,
const scalarField& dydx0,
......@@ -123,7 +122,6 @@ public:
//- Solve the ODE system and the update the state
void solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
......
......@@ -70,13 +70,12 @@ Foam::SIBS::SIBS(const ODESystem& ode, const dictionary& dict)
void Foam::SIBS::solve
(
const ODESystem& ode,
scalar& x,
scalarField& y,
scalar& dxTry
) const
{
ode.derivatives(x, y, dydx0_);
odes_.derivatives(x, y, dydx0_);
scalar h = dxTry;