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  scalar& x,
217  scalarField& y,
218  stepState& step
219 ) const
220 {
221  temp_[0] = GREAT;
222  scalar dx = step.dxTry;
223  y0_ = y;
224  dxOpt_[0] = mag(0.1*dx);
225 
226  if (step.first || step.prevReject)
227  {
228  theta_ = 2*jacRedo_;
229  }
230 
231  if (step.first)
232  {
233  // NOTE: the first element of relTol_ and absTol_ are used here.
234  scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
235  kTarg_ = max(1, min(kMaxx_ - 1, int(logTol)));
236  }
237 
238  forAll(scale_, i)
239  {
240  scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
241  }
242 
243  bool jacUpdated = false;
244 
245  if (theta_ > jacRedo_)
246  {
247  odes_.jacobian(x, y, dfdx_, dfdy_);
248  jacUpdated = true;
249  }
250 
251  int k;
252  scalar dxNew = mag(dx);
253  bool firstk = true;
254 
255  while (firstk || step.reject)
256  {
257  dx = step.forward ? dxNew : -dxNew;
258  firstk = false;
259  step.reject = false;
260 
261  if (mag(dx) <= mag(x)*sqr(SMALL))
262  {
264  << "step size underflow :" << dx << endl;
265  }
266 
267  scalar errOld = 0;
268 
269  for (k=0; k<=kTarg_+1; k++)
270  {
271  bool success = seul(x, y0_, dx, k, ySequence_, scale_);
272 
273  if (!success)
274  {
275  step.reject = true;
276  dxNew = mag(dx)*stepFactor5_;
277  break;
278  }
279 
280  if (k == 0)
281  {
282  y = ySequence_;
283  }
284  else
285  {
286  forAll(ySequence_, i)
287  {
288  table_[k-1][i] = ySequence_[i];
289  }
290  }
291 
292  if (k != 0)
293  {
294  extrapolate(k, table_, y);
295  scalar err = 0;
296  forAll(scale_, i)
297  {
298  scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
299  err += sqr((y[i] - table_(0, i))/scale_[i]);
300  }
301  err = sqrt(err/n_);
302  if (err > 1/SMALL || (k > 1 && err >= errOld))
303  {
304  step.reject = true;
305  dxNew = mag(dx)*stepFactor5_;
306  break;
307  }
308  errOld = min(4*err, 1);
309  scalar expo = 1.0/(k + 1);
310  scalar facmin = pow(stepFactor3_, expo);
311  scalar fac;
312  if (err == 0)
313  {
314  fac = 1/facmin;
315  }
316  else
317  {
318  fac = stepFactor2_/pow(err/stepFactor1_, expo);
319  fac = max(facmin/stepFactor4_, min(1/facmin, fac));
320  }
321  dxOpt_[k] = mag(dx*fac);
322  temp_[k] = cpu_[k]/dxOpt_[k];
323 
324  if ((step.first || step.last) && err <= 1)
325  {
326  break;
327  }
328 
329  if
330  (
331  k == kTarg_ - 1
332  && !step.prevReject
333  && !step.first && !step.last
334  )
335  {
336  if (err <= 1)
337  {
338  break;
339  }
340  else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4)
341  {
342  step.reject = true;
343  kTarg_ = k;
344  if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
345  {
346  kTarg_--;
347  }
348  dxNew = dxOpt_[kTarg_];
349  break;
350  }
351  }
352 
353  if (k == kTarg_)
354  {
355  if (err <= 1)
356  {
357  break;
358  }
359  else if (err > nSeq_[k + 1]*2)
360  {
361  step.reject = true;
362  if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
363  {
364  kTarg_--;
365  }
366  dxNew = dxOpt_[kTarg_];
367  break;
368  }
369  }
370 
371  if (k == kTarg_+1)
372  {
373  if (err > 1)
374  {
375  step.reject = true;
376  if
377  (
378  kTarg_ > 1
379  && temp_[kTarg_-1] < kFactor1_*temp_[kTarg_]
380  )
381  {
382  kTarg_--;
383  }
384  dxNew = dxOpt_[kTarg_];
385  }
386  break;
387  }
388  }
389  }
390  if (step.reject)
391  {
392  step.prevReject = true;
393  if (!jacUpdated)
394  {
395  theta_ = 2*jacRedo_;
396 
397  if (theta_ > jacRedo_ && !jacUpdated)
398  {
399  odes_.jacobian(x, y, dfdx_, dfdy_);
400  jacUpdated = true;
401  }
402  }
403  }
404  }
405 
406  jacUpdated = false;
407 
408  step.dxDid = dx;
409  x += dx;
410 
411  label kopt;
412  if (k == 1)
413  {
414  kopt = 2;
415  }
416  else if (k <= kTarg_)
417  {
418  kopt=k;
419  if (temp_[k-1] < kFactor1_*temp_[k])
420  {
421  kopt = k - 1;
422  }
423  else if (temp_[k] < kFactor2_*temp_[k - 1])
424  {
425  kopt = min(k + 1, kMaxx_ - 1);
426  }
427  }
428  else
429  {
430  kopt = k - 1;
431  if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1])
432  {
433  kopt = k - 2;
434  }
435  if (temp_[k] < kFactor2_*temp_[kopt])
436  {
437  kopt = min(k, kMaxx_ - 1);
438  }
439  }
440 
441  if (step.prevReject)
442  {
443  kTarg_ = min(kopt, k);
444  dxNew = min(mag(dx), dxOpt_[kTarg_]);
445  step.prevReject = false;
446  }
447  else
448  {
449  if (kopt <= k)
450  {
451  dxNew = dxOpt_[kopt];
452  }
453  else
454  {
455  if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1])
456  {
457  dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k];
458  }
459  else
460  {
461  dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k];
462  }
463  }
464  kTarg_ = kopt;
465  }
466 
467  step.dxTry = step.forward ? dxNew : -dxNew;
468 }
469 
470 
471 // ************************************************************************* //
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
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.
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
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 > &)
void solve(scalar &x, scalarField &y, stepState &step) const
Solve the ODE system and the update the state.
Definition: seulex.C:215
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:67
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:50
scalarField absTol_
Absolute convergence tolerance per step.
Definition: ODESolver.H:64
seulex(const ODESystem &ode, const dictionary &dict)
Construct from ODE.
Definition: seulex.C:50
dimensioned< scalar > mag(const dimensioned< Type > &)
label n_
Size of the ODESystem.
Definition: ODESolver.H:61
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.