songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

seulex ODE solver: further rationalisation

Browse Source 
Careful handing of cached state and use of stepState class to carry information between sub-steps.

Also careful debugging of the handling of kTarg to resolve SEGV when running at low tolerance.
This commit is contained in:
Henry 2013-11-10 23:52:14 +00:00
parent e026a79e39
commit ee7284d4bd
2 changed files with 87 additions and 87 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

117
src/ODE/ODESolvers/seulex/seulex.C
View File

@ -50,16 +50,15 @@ namespace Foam
Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict) Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
: :
ODESolver(ode, dict), ODESolver(ode, dict),
jacRedo_(min(1.0e-4, min(relTol_))),
nSeq_(iMaxx_), nSeq_(iMaxx_),
cpu_(iMaxx_), cpu_(iMaxx_),
table_(kMaxx_, n_), coeff_(iMaxx_, iMaxx_),
jacRedo_(min(1.0e-4, min(relTol_))),
theta_(2.0*jacRedo_), theta_(2.0*jacRedo_),
calcJac_(false), table_(kMaxx_, n_),
dfdx_(n_), dfdx_(n_),
dfdy_(n_), dfdy_(n_),
a_(n_), a_(n_),
coeff_(iMaxx_, iMaxx_),
pivotIndices_(n_), pivotIndices_(n_),
dxOpt_(iMaxx_), dxOpt_(iMaxx_),
temp_(iMaxx_), temp_(iMaxx_),
@ -70,6 +69,7 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
yTemp_(n_), yTemp_(n_),
dydx_(n_) dydx_(n_)
{ {
// 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; const scalar cpuFunc = 1.0, cpuJac = 5.0, cpuLU = 1.0, cpuSolve = 1.0;
nSeq_[0] = 2; nSeq_[0] = 2;
@ -86,10 +86,7 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict)
cpu_[k+1] = cpu_[k] + (nSeq_[k+1]-1)*(cpuFunc + cpuSolve) + cpuLU; cpu_[k+1] = cpu_[k] + (nSeq_[k+1]-1)*(cpuFunc + cpuSolve) + cpuLU;
} }
//NOTE: the first element of relTol_ and absTol_ are used here. // Set the extrapolation coefficients array
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++) for (int k=0; k<iMaxx_; k++)
{ {
for (int l=0; l<k; l++) for (int l=0; l<k; l++)
@ -142,7 +139,7 @@ bool Foam::seulex::seul
if (nn == 1 && k<=1) if (nn == 1 && k<=1)
{ {
scalar dy1 = 0.0; scalar dy1 = 0.0;
for (label i = 0; i < n_; i++) for (label i=0; i<n_; i++)
{ {
dy1 += sqr(dy_[i]/scale[i]); dy1 += sqr(dy_[i]/scale[i]);
} }
@ -157,7 +154,7 @@ bool Foam::seulex::seul
LUBacksubstitute(a_, pivotIndices_, dy_); LUBacksubstitute(a_, pivotIndices_, dy_);
scalar dy2 = 0.0; scalar dy2 = 0.0;
for (label i = 0; i <n_ ; i++) for (label i=0; i<n_; i++)
{ {
dy2 += sqr(dy_[i]/scale[i]); dy2 += sqr(dy_[i]/scale[i]);
} }
@ -210,30 +207,24 @@ void Foam::seulex::solve
( (
scalar& x, scalar& x,
scalarField& y, scalarField& y,
scalar& dxTry stepState& step
) const ) const
{ {
static bool firstStep = true, lastStep = false;
static bool reject = false, prevReject = false;
temp_[0] = GREAT; temp_[0] = GREAT;
scalar dx = dxTry; scalar dx = step.dxTry;
dxTry = -VGREAT;
bool forward = dx > 0 ? true : false;
y0_ = y; y0_ = y;
dxOpt_[0] = mag(0.1*dx);
if (dx != dxTry && !firstStep) if (step.first || step.prevReject)
{ {
lastStep = true; theta_ = 2.0*jacRedo_;
} }
if (reject) if (step.first)
{ {
prevReject = true; // NOTE: the first element of relTol_ and absTol_ are used here.
lastStep = false; scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
theta_=2.0*jacRedo_; kTarg_ = max(1, min(kMaxx_ - 1, int(logTol)));
} }
forAll(scale_, i) forAll(scale_, i)
@ -241,39 +232,39 @@ void Foam::seulex::solve
scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]); scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
} }
reject = false; bool jacUpdated = false;
bool firstk = true;
scalar dxNew = mag(dx);
if (theta_ > jacRedo_ && !calcJac_) if (theta_ > jacRedo_)
{ {
odes_.jacobian(x, y, dfdx_, dfdy_); odes_.jacobian(x, y, dfdx_, dfdy_);
calcJac_ = true; jacUpdated = true;
} }
int k; int k;
scalar errOld; scalar dxNew = mag(dx);
bool firstk = true;
while (firstk || reject) while (firstk || step.reject)
{ {
dx = forward ? dxNew : -dxNew; dx = step.forward ? dxNew : -dxNew;
firstk = false; firstk = false;
reject = false; step.reject = false;
if (mag(dx) <= mag(x)*SMALL) if (mag(dx) <= mag(x)*SMALL)
{ {
WarningIn("seulex::step(const scalar") WarningIn("seulex::solve(scalar& x, scalarField& y, stepState&")
<< "step size underflow in step :" << dx << endl; << "step size underflow :" << dx << endl;
} }
for (k = 0; k <= kTarg_ + 1; k++) scalar errOld;
for (k=0; k<=kTarg_+1; k++)
{ {
bool success = seul(x, y0_, dx, k, ySequence_, scale_); bool success = seul(x, y0_, dx, k, ySequence_, scale_);
if (!success) if (!success)
{ {
reject = true; step.reject = true;
dxNew = mag(dx)*stepFactor5_; dxNew = mag(dx)*stepFactor5_;
break; break;
} }
@ -289,6 +280,7 @@ void Foam::seulex::solve
table_[k-1][i] = ySequence_[i]; table_[k-1][i] = ySequence_[i];
} }
} }
if (k != 0) if (k != 0)
{ {
extrapolate(k, table_, y); extrapolate(k, table_, y);
@ -301,7 +293,7 @@ void Foam::seulex::solve
err = sqrt(err/n_); err = sqrt(err/n_);
if (err > 1.0/SMALL || (k > 1 && err >= errOld)) if (err > 1.0/SMALL || (k > 1 && err >= errOld))
{ {
reject = true; step.reject = true;
dxNew = mag(dx)*stepFactor5_; dxNew = mag(dx)*stepFactor5_;
break; break;
} }
@ -321,11 +313,17 @@ void Foam::seulex::solve
dxOpt_[k] = mag(dx*fac); dxOpt_[k] = mag(dx*fac);
temp_[k] = cpu_[k]/dxOpt_[k]; temp_[k] = cpu_[k]/dxOpt_[k];
if ((firstStep || lastStep) && err <= 1.0) if ((step.first || step.last) && err <= 1.0)
{ {
break; break;
} }
if (k == kTarg_-1 && !prevReject && !firstStep && !lastStep)
if
(
k == kTarg_ - 1
&& !step.prevReject
&& !step.first && !step.last
)
{ {
if (err <= 1.0) if (err <= 1.0)
{ {
@ -333,7 +331,7 @@ void Foam::seulex::solve
} }
else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0) else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0)
{ {
reject = true; step.reject = true;
kTarg_ = k; kTarg_ = k;
if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k]) if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
{ {
@ -343,6 +341,7 @@ void Foam::seulex::solve
break; break;
} }
} }
if (k == kTarg_) if (k == kTarg_)
{ {
if (err <= 1.0) if (err <= 1.0)
@ -351,7 +350,7 @@ void Foam::seulex::solve
} }
else if (err > nSeq_[k + 1]*2.0) else if (err > nSeq_[k + 1]*2.0)
{ {
reject = true; step.reject = true;
if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k]) if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
{ {
kTarg_--; kTarg_--;
@ -360,11 +359,12 @@ void Foam::seulex::solve
break; break;
} }
} }
if (k == kTarg_+1) if (k == kTarg_+1)
{ {
if (err > 1.0) if (err > 1.0)
{ {
reject = true; step.reject = true;
if if
( (
kTarg_ > 1 kTarg_ > 1
@ -379,26 +379,26 @@ void Foam::seulex::solve
} }
} }
} }
if (reject) if (step.reject)
{ {
prevReject = true; step.prevReject = true;
if (!calcJac_) if (!jacUpdated)
{ {
theta_ = 2.0*jacRedo_; theta_ = 2.0*jacRedo_;
if (theta_ > jacRedo_ && !calcJac_) if (theta_ > jacRedo_ && !jacUpdated)
{ {
odes_.jacobian(x, y, dfdx_, dfdy_); odes_.jacobian(x, y, dfdx_, dfdy_);
calcJac_ = true; jacUpdated = true;
} }
} }
} }
} }
calcJac_ = false; jacUpdated = false;
step.dxDid = dx;
x += dx; x += dx;
firstStep = false;
label kopt; label kopt;
if (k == 1) if (k == 1)
@ -430,11 +430,11 @@ void Foam::seulex::solve
} }
} }
if (prevReject) if (step.prevReject)
{ {
kTarg_ = min(kopt, k); kTarg_ = min(kopt, k);
dxNew = min(mag(dx), dxOpt_[kTarg_]); dxNew = min(mag(dx), dxOpt_[kTarg_]);
prevReject = false; step.prevReject = false;
} }
else else
{ {
@ -456,14 +456,7 @@ void Foam::seulex::solve
kTarg_ = kopt; kTarg_ = kopt;
} }
if (forward) step.dxTry = step.forward ? dxNew : -dxNew;
{
dxTry = dxNew;
}
else
{
dxTry = -dxNew;
}
} }

57
src/ODE/ODESolvers/seulex/seulex.H
View File

@ -28,7 +28,7 @@ Description
An extrapolation-algorithm, based on the linearly implicit Euler method An extrapolation-algorithm, based on the linearly implicit Euler method
with step size control and order selection. with step size control and order selection.
The implementation is based on the SEULEX code in The implementation is based on the SEULEX algorithm in
\verbatim \verbatim
"Solving Ordinary Differential Equations II: Stiff "Solving Ordinary Differential Equations II: Stiff
and Differential-Algebraic Problems, second edition", and Differential-Algebraic Problems, second edition",
@ -65,35 +65,42 @@ class seulex
{ {
// Private data // Private data
// SEULEX constants // Static constants
static const label kMaxx_ = 12;
static const label iMaxx_ = kMaxx_ + 1;
static const scalar static const label kMaxx_ = 12;
stepFactor1_, stepFactor2_, stepFactor3_, static const label iMaxx_ = kMaxx_ + 1;
stepFactor4_, stepFactor5_,
kFactor1_, kFactor2_;
// Temporary fields static const scalar
mutable label kTarg_; stepFactor1_, stepFactor2_, stepFactor3_,
mutable labelField nSeq_; stepFactor4_, stepFactor5_,
mutable scalarField cpu_; kFactor1_, kFactor2_;
mutable scalarRectangularMatrix table_;
mutable scalar jacRedo_, theta_; // Evaluated constants
mutable bool calcJac_;
mutable scalarField dfdx_; scalar jacRedo_;
mutable scalarSquareMatrix dfdy_; labelField nSeq_;
mutable scalarSquareMatrix a_, coeff_; scalarField cpu_;
mutable labelList pivotIndices_; scalarSquareMatrix coeff_;
// Fields space for "solve" function // Temporary storage
mutable scalarField dxOpt_, temp_; // held to avoid dynamic memory allocation between calls
mutable scalarField y0_, ySequence_, scale_; // and to transfer internal values between functions
// Fields used in "dy" function mutable scalar theta_;
mutable scalarField dy_, yTemp_, dydx_; mutable label kTarg_;
mutable scalarRectangularMatrix table_;
mutable scalarField dfdx_;
mutable scalarSquareMatrix dfdy_;
mutable scalarSquareMatrix a_;
mutable labelList pivotIndices_;
// Fields space for "solve" function
mutable scalarField dxOpt_, temp_;
mutable scalarField y0_, ySequence_, scale_;
// Fields used in "seul" function
mutable scalarField dy_, yTemp_, dydx_;
// Private Member Functions // Private Member Functions
@ -137,7 +144,7 @@ public:
( (
scalar& x, scalar& x,
scalarField& y, scalarField& y,
scalar& dxTry stepState& step
) const; ) const;
}; };