diff --git a/src/ODE/Make/files b/src/ODE/Make/files
index a22ec1d3be..34ef0d4cb1 100644
--- a/src/ODE/Make/files
+++ b/src/ODE/Make/files
@@ -14,5 +14,6 @@ ODESolvers/rodas43/rodas43.C
ODESolvers/SIBS/SIBS.C
ODESolvers/SIBS/SIMPR.C
ODESolvers/SIBS/polyExtrapolate.C
+ODESolvers/seulex/seulex.C
LIB = $(FOAM_LIBBIN)/libODE
diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C
new file mode 100644
index 0000000000..69227f1963
--- /dev/null
+++ b/src/ODE/ODESolvers/seulex/seulex.C
@@ -0,0 +1,501 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see .
+
+\*---------------------------------------------------------------------------*/
+
+#include "seulex.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(seulex, 0);
+
+ 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;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
+:
+ ODESolver(ode, dict),
+ ode_(ode),
+ nSeq_(IMAXX),
+ cost_(IMAXX),
+ 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_),
+ ySequence_(n_),
+ scale_(n_),
+ del_(n_),
+ yTemp_(n_),
+ dyTemp_(n_),
+ delTemp_(n_)
+{
+ const scalar costfunc = 1.0, costjac = 5.0, costlu = 1.0, costsolve = 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++)
+ {
+ nSeq_[i] = 2*nSeq_[i-2];
+ }
+ cost_[0] = costjac + costlu + nSeq_[0]*(costfunc + costsolve);
+
+ for (label k=0;k 0 ? true : false;
+
+ ySaved_ = y;
+
+ if (h != hnext_ && !first_step)
+ {
+ last_step = true;
+ }
+
+ if (reject)
+ {
+ prev_reject = true;
+ last_step = false;
+ theta_=2.0*jacRedo_;
+ }
+
+ for (i = 0; i < n_; i++)
+ {
+ scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
+ }
+
+ reject = false;
+ firstk = true;
+ hnew = mag(h);
+
+ if (theta_ > jacRedo_ && !calcJac_)
+ {
+ ode_.jacobian(x, y, dfdx_, dfdy_);
+ calcJac_ = true;
+ }
+
+ while (firstk || reject)
+ {
+ h = forward ? hnew : -hnew;
+ firstk = false;
+ reject = false;
+
+ if (mag(h) <= mag(x)*EPS)
+ {
+ WarningIn("seulex::step(const scalar")
+ << "step size underflow in step :" << h << endl;
+ }
+ label ipt=-1;
+
+ for (k = 0; k <= kTarg_ + 1; k++)
+ {
+ bool success = dy(x, ySaved_, h, k, ySequence_, ipt, scale_);
+
+ if (!success)
+ {
+ reject = true;
+ hnew = mag(h)*STEPFAC5;
+ break;
+ }
+
+ if (k == 0)
+ {
+ y = ySequence_;
+ }
+ else
+ {
+ for (i = 0; i < n_; i++)
+ table_[k-1][i] = ySequence_[i];
+ }
+ if (k != 0)
+ {
+ polyextr(k, table_, y);
+ err = 0.0;
+ for (i = 0; i < n_; i++)
+ {
+ scale_[i] = absTol_[i] + relTol_[i]*mag(ySaved_[i]);
+ err += sqr((y[i] - table_[0][i])/scale_[i]);
+ }
+ err = sqrt(err/n_);
+ if (err > 1.0/EPS || (k > 1 && err >= errold))
+ {
+ reject = true;
+ hnew = mag(h)*STEPFAC5;
+ break;
+ }
+ errold = min(4.0*err, 1.0);
+ scalar expo = 1.0/(k + 1);
+ scalar facmin = pow(STEPFAC3, expo);
+ if (err == 0.0)
+ {
+ fac = 1.0/facmin;
+ }
+ else
+ {
+ fac = STEPFAC2/pow(err/STEPFAC1, expo);
+ fac = max(facmin/STEPFAC4, min(1.0/facmin, fac));
+ }
+ hOpt_[k] = mag(h*fac);
+ work_[k] = cost_[k]/hOpt_[k];
+
+ if ((first_step || last_step) && err <= 1.0)
+ {
+ break;
+ }
+ if (k == kTarg_-1 && !prev_reject && !first_step && !last_step)
+ {
+ if (err <= 1.0)
+ {
+ break;
+ }
+ else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0)
+ {
+ reject = true;
+ kTarg_ = k;
+ if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
+ {
+ kTarg_--;
+ }
+ hnew = hOpt_[kTarg_];
+ break;
+ }
+ }
+ if (k == kTarg_)
+ {
+ if (err <= 1.0)
+ {
+ break;
+ }
+ else if (err > nSeq_[k + 1]*2.0)
+ {
+ reject = true;
+ if (kTarg_>1 && work_[k-1] < KFAC1*work_[k])
+ {
+ kTarg_--;
+ }
+ hnew = hOpt_[kTarg_];
+ break;
+ }
+ }
+ if (k == kTarg_+1)
+ {
+ if (err > 1.0)
+ {
+ reject = true;
+ if (kTarg_ > 1 && work_[kTarg_-1] < KFAC1*work_[kTarg_])
+ {
+ kTarg_--;
+ }
+ hnew = hOpt_[kTarg_];
+ }
+ break;
+ }
+ }
+ }
+ if (reject)
+ {
+ prev_reject = true;
+ if (!calcJac_)
+ {
+ theta_ = 2.0*jacRedo_;
+
+ if (theta_ > jacRedo_ && !calcJac_)
+ {
+ ode_.jacobian(x, y, dfdx_, dfdy_);
+ calcJac_ = true;
+ }
+ }
+ }
+ }
+
+ calcJac_ = false;
+
+ x += h;
+ hdid_ = h;
+ first_step = false;
+
+ label kopt;
+ if (k == 1)
+ {
+ kopt = 2;
+ }
+ else if (k <= kTarg_)
+ {
+ kopt=k;
+ if (work_[k-1] < KFAC1*work_[k])
+ {
+ kopt = k - 1;
+ }
+ else if (work_[k] < KFAC2*work_[k - 1])
+ {
+ kopt = min(k + 1, KMAXX - 1);
+ }
+ }
+ else
+ {
+ kopt = k - 1;
+ if (k > 2 && work_[k-2] < KFAC1*work_[k - 1])
+ {
+ kopt = k - 2;
+ }
+ if (work_[k] < KFAC2*work_[kopt])
+ {
+ kopt = min(k, KMAXX - 1);
+ }
+ }
+
+ if (prev_reject)
+ {
+ kTarg_ = min(kopt, k);
+ hnew = min(mag(h), hOpt_[kTarg_]);
+ prev_reject = false;
+ }
+ else
+ {
+ if (kopt <= k)
+ {
+ hnew = hOpt_[kopt];
+ }
+ else
+ {
+ if (k < kTarg_ && work_[k] < KFAC2*work_[k - 1])
+ {
+ hnew = hOpt_[k]*cost_[kopt + 1]/cost_[k];
+ }
+ else
+ {
+ hnew = hOpt_[k]*cost_[kopt]/cost_[k];
+ }
+ }
+ kTarg_ = kopt;
+ }
+
+ if (forward)
+ {
+ hnext_ = hnew;
+ }
+ else
+ {
+ hnext_ =- hnew;
+ }
+}
+
+
+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.
+
+Class
+ Foam::seulex
+
+Description
+ An extrapolation-algorithm, based on the linearly implicit Euler method
+ with step size control and order selection.
+
+ The implementation is based on the SEULEX code in
+ \verbatim
+ "Solving Ordinary Differential Equations II: Stiff
+ and Differential-Algebraic Problems, second edition",
+ Hairer, E.,
+ Nørsett, S.,
+ Wanner, G.,
+ Springer-Verlag, Berlin. 1996.
+ \endverbatim
+
+SourceFiles
+ seulex.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef seulex_H
+#define seulex_H
+
+#include "ODESolver.H"
+
+#include "scalarFieldField.H"
+#include "scalarMatrices.H"
+#include "labelField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class seulex Declaration
+\*---------------------------------------------------------------------------*/
+
+class seulex
+:
+ public ODESolver
+{
+
+ // Private data
+
+ //- seulex constants
+ static const scalar EPS;
+
+ static const label KMAXX=12, IMAXX = KMAXX + 1;
+
+ static const scalar
+ STEPFAC1, STEPFAC2, STEPFAC3, STEPFAC4, STEPFAC5, KFAC1, KFAC2;
+
+
+ //- Reference to ODESystem
+ const ODESystem& ode_;
+
+ // Temporary fields
+ mutable label kTarg_;
+ mutable labelField nSeq_;
+ mutable scalarField cost_, dfdx_, factrl_;
+ mutable scalarRectangularMatrix table_, fSave_;
+ mutable scalarSquareMatrix dfdy_;
+ mutable scalar jacRedo_, theta_;
+ mutable bool calcJac_, dense_;
+ mutable scalarSquareMatrix a_, coeff_;
+
+
+ mutable scalar hnext_, hdid_;
+
+ mutable labelList pivotIndices_;
+
+ // Fields space for "step" function
+ mutable scalarField hOpt_ ,work_, ySaved_, ySequence_, scale_;
+
+ // Fields used in "dy" function
+ mutable scalarField del_, yTemp_, dyTemp_, delTemp_;
+
+
+ // Private Member Functions
+
+ void step(const scalar& htry, scalar& x, scalarField& y) const;
+
+ bool dy
+ (
+ const scalar& x,
+ scalarField& y,
+ const scalar htot,
+ const label k,
+ scalarField& yend,
+ label& ipt,
+ scalarField& scale
+ ) const;
+
+
+ void polyextr
+ (
+ const label k,
+ scalarRectangularMatrix& table,
+ scalarField& last
+ ) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("seulex");
+
+
+ // Constructors
+
+ //- Construct from ODE
+ seulex(const ODESystem& ode, const dictionary& dict);
+
+
+ // Member Functions
+ void solve
+ (
+ const ODESystem& ode,
+ scalar& x,
+ scalarField& y,
+ scalar& dxTry
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //