36 const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
41 SIBS::redMax = 1.0e-5,
54 d_p_(n_, kMaxX_, 0.0),
102 bool exitflag =
false;
106 dxTry = xNew_ = -great;
107 scalar eps1 = safe1*
relTol_[0];
108 a_[0] = nSeq_[0] + 1;
112 a_[
k + 1] = a_[
k] + nSeq_[
k + 1];
115 for (
label iq = 1; iq<kMaxX_; iq++)
120 pow(eps1, (a_[
k + 1] - a_[iq + 1])
121 /((a_[iq + 1] - a_[0] + 1.0)*(2*
k + 3)));
125 epsOld_ = relTol_[0];
130 a_[
k + 1] = a_[
k] + nSeq_[
k + 1];
133 for (kOpt_ = 1; kOpt_<kMaxX_ - 1; kOpt_++)
135 if (a_[kOpt_ + 1] > a_[kOpt_]*alpha_[kOpt_ - 1][kOpt_])
149 if (x != xNew_ || h != dxTry)
157 scalar maxErr = small;
163 for (k=0; k <= kMax_; k++)
170 <<
"step size underflow" 174 SIMPR(x, yTemp_, li, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
175 scalar xest =
sqr(h/nSeq_[k]);
177 polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
191 err_[km] =
pow(maxErr/safe1, 1.0/(2*km + 3));
194 if (k != 0 && (k >= kOpt_ - 1 || first_))
202 if (k == kMax_ || k == kOpt_ + 1)
204 red = safe2/err_[km];
207 else if (k == kOpt_ && alpha_[kOpt_ - 1][kOpt_] < err_[km])
212 else if (kOpt_ == kMax_ && alpha_[km][kMax_ - 1] < err_[km])
214 red = alpha_[km][kMax_ - 1]*safe2/err_[km];
217 else if (alpha_[km][kOpt_] < err_[km])
219 red = alpha_[km][kOpt_ - 1]/err_[km];
230 red =
min(red, redMin);
231 red =
max(red, redMax);
238 scalar wrkmin = great;
241 for (
label kk=0; kk<=km; kk++)
243 scalar fact =
max(err_[kk], scaleMX);
244 scalar work = fact*a_[kk + 1];
255 if (kOpt_ >= k && kOpt_ != kMax_ && !reduct)
257 scalar fact =
max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX);
258 if (a_[kOpt_ + 1]*fact <= wrkmin)
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Abstract base class for the systems of ordinary differential equations.
virtual bool resize()=0
Resize the ODE solver.
errorManipArg< error, int > exit(error &err, const int errNo=1)
A list of keyword definitions, which are a keyword followed by any number of values (e...
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
virtual void derivatives(const scalar x, const scalarField &y, const label li, scalarField &dydx) const =0
Calculate the derivatives in dydx.
label k
Boltzmann constant.
An ODE solver for chemistry.
Macros for easy insertion into run-time selection tables.
void resizeMatrix(scalarSquareMatrix &m) const
const ODESystem & odes_
Reference to ODESystem.
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual void solve(scalar &x, scalarField &y, const label li, scalar &dxTry) const
Solve the ODE system from the current state xStart, y.
SIBS(const ODESystem &ode, const dictionary &dict)
Construct from ODE system.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
scalarField relTol_
Relative convergence tolerance per step.
Abstract base-class for ODE system solvers.
virtual void jacobian(const scalar x, const scalarField &y, const label li, scalarField &dfdx, scalarSquareMatrix &dfdy) const =0
Calculate the Jacobian of the system.
scalarField absTol_
Absolute convergence tolerance per step.
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionedScalar & h
Planck constant.
static void resizeField(UList< Type > &f, const label n)
virtual bool resize()
Resize the ODE solver.
label n_
Size of the ODESystem (adjustable)