diff --git a/applications/test/ODE/Test-ODE.C b/applications/test/ODE/Test-ODE.C index 4222fa4cc4..ae20374259 100644 --- a/applications/test/ODE/Test-ODE.C +++ b/applications/test/ODE/Test-ODE.C @@ -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; diff --git a/src/ODE/Make/files b/src/ODE/Make/files index 7c361626ed..912b36fecb 100644 --- a/src/ODE/Make/files +++ b/src/ODE/Make/files @@ -8,11 +8,11 @@ ODESolvers/Trapezoid/Trapezoid.C ODESolvers/RKF45/RKF45.C ODESolvers/RKCK45/RKCK45.C ODESolvers/RKDP45/RKDP45.C -ODESolvers/Rosenbrock21/Rosenbrock21.C -ODESolvers/Rosenbrock32/Rosenbrock32.C -ODESolvers/Rosenbrock43/Rosenbrock43.C -ODESolvers/rodas32/rodas32.C -ODESolvers/rodas43/rodas43.C +ODESolvers/Rosenbrock12/Rosenbrock12.C +ODESolvers/Rosenbrock23/Rosenbrock23.C +ODESolvers/Rosenbrock34/Rosenbrock34.C +ODESolvers/rodas23/rodas23.C +ODESolvers/rodas34/rodas34.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 index 1a6edf576c..1496099ae8 100644 --- a/src/ODE/ODESolvers/Euler/Euler.C +++ b/src/ODE/ODESolvers/Euler/Euler.C @@ -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); } diff --git a/src/ODE/ODESolvers/Euler/Euler.H b/src/ODE/ODESolvers/Euler/Euler.H index 90f67f4b4c..3cefc425be 100644 --- a/src/ODE/ODESolvers/Euler/Euler.H +++ b/src/ODE/ODESolvers/Euler/Euler.H @@ -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 diff --git a/src/ODE/ODESolvers/EulerSI/EulerSI.C b/src/ODE/ODESolvers/EulerSI/EulerSI.C index 1785101758..9b28f733ea 100644 --- a/src/ODE/ODESolvers/EulerSI/EulerSI.C +++ b/src/ODE/ODESolvers/EulerSI/EulerSI.C @@ -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("absTol", SMALL)), relTol_(n_, dict.lookupOrDefault("relTol", 1e-4)), @@ -52,6 +53,7 @@ Foam::ODESolver::ODESolver const scalarField& relTol ) : + odes_(ode), n_(ode.nEqns()), absTol_(absTol), relTol_(relTol), @@ -82,48 +84,70 @@ Foam::scalar Foam::ODESolver::normalizeError void Foam::ODESolver::solve ( - const ODESystem& ode, + scalar& x, + scalarField& y, + stepState& step +) const +{ + solve(x, y, step.dxTry); +} + + +void Foam::ODESolver::solve +( const scalar xStart, const scalar xEnd, scalarField& y, - scalar& dxEst + scalar& dxTry ) const { + stepState step(dxTry); scalar x = xStart; - bool truncated = false; for (label nStep=0; nStep 0) + step.reject = false; + + // Check if this is a truncated step and set dxTry to integrate to xEnd + if ((x + step.dxTry - xEnd)*(x + step.dxTry - xStart) > 0) { - truncated = true; - dxEst = xEnd - x; + step.last = true; + step.dxTry = xEnd - x; } - // Integrate as far as possible up to dxEst - solve(ode, x, y, dxEst); + // Integrate as far as possible up to step.dxTry + solve(x, y, step); // Check if reached xEnd if ((x - xEnd)*(xEnd - xStart) >= 0) { - if (nStep > 0 && truncated) + if (nStep > 0 && step.last) { - dxEst = dxEst0; + step.dxTry = dxTry0; } + dxTry = step.dxTry; + return; } + + step.first = false; + + // If the step.dxTry was reject set step.prevReject + if (step.reject) + { + step.prevReject = true; + } } FatalErrorIn ( "ODESolver::solve" - "(const ODESystem& ode, const scalar xStart, const scalar xEnd," - "scalarField& y, scalar& dxEst) const" + "(const scalar xStart, const scalar xEnd," + "scalarField& y, scalar& dxTry) 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 bd767a9063..de4ef399f8 100644 --- a/src/ODE/ODESolvers/ODESolver/ODESolver.H +++ b/src/ODE/ODESolvers/ODESolver/ODESolver.H @@ -55,6 +55,9 @@ protected: // Protected data + //- Reference to ODESystem + const ODESystem& odes_; + //- Size of the ODESystem label n_; @@ -90,6 +93,30 @@ public: //- Runtime type information TypeName("ODESolver"); + class stepState + { + public: + + const bool forward; + scalar dxTry; + scalar dxDid; + bool first; + bool last; + bool reject; + bool prevReject; + + stepState(const scalar dx) + : + forward(dx > 0 ? true : false), + dxTry(dx), + dxDid(0), + first(true), + last(false), + reject(false), + prevReject(false) + {} + }; + // Declare run-time constructor selection table @@ -150,17 +177,31 @@ 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 - ) const = 0; + ) const //= 0; + { + stepState step(dxTry); + solve(x, y, step); + dxTry = step.dxTry; + } + + //- 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 + ( + scalar& x, + scalarField& y, + stepState& step + ) const; //- 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, diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.C b/src/ODE/ODESolvers/RKCK45/RKCK45.C index 2de3cdf99e..30ae91de1a 100644 --- a/src/ODE/ODESolvers/RKCK45/RKCK45.C +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.C @@ -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); } diff --git a/src/ODE/ODESolvers/RKCK45/RKCK45.H b/src/ODE/ODESolvers/RKCK45/RKCK45.H index a6b551a475..5402c6efbd 100644 --- a/src/ODE/ODESolvers/RKCK45/RKCK45.H +++ b/src/ODE/ODESolvers/RKCK45/RKCK45.H @@ -34,10 +34,7 @@ Description Cash, J.R., Karp, A.H. 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., @@ -110,7 +107,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 +117,6 @@ public: //- Solve the ODE system and the update the state void solve ( - const ODESystem& ode, scalar& x, scalarField& y, scalar& dxTry diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.C b/src/ODE/ODESolvers/RKDP45/RKDP45.C index cf319dca8b..f9a3a78836 100644 --- a/src/ODE/ODESolvers/RKDP45/RKDP45.C +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.C @@ -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); } diff --git a/src/ODE/ODESolvers/RKDP45/RKDP45.H b/src/ODE/ODESolvers/RKDP45/RKDP45.H index 95ae212e66..e9ac93592a 100644 --- a/src/ODE/ODESolvers/RKDP45/RKDP45.H +++ b/src/ODE/ODESolvers/RKDP45/RKDP45.H @@ -34,10 +34,7 @@ Description 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., @@ -114,7 +111,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 +121,6 @@ public: //- Solve the ODE system and the update the state void solve ( - const ODESystem& ode, scalar& x, scalarField& y, scalar& dxTry diff --git a/src/ODE/ODESolvers/RKF45/RKF45.C b/src/ODE/ODESolvers/RKF45/RKF45.C index 661c618255..92804e0634 100644 --- a/src/ODE/ODESolvers/RKF45/RKF45.C +++ b/src/ODE/ODESolvers/RKF45/RKF45.C @@ -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); } diff --git a/src/ODE/ODESolvers/RKF45/RKF45.H b/src/ODE/ODESolvers/RKF45/RKF45.H index df91ed7808..7567adc959 100644 --- a/src/ODE/ODESolvers/RKF45/RKF45.H +++ b/src/ODE/ODESolvers/RKF45/RKF45.H @@ -34,8 +34,6 @@ Description Fehlberg, E., NASA Technical Report 315, 1969. - Based on code from: - \verbatim "Solving Ordinary Differential Equations I: Nonstiff Problems, second edition", Hairer, E., @@ -48,31 +46,6 @@ Description 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 @@ -138,7 +111,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 +121,6 @@ public: //- Solve the ODE system and the update the state void solve ( - const ODESystem& ode, scalar& x, scalarField& y, scalar& dxTry diff --git a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C similarity index 77% rename from src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C rename to src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C index 13371a8e90..f69e6a9cde 100644 --- a/src/ODE/ODESolvers/Rosenbrock21/Rosenbrock21.C +++ b/src/ODE/ODESolvers/Rosenbrock12/Rosenbrock12.C @@ -23,33 +23,33 @@ License \*---------------------------------------------------------------------------*/ -#include "Rosenbrock21.H" +#include "Rosenbrock12.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(Rosenbrock21, 0); - addToRunTimeSelectionTable(ODESolver, Rosenbrock21, dictionary); + defineTypeNameAndDebug(Rosenbrock12, 0); + addToRunTimeSelectionTable(ODESolver, Rosenbrock12, dictionary); const scalar - 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 = gamma, - Rosenbrock21::d2 = -gamma; + Rosenbrock12::gamma = 1 + 1.0/std::sqrt(2.0), + Rosenbrock12::a21 = 1.0/gamma, + Rosenbrock12::c2 = 1.0, + Rosenbrock12::c21 = -2.0/gamma, + Rosenbrock12::b1 = (3.0/2.0)/gamma, + Rosenbrock12::b2 = (1.0/2.0)/gamma, + Rosenbrock12::e1 = b1 - 1.0/gamma, + Rosenbrock12::e2 = b2, + Rosenbrock12::d1 = gamma, + Rosenbrock12::d2 = -gamma; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict) +Foam::Rosenbrock12::Rosenbrock12(const ODESystem& ode, const dictionary& dict) : ODESolver(ode, dict), adaptiveSolver(ode, dict), @@ -66,9 +66,8 @@ Foam::Rosenbrock21::Rosenbrock21(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::Rosenbrock21::solve +Foam::scalar Foam::Rosenbrock12::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. Class - Foam::Rosenbrock21 + Foam::Rosenbrock12 Description L-stable embedded Rosenbrock ODE solver of order (1)2. @@ -39,12 +39,12 @@ Description \endverbatim SourceFiles - Rosenbrock21.C + Rosenbrock12.C \*---------------------------------------------------------------------------*/ -#ifndef Rosenbrock21_H -#define Rosenbrock21_H +#ifndef Rosenbrock12_H +#define Rosenbrock12_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -55,10 +55,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class Rosenbrock21 Declaration + Class Rosenbrock12 Declaration \*---------------------------------------------------------------------------*/ -class Rosenbrock21 +class Rosenbrock12 : public ODESolver, public adaptiveSolver @@ -87,13 +87,13 @@ class Rosenbrock21 public: //- Runtime type information - TypeName("Rosenbrock21"); + TypeName("Rosenbrock12"); // Constructors //- Construct from ODE - Rosenbrock21(const ODESystem& ode, const dictionary& dict); + Rosenbrock12(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -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 diff --git a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C similarity index 69% rename from src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C rename to src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C index e471e3ffe1..035502714e 100644 --- a/src/ODE/ODESolvers/Rosenbrock32/Rosenbrock32.C +++ b/src/ODE/ODESolvers/Rosenbrock23/Rosenbrock23.C @@ -23,45 +23,45 @@ License \*---------------------------------------------------------------------------*/ -#include "Rosenbrock32.H" +#include "Rosenbrock23.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(Rosenbrock32, 0); - addToRunTimeSelectionTable(ODESolver, Rosenbrock32, dictionary); + defineTypeNameAndDebug(Rosenbrock23, 0); + addToRunTimeSelectionTable(ODESolver, Rosenbrock23, dictionary); const scalar - Rosenbrock32::a21 = 1, - Rosenbrock32::a31 = 1, - Rosenbrock32::a32 = 0, + Rosenbrock23::a21 = 1, + Rosenbrock23::a31 = 1, + Rosenbrock23::a32 = 0, - Rosenbrock32::c21 = -1.0156171083877702091975600115545, - Rosenbrock32::c31 = 4.0759956452537699824805835358067, - Rosenbrock32::c32 = 9.2076794298330791242156818474003, + Rosenbrock23::c21 = -1.0156171083877702091975600115545, + Rosenbrock23::c31 = 4.0759956452537699824805835358067, + Rosenbrock23::c32 = 9.2076794298330791242156818474003, - Rosenbrock32::b1 = 1, - Rosenbrock32::b2 = 6.1697947043828245592553615689730, - Rosenbrock32::b3 = -0.4277225654321857332623837380651, + Rosenbrock23::b1 = 1, + Rosenbrock23::b2 = 6.1697947043828245592553615689730, + Rosenbrock23::b3 = -0.4277225654321857332623837380651, - Rosenbrock32::e1 = 0.5, - Rosenbrock32::e2 = -2.9079558716805469821718236208017, - Rosenbrock32::e3 = 0.2235406989781156962736090927619, + Rosenbrock23::e1 = 0.5, + Rosenbrock23::e2 = -2.9079558716805469821718236208017, + Rosenbrock23::e3 = 0.2235406989781156962736090927619, - Rosenbrock32::gamma = 0.43586652150845899941601945119356, - Rosenbrock32::c2 = 0.43586652150845899941601945119356, + Rosenbrock23::gamma = 0.43586652150845899941601945119356, + Rosenbrock23::c2 = 0.43586652150845899941601945119356, - Rosenbrock32::d1 = 0.43586652150845899941601945119356, - Rosenbrock32::d2 = 0.24291996454816804366592249683314, - Rosenbrock32::d3 = 2.1851380027664058511513169485832; + Rosenbrock23::d1 = 0.43586652150845899941601945119356, + Rosenbrock23::d2 = 0.24291996454816804366592249683314, + Rosenbrock23::d3 = 2.1851380027664058511513169485832; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict) +Foam::Rosenbrock23::Rosenbrock23(const ODESystem& ode, const dictionary& dict) : ODESolver(ode, dict), adaptiveSolver(ode, dict), @@ -79,9 +79,8 @@ Foam::Rosenbrock32::Rosenbrock32(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::Rosenbrock32::solve +Foam::scalar Foam::Rosenbrock23::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. Class - Foam::Rosenbrock32 + Foam::Rosenbrock23 Description L-stable embedded Rosenbrock ODE solver of order (3)4. @@ -42,12 +42,12 @@ Description \endverbatim SourceFiles - Rosenbrock32.C + Rosenbrock23.C \*---------------------------------------------------------------------------*/ -#ifndef Rosenbrock32_H -#define Rosenbrock32_H +#ifndef Rosenbrock23_H +#define Rosenbrock23_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -58,10 +58,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class Rosenbrock32 Declaration + Class Rosenbrock23 Declaration \*---------------------------------------------------------------------------*/ -class Rosenbrock32 +class Rosenbrock23 : public ODESolver, public adaptiveSolver @@ -91,13 +91,13 @@ class Rosenbrock32 public: //- Runtime type information - TypeName("Rosenbrock32"); + TypeName("Rosenbrock23"); // Constructors //- Construct from ODE - Rosenbrock32(const ODESystem& ode, const dictionary& dict); + Rosenbrock23(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -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 diff --git a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C similarity index 60% rename from src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C rename to src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C index 061390fdbc..65bf3ced87 100644 --- a/src/ODE/ODESolvers/Rosenbrock43/Rosenbrock43.C +++ b/src/ODE/ODESolvers/Rosenbrock34/Rosenbrock34.C @@ -23,89 +23,89 @@ License \*---------------------------------------------------------------------------*/ -#include "Rosenbrock43.H" +#include "Rosenbrock34.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(Rosenbrock43, 0); - addToRunTimeSelectionTable(ODESolver, Rosenbrock43, dictionary); + defineTypeNameAndDebug(Rosenbrock34, 0); + addToRunTimeSelectionTable(ODESolver, Rosenbrock34, dictionary); const scalar - // L-Stable constants from Hairer et. al. - Rosenbrock43::a21 = 2, - Rosenbrock43::a31 = 1.867943637803922, - Rosenbrock43::a32 = 0.2344449711399156, - - Rosenbrock43::c21 = -7.137615036412310, - Rosenbrock43::c31 = 2.580708087951457, - Rosenbrock43::c32 = 0.6515950076447975, - Rosenbrock43::c41 = -2.137148994382534, - Rosenbrock43::c42 = -0.3214669691237626, - Rosenbrock43::c43 = -0.6949742501781779, - - Rosenbrock43::b1 = 2.255570073418735, - Rosenbrock43::b2 = 0.2870493262186792, - Rosenbrock43::b3 = 0.435317943184018, - Rosenbrock43::b4 = 1.093502252409163, - - Rosenbrock43::e1 = -0.2815431932141155, - Rosenbrock43::e2 = -0.0727619912493892, - Rosenbrock43::e3 = -0.1082196201495311, - Rosenbrock43::e4 = -1.093502252409163, - - Rosenbrock43::gamma = 0.57282, - Rosenbrock43::c2 = 1.14564, - Rosenbrock43::c3 = 0.65521686381559, - - Rosenbrock43::d1 = 0.57282, - Rosenbrock43::d2 = -1.769193891319233, - Rosenbrock43::d3 = 0.7592633437920482, - Rosenbrock43::d4 = -0.1049021087100450; - // Constants by Shampine // More accurate than the L-Stable coefficients for small step-size // but less stable for large step-size + Rosenbrock34::a21 = 2, + Rosenbrock34::a31 = 48.0/25.0, + Rosenbrock34::a32 = 6.0/25.0, + + Rosenbrock34::c21 = -8, + Rosenbrock34::c31 = 372.0/25.0, + Rosenbrock34::c32 = 12.0/5.0, + + Rosenbrock34::c41 = -112.0/125.0, + Rosenbrock34::c42 = -54.0/125.0, + Rosenbrock34::c43 = -2.0/5.0, + + Rosenbrock34::b1 = 19.0/9.0, + Rosenbrock34::b2 = 1.0/2.0, + Rosenbrock34::b3 = 25.0/108.0, + Rosenbrock34::b4 = 125.0/108.0, + + Rosenbrock34::e1 = 34.0/108.0, + Rosenbrock34::e2 = 7.0/36.0, + Rosenbrock34::e3 = 0, + Rosenbrock34::e4 = 125.0/108.0, + + Rosenbrock34::gamma = 1.0/2.0, + Rosenbrock34::c2 = 1, + Rosenbrock34::c3 = 3.0/5.0, + + Rosenbrock34::d1 = 1.0/2.0, + Rosenbrock34::d2 = -3.0/2.0, + Rosenbrock34::d3 = 605.0/250.0, + Rosenbrock34::d4 = 29.0/250.0; + /* - Rosenbrock43::a21 = 2, - Rosenbrock43::a31 = 48.0/25.0, - Rosenbrock43::a32 = 6.0/25.0, + // L-Stable constants from Hairer et. al. + Rosenbrock34::a21 = 2, + Rosenbrock34::a31 = 1.867943637803922, + Rosenbrock34::a32 = 0.2344449711399156, - Rosenbrock43::c21 = -8, - Rosenbrock43::c31 = 372.0/25.0, - Rosenbrock43::c32 = 12.0/5.0, + Rosenbrock34::c21 = -7.137615036412310, + Rosenbrock34::c31 = 2.580708087951457, + Rosenbrock34::c32 = 0.6515950076447975, + Rosenbrock34::c41 = -2.137148994382534, + Rosenbrock34::c42 = -0.3214669691237626, + Rosenbrock34::c43 = -0.6949742501781779, - Rosenbrock43::c41 = -112.0/125.0, - Rosenbrock43::c42 = -54.0/125.0, - Rosenbrock43::c43 = -2.0/5.0, + Rosenbrock34::b1 = 2.255570073418735, + Rosenbrock34::b2 = 0.2870493262186792, + Rosenbrock34::b3 = 0.435317943184018, + Rosenbrock34::b4 = 1.093502252409163, - Rosenbrock43::b1 = 19.0/9.0, - Rosenbrock43::b2 = 1.0/2.0, - Rosenbrock43::b3 = 25.0/108.0, - Rosenbrock43::b4 = 125.0/108.0, + Rosenbrock34::e1 = -0.2815431932141155, + Rosenbrock34::e2 = -0.0727619912493892, + Rosenbrock34::e3 = -0.1082196201495311, + Rosenbrock34::e4 = -1.093502252409163, - Rosenbrock43::e1 = 34.0/108.0, - Rosenbrock43::e2 = 7.0/36.0, - Rosenbrock43::e3 = 0, - Rosenbrock43::e4 = 125.0/108.0, + Rosenbrock34::gamma = 0.57282, + Rosenbrock34::c2 = 1.14564, + Rosenbrock34::c3 = 0.65521686381559, - Rosenbrock43::gamma = 1.0/2.0, - Rosenbrock43::c2 = 1, - 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; + Rosenbrock34::d1 = 0.57282, + Rosenbrock34::d2 = -1.769193891319233, + Rosenbrock34::d3 = 0.7592633437920482, + Rosenbrock34::d4 = -0.1049021087100450; */ } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict) +Foam::Rosenbrock34::Rosenbrock34(const ODESystem& ode, const dictionary& dict) : ODESolver(ode, dict), adaptiveSolver(ode, dict), @@ -124,9 +124,8 @@ Foam::Rosenbrock43::Rosenbrock43(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::Rosenbrock43::solve +Foam::scalar Foam::Rosenbrock34::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. Class - Foam::Rosenbrock43 + Foam::Rosenbrock34 Description L-stable embedded Rosenbrock ODE solver of order (3)4. - Based on code from: \verbatim "Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, second edition", @@ -37,7 +36,7 @@ Description Springer-Verlag, Berlin. 1996. \endverbatim - Also provided are the constants from: + The default constants are from: \verbatim "Implementation of Rosenbrock Methods" Shampine, L. F., @@ -46,13 +45,15 @@ Description with which the scheme is more accurate than with the L-Stable coefficients for small step-size but less stable for large step-size. + The L-Stable scheme constants are provided commented-out in Rosenbrock34.C + SourceFiles - Rosenbrock43.C + Rosenbrock34.C \*---------------------------------------------------------------------------*/ -#ifndef Rosenbrock43_H -#define Rosenbrock43_H +#ifndef Rosenbrock34_H +#define Rosenbrock34_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -63,10 +64,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class Rosenbrock43 Declaration + Class Rosenbrock34 Declaration \*---------------------------------------------------------------------------*/ -class Rosenbrock43 +class Rosenbrock34 : public ODESolver, public adaptiveSolver @@ -98,13 +99,13 @@ class Rosenbrock43 public: //- Runtime type information - TypeName("Rosenbrock43"); + TypeName("Rosenbrock34"); // Constructors //- Construct from ODE - Rosenbrock43(const ODESystem& ode, const dictionary& dict); + Rosenbrock34(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -112,7 +113,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 +123,6 @@ public: //- Solve the ODE system and the update the state void solve ( - const ODESystem& ode, scalar& x, scalarField& y, scalar& dxTry diff --git a/src/ODE/ODESolvers/SIBS/SIBS.C b/src/ODE/ODESolvers/SIBS/SIBS.C index eae8552a47..b6b6cf2c1d 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.C +++ b/src/ODE/ODESolvers/SIBS/SIBS.C @@ -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; bool exitflag = false; @@ -124,7 +123,7 @@ void Foam::SIBS::solve label k = 0; yTemp_ = y; - ode.jacobian(x, y, dfdx_, dfdy_); + odes_.jacobian(x, y, dfdx_, dfdy_); if (x != xNew_ || h != dxTry) { @@ -151,7 +150,7 @@ void Foam::SIBS::solve << exit(FatalError); } - SIMPR(ode, x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_); + SIMPR(x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_); scalar xest = sqr(h/nSeq_[k]); polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_); diff --git a/src/ODE/ODESolvers/SIBS/SIBS.H b/src/ODE/ODESolvers/SIBS/SIBS.H index 03664fd53a..0ef0a1d40c 100644 --- a/src/ODE/ODESolvers/SIBS/SIBS.H +++ b/src/ODE/ODESolvers/SIBS/SIBS.H @@ -84,7 +84,6 @@ class SIBS void SIMPR ( - const ODESystem& ode, const scalar xStart, const scalarField& y, const scalarField& dydx, @@ -123,7 +122,6 @@ public: void solve ( - const ODESystem& ode, scalar& x, scalarField& y, scalar& dxTry diff --git a/src/ODE/ODESolvers/SIBS/SIMPR.C b/src/ODE/ODESolvers/SIBS/SIMPR.C index 36402a202d..d21a931e1d 100644 --- a/src/ODE/ODESolvers/SIBS/SIMPR.C +++ b/src/ODE/ODESolvers/SIBS/SIMPR.C @@ -29,7 +29,6 @@ License void Foam::SIBS::SIMPR ( - const ODESystem& ode, const scalar xStart, const scalarField& y, const scalarField& dydx, @@ -72,7 +71,7 @@ void Foam::SIBS::SIMPR scalar x = xStart + h; - ode.derivatives(x, ytemp, yEnd); + odes_.derivatives(x, ytemp, yEnd); for (register label nn=2; nn<=nSteps; nn++) { @@ -90,7 +89,7 @@ void Foam::SIBS::SIMPR x += h; - ode.derivatives(x, ytemp, yEnd); + odes_.derivatives(x, ytemp, yEnd); } for (register label i=0; i 1) diff --git a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H index d86bb88ce1..0d7dcb035c 100644 --- a/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H +++ b/src/ODE/ODESolvers/adaptiveSolver/adaptiveSolver.H @@ -77,7 +77,6 @@ public: //- Solve a single step dx and return the error virtual scalar solve ( - const ODESystem& ode, const scalar x0, const scalarField& y0, const scalarField& dydx0, diff --git a/src/ODE/ODESolvers/rodas32/rodas32.C b/src/ODE/ODESolvers/rodas23/rodas23.C similarity index 80% rename from src/ODE/ODESolvers/rodas32/rodas32.C rename to src/ODE/ODESolvers/rodas23/rodas23.C index 6221788ea6..e37bc3b00e 100644 --- a/src/ODE/ODESolvers/rodas32/rodas32.C +++ b/src/ODE/ODESolvers/rodas23/rodas23.C @@ -23,35 +23,35 @@ License \*---------------------------------------------------------------------------*/ -#include "rodas32.H" +#include "rodas23.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(rodas32, 0); - addToRunTimeSelectionTable(ODESolver, rodas32, dictionary); + defineTypeNameAndDebug(rodas23, 0); + addToRunTimeSelectionTable(ODESolver, rodas23, dictionary); const scalar - rodas32::c3 = 1, - rodas32::d1 = 1.0/2.0, - rodas32::d2 = 3.0/2.0, - rodas32::a31 = 2, - rodas32::a41 = 2, - rodas32::c21 = 4, - rodas32::c31 = 1, - rodas32::c32 = -1, - rodas32::c41 = 1, - rodas32::c42 = -1, - rodas32::c43 = -8.0/3.0, - rodas32::gamma = 1.0/2.0; + rodas23::c3 = 1, + rodas23::d1 = 1.0/2.0, + rodas23::d2 = 3.0/2.0, + rodas23::a31 = 2, + rodas23::a41 = 2, + rodas23::c21 = 4, + rodas23::c31 = 1, + rodas23::c32 = -1, + rodas23::c41 = 1, + rodas23::c42 = -1, + rodas23::c43 = -8.0/3.0, + rodas23::gamma = 1.0/2.0; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict) +Foam::rodas23::rodas23(const ODESystem& ode, const dictionary& dict) : ODESolver(ode, dict), adaptiveSolver(ode, dict), @@ -70,9 +70,8 @@ Foam::rodas32::rodas32(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::rodas32::solve +Foam::scalar Foam::rodas23::solve ( - const ODESystem& ode, const scalar x0, const scalarField& y0, const scalarField& dydx0, @@ -80,7 +79,7 @@ Foam::scalar Foam::rodas32::solve scalarField& y ) const { - ode.jacobian(x0, y0, dfdx_, dfdy_); + odes_.jacobian(x0, y0, dfdx_, dfdy_); for (register label i=0; i. Class - Foam::rodas32 + Foam::rodas23 Description L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (2)3. @@ -42,12 +42,12 @@ Description \endverbatim SourceFiles - rodas32.C + rodas23.C \*---------------------------------------------------------------------------*/ -#ifndef rodas32_H -#define rodas32_H +#ifndef rodas23_H +#define rodas23_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -58,10 +58,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class rodas32 Declaration + Class rodas23 Declaration \*---------------------------------------------------------------------------*/ -class rodas32 +class rodas23 : public ODESolver, public adaptiveSolver @@ -91,13 +91,13 @@ class rodas32 public: //- Runtime type information - TypeName("rodas32"); + TypeName("rodas23"); // Constructors //- Construct from ODE - rodas32(const ODESystem& ode, const dictionary& dict); + rodas23(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -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 diff --git a/src/ODE/ODESolvers/rodas43/rodas43.C b/src/ODE/ODESolvers/rodas34/rodas34.C similarity index 67% rename from src/ODE/ODESolvers/rodas43/rodas43.C rename to src/ODE/ODESolvers/rodas34/rodas34.C index 92b78a7ab4..5d3ca9ac35 100644 --- a/src/ODE/ODESolvers/rodas43/rodas43.C +++ b/src/ODE/ODESolvers/rodas34/rodas34.C @@ -23,56 +23,56 @@ License \*---------------------------------------------------------------------------*/ -#include "rodas43.H" +#include "rodas34.H" #include "addToRunTimeSelectionTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(rodas43, 0); - addToRunTimeSelectionTable(ODESolver, rodas43, dictionary); + defineTypeNameAndDebug(rodas34, 0); + addToRunTimeSelectionTable(ODESolver, rodas34, dictionary); const scalar - rodas43::c2 = 0.386, - rodas43::c3 = 0.21, - rodas43::c4 = 0.63, - 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; + rodas34::c2 = 0.386, + rodas34::c3 = 0.21, + rodas34::c4 = 0.63, + rodas34::d1 = 0.25, + rodas34::d2 = -0.1043, + rodas34::d3 = 0.1035, + rodas34::d4 = -0.3620000000000023e-01, + rodas34::a21 = 0.1544e1, + rodas34::a31 = 0.9466785280815826, + rodas34::a32 = 0.2557011698983284, + rodas34::a41 = 0.3314825187068521e1, + rodas34::a42 = 0.2896124015972201e1, + rodas34::a43 = 0.9986419139977817, + rodas34::a51 = 0.1221224509226641e1, + rodas34::a52 = 0.6019134481288629e1, + rodas34::a53 = 0.1253708332932087e2, + rodas34::a54 = -0.6878860361058950, + rodas34::c21 = -0.56688e1, + rodas34::c31 = -0.2430093356833875e1, + rodas34::c32 = -0.2063599157091915, + rodas34::c41 = -0.1073529058151375, + rodas34::c42 = -0.9594562251023355e1, + rodas34::c43 = -0.2047028614809616e2, + rodas34::c51 = 0.7496443313967647e1, + rodas34::c52 = -0.1024680431464352e2, + rodas34::c53 = -0.3399990352819905e2, + rodas34::c54 = 0.1170890893206160e2, + rodas34::c61 = 0.8083246795921522e1, + rodas34::c62 = -0.7981132988064893e1, + rodas34::c63 = -0.3152159432874371e2, + rodas34::c64 = 0.1631930543123136e2, + rodas34::c65 = -0.6058818238834054e1, + rodas34::gamma = 0.25; } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict) +Foam::rodas34::rodas34(const ODESystem& ode, const dictionary& dict) : ODESolver(ode, dict), adaptiveSolver(ode, dict), @@ -93,9 +93,8 @@ Foam::rodas43::rodas43(const ODESystem& ode, const dictionary& dict) // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -Foam::scalar Foam::rodas43::solve +Foam::scalar Foam::rodas34::solve ( - const ODESystem& ode, const scalar x0, const scalarField& y0, const scalarField& dydx0, @@ -103,7 +102,7 @@ Foam::scalar Foam::rodas43::solve scalarField& y ) const { - ode.jacobian(x0, y0, dfdx_, dfdy_); + odes_.jacobian(x0, y0, dfdx_, dfdy_); for (register label i=0; i. Class - Foam::rodas43 + Foam::rodas34 Description L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (3)4. @@ -38,12 +38,12 @@ Description \endverbatim SourceFiles - rodas43.C + rodas34.C \*---------------------------------------------------------------------------*/ -#ifndef rodas43_H -#define rodas43_H +#ifndef rodas34_H +#define rodas34_H #include "ODESolver.H" #include "adaptiveSolver.H" @@ -54,10 +54,10 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Class rodas43 Declaration + Class rodas34 Declaration \*---------------------------------------------------------------------------*/ -class rodas43 +class rodas34 : public ODESolver, public adaptiveSolver @@ -92,13 +92,13 @@ class rodas43 public: //- Runtime type information - TypeName("rodas43"); + TypeName("rodas34"); // Constructors //- Construct from ODE - rodas43(const ODESystem& ode, const dictionary& dict); + rodas34(const ODESystem& ode, const dictionary& dict); // Member Functions @@ -106,7 +106,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, @@ -117,7 +116,6 @@ public: //- Solve the ODE system and the update the state void solve ( - const ODESystem& ode, scalar& x, scalarField& y, scalar& dxTry diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C index 69227f1963..cc788f9308 100644 --- a/src/ODE/ODESolvers/seulex/seulex.C +++ b/src/ODE/ODESolvers/seulex/seulex.C @@ -35,14 +35,13 @@ namespace Foam addToRunTimeSelectionTable(ODESolver, seulex, dictionary); const scalar - seulex::EPS = VSMALL, - seulex::STEPFAC1 = 0.6, - seulex::STEPFAC2 = 0.93, - seulex::STEPFAC3 = 0.1, - seulex::STEPFAC4 = 4.0, - seulex::STEPFAC5 = 0.5, - seulex::KFAC1 = 0.7, - seulex::KFAC2 = 0.9; + seulex::stepFactor1_ = 0.6, + seulex::stepFactor2_ = 0.93, + seulex::stepFactor3_ = 0.1, + seulex::stepFactor4_ = 4.0, + seulex::stepFactor5_ = 0.5, + seulex::kFactor1_ = 0.7, + seulex::kFactor2_ = 0.9; } @@ -51,150 +50,222 @@ namespace Foam Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict) : ODESolver(ode, dict), - ode_(ode), - nSeq_(IMAXX), - cost_(IMAXX), + jacRedo_(min(1.0e-4, min(relTol_))), + nSeq_(iMaxx_), + cpu_(iMaxx_), + coeff_(iMaxx_, iMaxx_), + theta_(2.0*jacRedo_), + table_(kMaxx_, n_), dfdx_(n_), - factrl_(IMAXX), - table_(KMAXX,n_), - fSave_((IMAXX-1)*(IMAXX+1)/2 + 2,n_), dfdy_(n_), - calcJac_(false), - dense_(false), a_(n_), - coeff_(IMAXX,IMAXX), - pivotIndices_(n_, 0.0), - hOpt_(IMAXX), - work_(IMAXX), - ySaved_(n_), + pivotIndices_(n_), + dxOpt_(iMaxx_), + temp_(iMaxx_), + y0_(n_), ySequence_(n_), scale_(n_), - del_(n_), + dy_(n_), yTemp_(n_), - dyTemp_(n_), - delTemp_(n_) + dydx_(n_) { - const scalar costfunc = 1.0, costjac = 5.0, costlu = 1.0, costsolve = 1.0; + // The CPU time factors for the major parts of the algorithm + const scalar cpuFunc = 1.0, cpuJac = 5.0, cpuLU = 1.0, cpuSolve = 1.0; - jacRedo_ = min(1.0e-4, min(relTol_)); - theta_ = 2.0*jacRedo_; nSeq_[0] = 2; nSeq_[1] = 3; - for (label i = 2; i < IMAXX; i++) + for (int i=2; i 1.0) + { + return false; + } + } + + odes_.derivatives(xnew, yTemp_, dy_); + LUBacksubstitute(a_, pivotIndices_, dy_); + } + + for (label i=0; i 0 ? true : false; - - ySaved_ = y; - - if (h != hnext_ && !first_step) + for (int j=k-1; j>0; j--) { - last_step = true; + for (label i=0; i jacRedo_ && !calcJac_) + if (theta_ > jacRedo_) { - ode_.jacobian(x, y, dfdx_, dfdy_); - calcJac_ = true; + odes_.jacobian(x, y, dfdx_, dfdy_); + jacUpdated = true; } - while (firstk || reject) + int k; + scalar dxNew = mag(dx); + bool firstk = true; + + while (firstk || step.reject) { - h = forward ? hnew : -hnew; + dx = step.forward ? dxNew : -dxNew; firstk = false; - reject = false; + step.reject = false; - if (mag(h) <= mag(x)*EPS) + if (mag(dx) <= mag(x)*SMALL) { - WarningIn("seulex::step(const scalar") - << "step size underflow in step :" << h << endl; + WarningIn("seulex::solve(scalar& x, scalarField& y, stepState&") + << "step size underflow :" << dx << endl; } - label ipt=-1; - for (k = 0; k <= kTarg_ + 1; k++) + scalar errOld; + + for (k=0; k<=kTarg_+1; k++) { - bool success = dy(x, ySaved_, h, k, ySequence_, ipt, scale_); + bool success = seul(x, y0_, dx, k, ySequence_, scale_); if (!success) { - reject = true; - hnew = mag(h)*STEPFAC5; + step.reject = true; + dxNew = mag(dx)*stepFactor5_; break; } @@ -204,45 +275,55 @@ void Foam::seulex::step(const scalar& htry, scalar& x, scalarField& y) const } else { - for (i = 0; i < n_; i++) + forAll(ySequence_, i) + { table_[k-1][i] = ySequence_[i]; + } } + if (k != 0) { - polyextr(k, table_, y); - err = 0.0; - for (i = 0; i < n_; i++) + extrapolate(k, table_, y); + scalar err = 0.0; + forAll(scale_, i) { - scale_[i] = absTol_[i] + relTol_[i]*mag(ySaved_[i]); + scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]); err += sqr((y[i] - table_[0][i])/scale_[i]); } err = sqrt(err/n_); - if (err > 1.0/EPS || (k > 1 && err >= errold)) + if (err > 1.0/SMALL || (k > 1 && err >= errOld)) { - reject = true; - hnew = mag(h)*STEPFAC5; + step.reject = true; + dxNew = mag(dx)*stepFactor5_; break; } - errold = min(4.0*err, 1.0); + errOld = min(4*err, 1); scalar expo = 1.0/(k + 1); - scalar facmin = pow(STEPFAC3, expo); + scalar facmin = pow(stepFactor3_, expo); + scalar fac; if (err == 0.0) { fac = 1.0/facmin; } else { - fac = STEPFAC2/pow(err/STEPFAC1, expo); - fac = max(facmin/STEPFAC4, min(1.0/facmin, fac)); + fac = stepFactor2_/pow(err/stepFactor1_, expo); + fac = max(facmin/stepFactor4_, min(1.0/facmin, fac)); } - hOpt_[k] = mag(h*fac); - work_[k] = cost_[k]/hOpt_[k]; + dxOpt_[k] = mag(dx*fac); + temp_[k] = cpu_[k]/dxOpt_[k]; - if ((first_step || last_step) && err <= 1.0) + if ((step.first || step.last) && err <= 1.0) { break; } - if (k == kTarg_-1 && !prev_reject && !first_step && !last_step) + + if + ( + k == kTarg_ - 1 + && !step.prevReject + && !step.first && !step.last + ) { if (err <= 1.0) { @@ -250,16 +331,17 @@ void Foam::seulex::step(const scalar& htry, scalar& x, scalarField& y) const } else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0) { - reject = true; + step.reject = true; kTarg_ = k; - if (kTarg_>1 && work_[k-1] < KFAC1*work_[k]) + if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k]) { kTarg_--; } - hnew = hOpt_[kTarg_]; + dxNew = dxOpt_[kTarg_]; break; } } + if (k == kTarg_) { if (err <= 1.0) @@ -268,51 +350,55 @@ void Foam::seulex::step(const scalar& htry, scalar& x, scalarField& y) const } else if (err > nSeq_[k + 1]*2.0) { - reject = true; - if (kTarg_>1 && work_[k-1] < KFAC1*work_[k]) + step.reject = true; + if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k]) { kTarg_--; } - hnew = hOpt_[kTarg_]; + dxNew = dxOpt_[kTarg_]; break; } } + if (k == kTarg_+1) { if (err > 1.0) { - reject = true; - if (kTarg_ > 1 && work_[kTarg_-1] < KFAC1*work_[kTarg_]) + step.reject = true; + if + ( + kTarg_ > 1 + && temp_[kTarg_-1] < kFactor1_*temp_[kTarg_] + ) { kTarg_--; } - hnew = hOpt_[kTarg_]; + dxNew = dxOpt_[kTarg_]; } break; } } } - if (reject) + if (step.reject) { - prev_reject = true; - if (!calcJac_) + step.prevReject = true; + if (!jacUpdated) { theta_ = 2.0*jacRedo_; - if (theta_ > jacRedo_ && !calcJac_) + if (theta_ > jacRedo_ && !jacUpdated) { - ode_.jacobian(x, y, dfdx_, dfdy_); - calcJac_ = true; + odes_.jacobian(x, y, dfdx_, dfdy_); + jacUpdated = true; } } } } - calcJac_ = false; + jacUpdated = false; - x += h; - hdid_ = h; - first_step = false; + step.dxDid = dx; + x += dx; label kopt; if (k == 1) @@ -322,180 +408,56 @@ void Foam::seulex::step(const scalar& htry, scalar& x, scalarField& y) const else if (k <= kTarg_) { kopt=k; - if (work_[k-1] < KFAC1*work_[k]) + if (temp_[k-1] < kFactor1_*temp_[k]) { kopt = k - 1; } - else if (work_[k] < KFAC2*work_[k - 1]) + else if (temp_[k] < kFactor2_*temp_[k - 1]) { - kopt = min(k + 1, KMAXX - 1); + kopt = min(k + 1, kMaxx_ - 1); } } else { kopt = k - 1; - if (k > 2 && work_[k-2] < KFAC1*work_[k - 1]) + if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1]) { kopt = k - 2; } - if (work_[k] < KFAC2*work_[kopt]) + if (temp_[k] < kFactor2_*temp_[kopt]) { - kopt = min(k, KMAXX - 1); + kopt = min(k, kMaxx_ - 1); } } - if (prev_reject) + if (step.prevReject) { kTarg_ = min(kopt, k); - hnew = min(mag(h), hOpt_[kTarg_]); - prev_reject = false; + dxNew = min(mag(dx), dxOpt_[kTarg_]); + step.prevReject = false; } else { if (kopt <= k) { - hnew = hOpt_[kopt]; + dxNew = dxOpt_[kopt]; } else { - if (k < kTarg_ && work_[k] < KFAC2*work_[k - 1]) + if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1]) { - hnew = hOpt_[k]*cost_[kopt + 1]/cost_[k]; + dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k]; } else { - hnew = hOpt_[k]*cost_[kopt]/cost_[k]; + dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k]; } } kTarg_ = kopt; } - if (forward) - { - hnext_ = hnew; - } - else - { - hnext_ =- hnew; - } + step.dxTry = step.forward ? dxNew : -dxNew; } -bool Foam::seulex::dy -( - const scalar& x, - scalarField& y, - const scalar htot, - const label k, - scalarField& yend, - label& ipt, - scalarField& scale -) const -{ - label nstep = nSeq_[k]; - scalar h = htot/nstep; - for (label i = 0; i < n_; i++) - { - for (label j = 0; j < n_; j++) a_[i][j] = -dfdy_[i][j]; - a_[i][i] += 1.0/h; - } - LUDecompose(a_, pivotIndices_); - scalar xnew = x + h; - - ode_.derivatives(xnew, y, del_); - - yTemp_ = y; - - LUBacksubstitute(a_, pivotIndices_, del_); - - if (dense_ && nstep == k + 1) - { - ipt++; - for (label i = 0; i < n_; i++) - { - fSave_[ipt][i] = del_[i]; - } - } - for (label nn = 1; nn < nstep; nn++) - { - for (label i=0;i 1.0) - { - return false; - } - } - - delTemp_ = yend; - LUBacksubstitute(a_, pivotIndices_, delTemp_); - del_ = delTemp_; - - if (dense_ && nn >= nstep-k-1) - { - ipt++; - for (label i=0;i 0; j--) - { - for (label i=0; i::solve cTp_[nSpecie] = T; cTp_[nSpecie+1] = p; - odeSolver_->solve - ( - *this, - 0, - deltaT, - cTp_, - subDeltaT - ); + odeSolver_->solve(0, deltaT, cTp_, subDeltaT); for (register int i=0; i