39 const label SIBS::nSeq_[iMaxX_] = {2, 6, 10, 14, 22, 34, 50, 70};
44 SIBS::redMax = 1.0e-5,
57 d_p_(n_, kMaxX_,
Zero),
99 odes_.derivatives(
x,
y, dydx0_);
102 bool exitflag =
false;
104 if (relTol_[0] != epsOld_)
106 dxTry = xNew_ = -GREAT;
107 scalar eps1 = safe1*relTol_[0];
108 a_[0] = nSeq_[0] + 1;
110 for (label
k=0;
k<kMaxX_;
k++)
112 a_[
k + 1] = a_[
k] + nSeq_[
k + 1];
115 for (label iq = 1; iq<kMaxX_; iq++)
117 for (label
k=0;
k<iq;
k++)
120 pow(eps1, (a_[
k + 1] - a_[iq + 1])
121 /((a_[iq + 1] - a_[0] + 1.0)*(2*
k + 3)));
125 epsOld_ = relTol_[0];
128 for (label
k=0;
k<kMaxX_;
k++)
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_])
147 odes_.jacobian(
x,
y, dfdx_, dfdy_);
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_, 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_);
182 for (label i=0; i<n_; i++)
187 mag(yErr_[i])/(absTol_[i] + relTol_[i]*
mag(yTemp_[i]))
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)
virtual void solve(scalar &x, scalarField &y, scalar &dxTry) const
Solve the ODE system as far as possible up to dxTry.
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)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
label k
Boltzmann constant.
An ODE solver for chemistry.
Macros for easy insertion into run-time selection tables.
void resizeMatrix(scalarSquareMatrix &m) const
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
defineTypeNameAndDebug(combustionModel, 0)
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
const dimensionedScalar h
Planck constant.
SIBS(const ODESystem &ode, const dictionary &dict)
Construct from ODE system.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Abstract base-class for ODE system solvers.
static void resizeField(UList< Type > &f, const label n)
virtual bool resize()
Resize the ODE solver.
static constexpr const zero Zero
Global zero (0)