seulex.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2013-2016 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "seulex.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineTypeNameAndDebug(seulex, 0);
34 
35  addToRunTimeSelectionTable(ODESolver, seulex, dictionary);
36 
37  const scalar
38  seulex::stepFactor1_ = 0.6,
39  seulex::stepFactor2_ = 0.93,
40  seulex::stepFactor3_ = 0.1,
41  seulex::stepFactor4_ = 4,
42  seulex::stepFactor5_ = 0.5,
43  seulex::kFactor1_ = 0.7,
44  seulex::kFactor2_ = 0.9;
45 }
46 
47 
48 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
49 
51 :
52  ODESolver(ode, dict),
53  jacRedo_(min(1e-4, min(relTol_))),
54  nSeq_(iMaxx_),
55  cpu_(iMaxx_),
56  coeff_(iMaxx_, iMaxx_),
57  theta_(2*jacRedo_),
58  table_(kMaxx_, n_),
59  dfdx_(n_),
60  dfdy_(n_),
61  a_(n_),
62  pivotIndices_(n_),
63  dxOpt_(iMaxx_),
64  temp_(iMaxx_),
65  y0_(n_),
66  ySequence_(n_),
67  scale_(n_),
68  dy_(n_),
69  yTemp_(n_),
70  dydx_(n_)
71 {
72  // The CPU time factors for the major parts of the algorithm
73  const scalar cpuFunc = 1, cpuJac = 5, cpuLU = 1, cpuSolve = 1;
74 
75  nSeq_[0] = 2;
76  nSeq_[1] = 3;
77 
78  for (int i=2; i<iMaxx_; i++)
79  {
80  nSeq_[i] = 2*nSeq_[i-2];
81  }
82  cpu_[0] = cpuJac + cpuLU + nSeq_[0]*(cpuFunc + cpuSolve);
83 
84  for (int k=0; k<kMaxx_; k++)
85  {
86  cpu_[k+1] = cpu_[k] + (nSeq_[k+1]-1)*(cpuFunc + cpuSolve) + cpuLU;
87  }
88 
89  // Set the extrapolation coefficients array
90  for (int k=0; k<iMaxx_; k++)
91  {
92  for (int l=0; l<k; l++)
93  {
94  scalar ratio = scalar(nSeq_[k])/nSeq_[l];
95  coeff_(k, l) = 1/(ratio - 1);
96  }
97  }
98 }
99 
100 
101 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
102 
103 bool Foam::seulex::seul
104 (
105  const scalar x0,
106  const scalarField& y0,
107  const scalar dxTot,
108  const label k,
109  scalarField& y,
110  const scalarField& scale
111 ) const
112 {
113  label nSteps = nSeq_[k];
114  scalar dx = dxTot/nSteps;
115 
116  for (label i=0; i<n_; i++)
117  {
118  for (label j=0; j<n_; j++)
119  {
120  a_(i, j) = -dfdy_(i, j);
121  }
122 
123  a_(i, i) += 1/dx;
124  }
125 
126  LUDecompose(a_, pivotIndices_);
127 
128  scalar xnew = x0 + dx;
129  odes_.derivatives(xnew, y0, dy_);
130  LUBacksubstitute(a_, pivotIndices_, dy_);
131 
132  yTemp_ = y0;
133 
134  for (label nn=1; nn<nSteps; nn++)
135  {
136  yTemp_ += dy_;
137  xnew += dx;
138 
139  if (nn == 1 && k<=1)
140  {
141  scalar dy1 = 0;
142  for (label i=0; i<n_; i++)
143  {
144  dy1 += sqr(dy_[i]/scale[i]);
145  }
146  dy1 = sqrt(dy1);
147 
148  odes_.derivatives(x0 + dx, yTemp_, dydx_);
149  for (label i=0; i<n_; i++)
150  {
151  dy_[i] = dydx_[i] - dy_[i]/dx;
152  }
153 
154  LUBacksubstitute(a_, pivotIndices_, dy_);
155 
156  const scalar denom = min(1, dy1 + SMALL);
157  scalar dy2 = 0;
158  for (label i=0; i<n_; i++)
159  {
160  // Test of dy_[i] to avoid overflow
161  if (mag(dy_[i]) > scale[i]*denom)
162  {
163  theta_ = 1;
164  return false;
165  }
166 
167  dy2 += sqr(dy_[i]/scale[i]);
168  }
169  dy2 = sqrt(dy2);
170  theta_ = dy2/denom;
171 
172  if (theta_ > 1)
173  {
174  return false;
175  }
176  }
177 
178  odes_.derivatives(xnew, yTemp_, dy_);
179  LUBacksubstitute(a_, pivotIndices_, dy_);
180  }
181 
182  for (label i=0; i<n_; i++)
183  {
184  y[i] = yTemp_[i] + dy_[i];
185  }
186 
187  return true;
188 }
189 
190 
191 void Foam::seulex::extrapolate
192 (
193  const label k,
195  scalarField& y
196 ) const
197 {
198  for (int j=k-1; j>0; j--)
199  {
200  for (label i=0; i<n_; i++)
201  {
202  table[j-1][i] =
203  table(j, i) + coeff_(k, j)*(table(j, i) - table[j-1][i]);
204  }
205  }
206 
207  for (int i=0; i<n_; i++)
208  {
209  y[i] = table(0, i) + coeff_(k, 0)*(table(0, i) - y[i]);
210  }
211 }
212 
213 
215 {
216  if (ODESolver::resize())
217  {
218  table_.shallowResize(kMaxx_, n_);
219  resizeField(dfdx_);
220  resizeMatrix(dfdy_);
221  resizeMatrix(a_);
222  resizeField(pivotIndices_);
223  resizeField(y0_);
224  resizeField(ySequence_);
225  resizeField(scale_);
226  resizeField(dy_);
227  resizeField(yTemp_);
228  resizeField(dydx_);
229 
230  return true;
231  }
232  else
233  {
234  return false;
235  }
236 }
237 
238 
240 (
241  scalar& x,
242  scalarField& y,
243  stepState& step
244 ) const
245 {
246  temp_[0] = GREAT;
247  scalar dx = step.dxTry;
248  y0_ = y;
249  dxOpt_[0] = mag(0.1*dx);
250 
251  if (step.first || step.prevReject)
252  {
253  theta_ = 2*jacRedo_;
254  }
255 
256  if (step.first)
257  {
258  // NOTE: the first element of relTol_ and absTol_ are used here.
259  scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
260  kTarg_ = max(1, min(kMaxx_ - 1, int(logTol)));
261  }
262 
263  forAll(scale_, i)
264  {
265  scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
266  }
267 
268  bool jacUpdated = false;
269 
270  if (theta_ > jacRedo_)
271  {
272  odes_.jacobian(x, y, dfdx_, dfdy_);
273  jacUpdated = true;
274  }
275 
276  int k;
277  scalar dxNew = mag(dx);
278  bool firstk = true;
279 
280  while (firstk || step.reject)
281  {
282  dx = step.forward ? dxNew : -dxNew;
283  firstk = false;
284  step.reject = false;
285 
286  if (mag(dx) <= mag(x)*sqr(SMALL))
287  {
289  << "step size underflow :" << dx << endl;
290  }
291 
292  scalar errOld = 0;
293 
294  for (k=0; k<=kTarg_+1; k++)
295  {
296  bool success = seul(x, y0_, dx, k, ySequence_, scale_);
297 
298  if (!success)
299  {
300  step.reject = true;
301  dxNew = mag(dx)*stepFactor5_;
302  break;
303  }
304 
305  if (k == 0)
306  {
307  y = ySequence_;
308  }
309  else
310  {
311  forAll(ySequence_, i)
312  {
313  table_[k-1][i] = ySequence_[i];
314  }
315  }
316 
317  if (k != 0)
318  {
319  extrapolate(k, table_, y);
320  scalar err = 0;
321  forAll(scale_, i)
322  {
323  scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
324  err += sqr((y[i] - table_(0, i))/scale_[i]);
325  }
326  err = sqrt(err/n_);
327  if (err > 1/SMALL || (k > 1 && err >= errOld))
328  {
329  step.reject = true;
330  dxNew = mag(dx)*stepFactor5_;
331  break;
332  }
333  errOld = min(4*err, 1);
334  scalar expo = 1.0/(k + 1);
335  scalar facmin = pow(stepFactor3_, expo);
336  scalar fac;
337  if (err == 0)
338  {
339  fac = 1/facmin;
340  }
341  else
342  {
343  fac = stepFactor2_/pow(err/stepFactor1_, expo);
344  fac = max(facmin/stepFactor4_, min(1/facmin, fac));
345  }
346  dxOpt_[k] = mag(dx*fac);
347  temp_[k] = cpu_[k]/dxOpt_[k];
348 
349  if ((step.first || step.last) && err <= 1)
350  {
351  break;
352  }
353 
354  if
355  (
356  k == kTarg_ - 1
357  && !step.prevReject
358  && !step.first && !step.last
359  )
360  {
361  if (err <= 1)
362  {
363  break;
364  }
365  else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4)
366  {
367  step.reject = true;
368  kTarg_ = k;
369  if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
370  {
371  kTarg_--;
372  }
373  dxNew = dxOpt_[kTarg_];
374  break;
375  }
376  }
377 
378  if (k == kTarg_)
379  {
380  if (err <= 1)
381  {
382  break;
383  }
384  else if (err > nSeq_[k + 1]*2)
385  {
386  step.reject = true;
387  if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
388  {
389  kTarg_--;
390  }
391  dxNew = dxOpt_[kTarg_];
392  break;
393  }
394  }
395 
396  if (k == kTarg_+1)
397  {
398  if (err > 1)
399  {
400  step.reject = true;
401  if
402  (
403  kTarg_ > 1
404  && temp_[kTarg_-1] < kFactor1_*temp_[kTarg_]
405  )
406  {
407  kTarg_--;
408  }
409  dxNew = dxOpt_[kTarg_];
410  }
411  break;
412  }
413  }
414  }
415  if (step.reject)
416  {
417  step.prevReject = true;
418  if (!jacUpdated)
419  {
420  theta_ = 2*jacRedo_;
421 
422  if (theta_ > jacRedo_ && !jacUpdated)
423  {
424  odes_.jacobian(x, y, dfdx_, dfdy_);
425  jacUpdated = true;
426  }
427  }
428  }
429  }
430 
431  jacUpdated = false;
432 
433  step.dxDid = dx;
434  x += dx;
435 
436  label kopt;
437  if (k == 1)
438  {
439  kopt = 2;
440  }
441  else if (k <= kTarg_)
442  {
443  kopt=k;
444  if (temp_[k-1] < kFactor1_*temp_[k])
445  {
446  kopt = k - 1;
447  }
448  else if (temp_[k] < kFactor2_*temp_[k - 1])
449  {
450  kopt = min(k + 1, kMaxx_ - 1);
451  }
452  }
453  else
454  {
455  kopt = k - 1;
456  if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1])
457  {
458  kopt = k - 2;
459  }
460  if (temp_[k] < kFactor2_*temp_[kopt])
461  {
462  kopt = min(k, kMaxx_ - 1);
463  }
464  }
465 
466  if (step.prevReject)
467  {
468  kTarg_ = min(kopt, k);
469  dxNew = min(mag(dx), dxOpt_[kTarg_]);
470  step.prevReject = false;
471  }
472  else
473  {
474  if (kopt <= k)
475  {
476  dxNew = dxOpt_[kopt];
477  }
478  else
479  {
480  if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1])
481  {
482  dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k];
483  }
484  else
485  {
486  dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k];
487  }
488  }
489  kTarg_ = kopt;
490  }
491 
492  step.dxTry = step.forward ? dxNew : -dxNew;
493 }
494 
495 
496 // ************************************************************************* //
void shallowResize(const label m, const label n)
Resize the matrix without reallocating storage (unsafe)
Definition: MatrixI.H:284
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:46
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:89
const double e
Elementary charge.
Definition: doubleFloat.H:78
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
dimensionedScalar y0(const dimensionedScalar &ds)
label k
Boltzmann constant.
virtual void solve(scalar &x, scalarField &y, stepState &step) const
Solve the ODE system and the update the state.
Definition: seulex.C:240
An ODE solver for chemistry.
Definition: ode.H:50
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.
scalar y
virtual bool resize()
Resize the ODE solver.
Definition: seulex.C:214
void resizeMatrix(scalarSquareMatrix &m) const
Definition: ODESolverI.H:61
const ODESystem & odes_
Reference to ODESystem.
Definition: ODESolver.H:58
addToRunTimeSelectionTable(ensightPart, ensightPartCells, istream)
bool success
defineTypeNameAndDebug(combustionModel, 0)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
#define WarningInFunction
Report a warning using Foam::Warning.
scalarField relTol_
Relative convergence tolerance per step.
Definition: ODESolver.H:70
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:50
scalarField absTol_
Absolute convergence tolerance per step.
Definition: ODESolver.H:67
seulex(const ODESystem &ode, const dictionary &dict)
Construct from ODESystem.
Definition: seulex.C:50
dimensioned< scalar > mag(const dimensioned< Type > &)
static void resizeField(UList< Type > &f, const label n)
Definition: ODESolverI.H:48
label n_
Size of the ODESystem (adjustable)
Definition: ODESolver.H:64
dimensionedScalar log10(const dimensionedScalar &ds)
virtual void derivatives(const scalar x, const scalarField &y, scalarField &dydx) const =0
Calculate the derivatives in dydx.
Namespace for OpenFOAM.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution.