From ee7284d4bdf1aae61a744ff46b9fb416e7149034 Mon Sep 17 00:00:00 2001 From: Henry Date: Sun, 10 Nov 2013 23:52:14 +0000 Subject: [PATCH] seulex ODE solver: further rationalisation 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. --- src/ODE/ODESolvers/seulex/seulex.C | 117 ++++++++++++++--------------- src/ODE/ODESolvers/seulex/seulex.H | 57 ++++++++------ 2 files changed, 87 insertions(+), 87 deletions(-) diff --git a/src/ODE/ODESolvers/seulex/seulex.C b/src/ODE/ODESolvers/seulex/seulex.C index 19abaf6f9b..cc788f9308 100644 --- a/src/ODE/ODESolvers/seulex/seulex.C +++ b/src/ODE/ODESolvers/seulex/seulex.C @@ -50,16 +50,15 @@ namespace Foam Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict) : ODESolver(ode, dict), + jacRedo_(min(1.0e-4, min(relTol_))), nSeq_(iMaxx_), cpu_(iMaxx_), - table_(kMaxx_, n_), - jacRedo_(min(1.0e-4, min(relTol_))), + coeff_(iMaxx_, iMaxx_), theta_(2.0*jacRedo_), - calcJac_(false), + table_(kMaxx_, n_), dfdx_(n_), dfdy_(n_), a_(n_), - coeff_(iMaxx_, iMaxx_), pivotIndices_(n_), dxOpt_(iMaxx_), temp_(iMaxx_), @@ -70,6 +69,7 @@ Foam::seulex::seulex(const ODESystem& ode, const dictionary& dict) yTemp_(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; 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; } - //NOTE: the first element of relTol_ and absTol_ are used here. - scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5; - kTarg_ = min(1, min(kMaxx_ - 1, int(logTol))); - + // Set the extrapolation coefficients array for (int k=0; k 0 ? true : false; - + scalar dx = step.dxTry; 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; - lastStep = false; - theta_=2.0*jacRedo_; + // NOTE: the first element of relTol_ and absTol_ are used here. + scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5; + kTarg_ = max(1, min(kMaxx_ - 1, int(logTol))); } forAll(scale_, i) @@ -241,39 +232,39 @@ void Foam::seulex::solve scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]); } - reject = false; - bool firstk = true; - scalar dxNew = mag(dx); + bool jacUpdated = false; - if (theta_ > jacRedo_ && !calcJac_) + if (theta_ > jacRedo_) { odes_.jacobian(x, y, dfdx_, dfdy_); - calcJac_ = true; + jacUpdated = true; } - 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; - reject = false; + step.reject = false; if (mag(dx) <= mag(x)*SMALL) { - WarningIn("seulex::step(const scalar") - << "step size underflow in step :" << dx << endl; + WarningIn("seulex::solve(scalar& x, scalarField& y, stepState&") + << "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_); if (!success) { - reject = true; + step.reject = true; dxNew = mag(dx)*stepFactor5_; break; } @@ -289,6 +280,7 @@ void Foam::seulex::solve table_[k-1][i] = ySequence_[i]; } } + if (k != 0) { extrapolate(k, table_, y); @@ -301,7 +293,7 @@ void Foam::seulex::solve err = sqrt(err/n_); if (err > 1.0/SMALL || (k > 1 && err >= errOld)) { - reject = true; + step.reject = true; dxNew = mag(dx)*stepFactor5_; break; } @@ -321,11 +313,17 @@ void Foam::seulex::solve dxOpt_[k] = mag(dx*fac); temp_[k] = cpu_[k]/dxOpt_[k]; - if ((firstStep || lastStep) && err <= 1.0) + if ((step.first || step.last) && err <= 1.0) { break; } - if (k == kTarg_-1 && !prevReject && !firstStep && !lastStep) + + if + ( + k == kTarg_ - 1 + && !step.prevReject + && !step.first && !step.last + ) { if (err <= 1.0) { @@ -333,7 +331,7 @@ void Foam::seulex::solve } else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4.0) { - reject = true; + step.reject = true; kTarg_ = k; if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k]) { @@ -343,6 +341,7 @@ void Foam::seulex::solve break; } } + if (k == kTarg_) { if (err <= 1.0) @@ -351,7 +350,7 @@ void Foam::seulex::solve } else if (err > nSeq_[k + 1]*2.0) { - reject = true; + step.reject = true; if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k]) { kTarg_--; @@ -360,11 +359,12 @@ void Foam::seulex::solve break; } } + if (k == kTarg_+1) { if (err > 1.0) { - reject = true; + step.reject = true; if ( kTarg_ > 1 @@ -379,26 +379,26 @@ void Foam::seulex::solve } } } - if (reject) + if (step.reject) { - prevReject = true; - if (!calcJac_) + step.prevReject = true; + if (!jacUpdated) { theta_ = 2.0*jacRedo_; - if (theta_ > jacRedo_ && !calcJac_) + if (theta_ > jacRedo_ && !jacUpdated) { odes_.jacobian(x, y, dfdx_, dfdy_); - calcJac_ = true; + jacUpdated = true; } } } } - calcJac_ = false; + jacUpdated = false; + step.dxDid = dx; x += dx; - firstStep = false; label kopt; if (k == 1) @@ -430,11 +430,11 @@ void Foam::seulex::solve } } - if (prevReject) + if (step.prevReject) { kTarg_ = min(kopt, k); dxNew = min(mag(dx), dxOpt_[kTarg_]); - prevReject = false; + step.prevReject = false; } else { @@ -456,14 +456,7 @@ void Foam::seulex::solve kTarg_ = kopt; } - if (forward) - { - dxTry = dxNew; - } - else - { - dxTry = -dxNew; - } + step.dxTry = step.forward ? dxNew : -dxNew; } diff --git a/src/ODE/ODESolvers/seulex/seulex.H b/src/ODE/ODESolvers/seulex/seulex.H index 8c01052331..c3cd540198 100644 --- a/src/ODE/ODESolvers/seulex/seulex.H +++ b/src/ODE/ODESolvers/seulex/seulex.H @@ -28,7 +28,7 @@ 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 + The implementation is based on the SEULEX algorithm in \verbatim "Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, second edition", @@ -65,35 +65,42 @@ class seulex { // Private data - // SEULEX constants - static const label kMaxx_ = 12; - static const label iMaxx_ = kMaxx_ + 1; + // Static constants - static const scalar - stepFactor1_, stepFactor2_, stepFactor3_, - stepFactor4_, stepFactor5_, - kFactor1_, kFactor2_; + static const label kMaxx_ = 12; + static const label iMaxx_ = kMaxx_ + 1; - // Temporary fields - mutable label kTarg_; - mutable labelField nSeq_; - mutable scalarField cpu_; - mutable scalarRectangularMatrix table_; + static const scalar + stepFactor1_, stepFactor2_, stepFactor3_, + stepFactor4_, stepFactor5_, + kFactor1_, kFactor2_; - mutable scalar jacRedo_, theta_; - mutable bool calcJac_; + // Evaluated constants - mutable scalarField dfdx_; - mutable scalarSquareMatrix dfdy_; - mutable scalarSquareMatrix a_, coeff_; - mutable labelList pivotIndices_; + scalar jacRedo_; + labelField nSeq_; + scalarField cpu_; + scalarSquareMatrix coeff_; - // Fields space for "solve" function - mutable scalarField dxOpt_, temp_; - mutable scalarField y0_, ySequence_, scale_; + // Temporary storage + // held to avoid dynamic memory allocation between calls + // and to transfer internal values between functions - // Fields used in "dy" function - mutable scalarField dy_, yTemp_, dydx_; + mutable scalar theta_; + 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 @@ -137,7 +144,7 @@ public: ( scalar& x, scalarField& y, - scalar& dxTry + stepState& step ) const; };