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),
101 bool exitflag =
false;
105 dxTry = xNew_ = -great;
106 scalar eps1 = safe1*
relTol_[0];
107 a_[0] = nSeq_[0] + 1;
111 a_[
k + 1] = a_[
k] + nSeq_[
k + 1];
114 for (
label iq = 1; iq<kMaxX_; iq++)
119 pow(eps1, (a_[
k + 1] - a_[iq + 1])
120 /((a_[iq + 1] - a_[0] + 1.0)*(2*
k + 3)));
124 epsOld_ = relTol_[0];
129 a_[
k + 1] = a_[
k] + nSeq_[
k + 1];
132 for (kOpt_ = 1; kOpt_<kMaxX_ - 1; kOpt_++)
134 if (a_[kOpt_ + 1] > a_[kOpt_]*alpha_[kOpt_ - 1][kOpt_])
148 if (x != xNew_ || h != dxTry)
156 scalar maxErr = small;
162 for (k=0; k <= kMax_; k++)
169 <<
"step size underflow" 173 SIMPR(x, yTemp_, dydx0_, dfdx_, dfdy_, h, nSeq_[k], ySeq_);
174 scalar xest =
sqr(h/nSeq_[k]);
176 polyExtrapolate(k, xest, ySeq_, y, yErr_, x_p_, d_p_);
190 err_[km] =
pow(maxErr/safe1, 1.0/(2*km + 3));
193 if (k != 0 && (k >= kOpt_ - 1 || first_))
201 if (k == kMax_ || k == kOpt_ + 1)
203 red = safe2/err_[km];
206 else if (k == kOpt_ && alpha_[kOpt_ - 1][kOpt_] < err_[km])
211 else if (kOpt_ == kMax_ && alpha_[km][kMax_ - 1] < err_[km])
213 red = alpha_[km][kMax_ - 1]*safe2/err_[km];
216 else if (alpha_[km][kOpt_] < err_[km])
218 red = alpha_[km][kOpt_ - 1]/err_[km];
229 red =
min(red, redMin);
230 red =
max(red, redMax);
237 scalar wrkmin = great;
240 for (
label kk=0; kk<=km; kk++)
242 scalar fact =
max(err_[kk], scaleMX);
243 scalar work = fact*a_[kk + 1];
254 if (kOpt_ >= k && kOpt_ != kMax_ && !reduct)
256 scalar fact =
max(scale/alpha_[kOpt_ - 1][kOpt_], scaleMX);
257 if (a_[kOpt_ + 1]*fact <= wrkmin)
virtual void solve(scalar &x, scalarField &y, scalar &dxTry) const
Solve the ODE system as far as possible up to dxTry.
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)
label k
Boltzmann constant.
An ODE solver for chemistry.
Macros for easy insertion into run-time selection tables.
virtual void jacobian(const scalar x, const scalarField &y, scalarField &dfdx, scalarSquareMatrix &dfdy) const =0
Calculate the Jacobian of the system.
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 > &)
SIBS(const ODESystem &ode, const dictionary &dict)
Construct from ODE system.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
volScalarField & h
Planck constant.
scalarField relTol_
Relative convergence tolerance per step.
Abstract base-class for ODE system solvers.
scalarField absTol_
Absolute convergence tolerance per step.
dimensioned< scalar > mag(const dimensioned< Type > &)
static void resizeField(UList< Type > &f, const label n)
virtual bool resize()
Resize the ODE solver.
label n_
Size of the ODESystem (adjustable)
virtual void derivatives(const scalar x, const scalarField &y, scalarField &dydx) const =0
Calculate the derivatives in dydx.