diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C
index 5b510c723c79decf995824c7bb5efc691bd8ae09..19abaf6f9bc661274733378b2537c6a5d1e2c8c9 100644
--- a/src/ODE/ODESolvers/seulex/seulex.C
+++ b/src/ODE/ODESolvers/seulex/seulex.C
@@ -51,31 +51,27 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
:
ODESolver(ode, dict),
nSeq_(iMaxx_),
- cost_(iMaxx_),
- dfdx_(n_),
- factrl_(iMaxx_),
+ cpu_(iMaxx_),
table_(kMaxx_, n_),
- fSave_((iMaxx_ - 1)*(iMaxx_ + 1)/2 + 2, n_),
- dfdy_(n_),
+ jacRedo_(min(1.0e-4, min(relTol_))),
+ theta_(2.0*jacRedo_),
calcJac_(false),
- dense_(false),
+ dfdx_(n_),
+ dfdy_(n_),
a_(n_),
coeff_(iMaxx_, iMaxx_),
- pivotIndices_(n_, 0.0),
+ pivotIndices_(n_),
dxOpt_(iMaxx_),
- work_(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;
+ 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;
@@ -83,17 +79,16 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
{
nSeq_[i] = 2*nSeq_[i-2];
}
- cost_[0] = costjac + costlu + nSeq_[0]*(costfunc + costsolve);
+ cpu_[0] = cpuJac + cpuLU + nSeq_[0]*(cpuFunc + cpuSolve);
for (int k=0; k<kMaxx_; k++)
{
- cost_[k+1] = cost_[k] + (nSeq_[k+1]-1)*(costfunc + costsolve) + costlu;
+ cpu_[k+1] = cpu_[k] + (nSeq_[k+1]-1)*(cpuFunc + cpuSolve) + cpuLU;
}
//NOTE: the first element of relTol_ and absTol_ are used here.
- scalar logfact =- log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
-
- kTarg_ = min(1, min(kMaxx_ - 1, label(logfact)));
+ scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
+ kTarg_ = min(1, min(kMaxx_ - 1, int(logTol)));
for (int k=0; k<iMaxx_; k++)
{
@@ -103,30 +98,23 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
coeff_[k][l] = 1.0/(ratio - 1.0);
}
}
-
- factrl_[0] = 1.0;
- for (int k=0; k<iMaxx_ - 1; k++)
- {
- factrl_[k+1] = (k+1)*factrl_[k];
- }
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-bool Foam::seulex::dy
+bool Foam::seulex::seul
(
- const scalar x,
- scalarField& y,
+ const scalar x0,
+ const scalarField& y0,
const scalar dxTot,
const label k,
- scalarField& yend,
- label& ipt,
- scalarField& scale
+ scalarField& y,
+ const scalarField& scale
) const
{
- label nstep = nSeq_[k];
- scalar dx = dxTot/nstep;
+ label nSteps = nSeq_[k];
+ scalar dx = dxTot/nSteps;
for (label i=0; i<n_; i++)
{
@@ -140,56 +128,41 @@ bool Foam::seulex::dy
LUDecompose(a_, pivotIndices_);
- scalar xnew = x + dx;
-
- odes_.derivatives(xnew, y, del_);
+ scalar xnew = x0 + dx;
+ odes_.derivatives(xnew, y0, dy_);
+ LUBacksubstitute(a_, pivotIndices_, dy_);
- yTemp_ = y;
+ yTemp_ = y0;
- 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 nn=1; nn<nSteps; nn++)
{
- for (label i=0;i<n_;i++)
- {
- yTemp_[i] += del_[i];
- }
+ yTemp_ += dy_;
xnew += dx;
- odes_.derivatives(xnew, yTemp_, yend);
-
if (nn == 1 && k<=1)
{
- scalar del1=0.0;
+ scalar dy1 = 0.0;
for (label i = 0; i < n_; i++)
{
- del1 += sqr(del_[i]/scale[i]);
+ dy1 += sqr(dy_[i]/scale[i]);
}
- del1 = sqrt(del1);
+ dy1 = sqrt(dy1);
- odes_.derivatives(x + dx, yTemp_, dyTemp_);
- for (label i=0;i<n_;i++)
+ odes_.derivatives(x0 + dx, yTemp_, dydx_);
+ for (label i=0; i<n_; i++)
{
- del_[i] = dyTemp_[i] - del_[i]/dx;
+ dy_[i] = dydx_[i] - dy_[i]/dx;
}
- LUBacksubstitute(a_, pivotIndices_, del_);
+ LUBacksubstitute(a_, pivotIndices_, dy_);
- scalar del2 = 0.0;
+ scalar dy2 = 0.0;
for (label i = 0; i <n_ ; i++)
{
- del2 += sqr(del_[i]/scale[i]);
+ dy2 += sqr(dy_[i]/scale[i]);
}
- del2 = sqrt(del2);
- theta_ = del2 / min(1.0, del1 + SMALL);
+ dy2 = sqrt(dy2);
+ theta_ = dy2/min(1.0, dy1 + SMALL);
if (theta_ > 1.0)
{
@@ -197,46 +170,38 @@ bool Foam::seulex::dy
}
}
- delTemp_ = yend;
- LUBacksubstitute(a_, pivotIndices_, delTemp_);
- del_ = delTemp_;
-
- if (dense_ && nn >= nstep-k-1)
- {
- ipt++;
- for (label i=0;i<n_;i++)
- {
- fSave_[ipt][i]=del_[i];
- }
- }
+ odes_.derivatives(xnew, yTemp_, dy_);
+ LUBacksubstitute(a_, pivotIndices_, dy_);
}
- yend = yTemp_ + del_;
+ for (label i=0; i<n_; i++)
+ {
+ y[i] = yTemp_[i] + dy_[i];
+ }
return true;
}
-void Foam::seulex::polyextr
+void Foam::seulex::extrapolate
(
const label k,
scalarRectangularMatrix& table,
- scalarField& last
+ scalarField& y
) const
{
- int l = last.size();
- for (int j=k - 1; j>0; j--)
+ for (int j=k-1; j>0; j--)
{
- for (label i=0; i<l; i++)
+ for (label i=0; i<n_; i++)
{
table[j-1][i] =
table[j][i] + coeff_[k][j]*(table[j][i] - table[j-1][i]);
}
}
- for (int i=0; i<l; i++)
+ for (int i=0; i<n_; i++)
{
- last[i] = table[0][i] + coeff_[k][0]*(table[0][i] - last[i]);
+ y[i] = table[0][i] + coeff_[k][0]*(table[0][i] - y[i]);
}
}
@@ -250,9 +215,8 @@ void Foam::seulex::solve
{
static bool firstStep = true, lastStep = false;
static bool reject = false, prevReject = false;
- static scalar errold;
- work_[0] = GREAT;
+ temp_[0] = GREAT;
scalar dx = dxTry;
dxTry = -VGREAT;
@@ -289,6 +253,8 @@ void Foam::seulex::solve
int k;
+ scalar errOld;
+
while (firstk || reject)
{
dx = forward ? dxNew : -dxNew;
@@ -300,11 +266,10 @@ void Foam::seulex::solve
WarningIn("seulex::step(const scalar")
<< "step size underflow in step :" << dx << endl;
}
- label ipt = -1;
for (k = 0; k <= kTarg_ + 1; k++)
{
- bool success = dy(x, y0_, dx, k, ySequence_, ipt, scale_);
+ bool success = seul(x, y0_, dx, k, ySequence_, scale_);
if (!success)
{
@@ -326,7 +291,7 @@ void Foam::seulex::solve
}
if (k != 0)
{
- polyextr(k, table_, y);
+ extrapolate(k, table_, y);
scalar err = 0.0;
forAll(scale_, i)
{
@@ -334,13 +299,13 @@ void Foam::seulex::solve
err += sqr((y[i] - table_[0][i])/scale_[i]);
}
err = sqrt(err/n_);
- if (err > 1.0/SMALL || (k > 1 && err >= errold))
+ if (err > 1.0/SMALL || (k > 1 && err >= errOld))
{
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(stepFactor3_, expo);
scalar fac;
@@ -354,7 +319,7 @@ void Foam::seulex::solve
fac = max(facmin/stepFactor4_, min(1.0/facmin, fac));
}
dxOpt_[k] = mag(dx*fac);
- work_[k] = cost_[k]/dxOpt_[k];
+ temp_[k] = cpu_[k]/dxOpt_[k];
if ((firstStep || lastStep) && err <= 1.0)
{
@@ -370,7 +335,7 @@ void Foam::seulex::solve
{
reject = true;
kTarg_ = k;
- if (kTarg_>1 && work_[k-1] < kFactor1_*work_[k])
+ if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
{
kTarg_--;
}
@@ -387,7 +352,7 @@ void Foam::seulex::solve
else if (err > nSeq_[k + 1]*2.0)
{
reject = true;
- if (kTarg_>1 && work_[k-1] < kFactor1_*work_[k])
+ if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
{
kTarg_--;
}
@@ -403,7 +368,7 @@ void Foam::seulex::solve
if
(
kTarg_ > 1
- && work_[kTarg_-1] < kFactor1_*work_[kTarg_]
+ && temp_[kTarg_-1] < kFactor1_*temp_[kTarg_]
)
{
kTarg_--;
@@ -443,11 +408,11 @@ void Foam::seulex::solve
else if (k <= kTarg_)
{
kopt=k;
- if (work_[k-1] < kFactor1_*work_[k])
+ if (temp_[k-1] < kFactor1_*temp_[k])
{
kopt = k - 1;
}
- else if (work_[k] < kFactor2_*work_[k - 1])
+ else if (temp_[k] < kFactor2_*temp_[k - 1])
{
kopt = min(k + 1, kMaxx_ - 1);
}
@@ -455,11 +420,11 @@ void Foam::seulex::solve
else
{
kopt = k - 1;
- if (k > 2 && work_[k-2] < kFactor1_*work_[k - 1])
+ if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1])
{
kopt = k - 2;
}
- if (work_[k] < kFactor2_*work_[kopt])
+ if (temp_[k] < kFactor2_*temp_[kopt])
{
kopt = min(k, kMaxx_ - 1);
}
@@ -479,13 +444,13 @@ void Foam::seulex::solve
}
else
{
- if (k < kTarg_ && work_[k] < kFactor2_*work_[k - 1])
+ if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1])
{
- dxNew = dxOpt_[k]*cost_[kopt + 1]/cost_[k];
+ dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k];
}
else
{
- dxNew = dxOpt_[k]*cost_[kopt]/cost_[k];
+ dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k];
}
}
kTarg_ = kopt;
diff --git a/src/ODE/ODESolvers/seulex/seulex.H b/src/ODE/ODESolvers/seulex/seulex.H
index a94054cb474cb534191307800b52a03dbbf666e8..8c01052331b1c1853215edf5e68fef8853d58572 100644
--- a/src/ODE/ODESolvers/seulex/seulex.H
+++ b/src/ODE/ODESolvers/seulex/seulex.H
@@ -77,40 +77,44 @@ class seulex
// Temporary fields
mutable label kTarg_;
mutable labelField nSeq_;
- mutable scalarField cost_, dfdx_, factrl_;
- mutable scalarRectangularMatrix table_, fSave_;
- mutable scalarSquareMatrix dfdy_;
+ mutable scalarField cpu_;
+ mutable scalarRectangularMatrix table_;
+
mutable scalar jacRedo_, theta_;
- mutable bool calcJac_, dense_;
+ mutable bool calcJac_;
+
+ mutable scalarField dfdx_;
+ mutable scalarSquareMatrix dfdy_;
mutable scalarSquareMatrix a_, coeff_;
mutable labelList pivotIndices_;
// Fields space for "solve" function
- mutable scalarField dxOpt_, work_, y0_, ySequence_, scale_;
+ mutable scalarField dxOpt_, temp_;
+ mutable scalarField y0_, ySequence_, scale_;
// Fields used in "dy" function
- mutable scalarField del_, yTemp_, dyTemp_, delTemp_;
+ mutable scalarField dy_, yTemp_, dydx_;
// Private Member Functions
- bool dy
+ //- Computes the j-th line of the extrapolation table
+ bool seul
(
- const scalar x,
- scalarField& y,
+ const scalar x0,
+ const scalarField& y0,
const scalar dxTot,
const label k,
- scalarField& yend,
- label& ipt,
- scalarField& scale
+ scalarField& y,
+ const scalarField& scale
) const;
-
- void polyextr
+ //- Polynomial extrpolation
+ void extrapolate
(
const label k,
scalarRectangularMatrix& table,
- scalarField& last
+ scalarField& y
) const;
@@ -127,12 +131,14 @@ public:
// Member Functions
- void solve
- (
- scalar& x,
- scalarField& y,
- scalar& dxTry
- ) const;
+
+ //- Solve the ODE system and the update the state
+ void solve
+ (
+ scalar& x,
+ scalarField& y,
+ scalar& dxTry
+ ) const;
};
diff --git a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
index afb86830866626eeac5c944184a5436569217113..3ca485902ed18b062415a5a63b711218b657dcb5 100644
--- a/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
+++ b/src/thermophysicalModels/chemistryModel/chemistrySolver/ode/ode.C
@@ -70,14 +70,7 @@ void Foam::ode<ChemistryModel>::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<nSpecie; i++)
{