seulex.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2013-2019 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 {
34 
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 label li,
108  const scalar dxTot,
109  const label k,
110  scalarField& y,
111  const scalarField& scale
112 ) const
113 {
114  label nSteps = nSeq_[k];
115  scalar dx = dxTot/nSteps;
116 
117  for (label i=0; i<n_; i++)
118  {
119  for (label j=0; j<n_; j++)
120  {
121  a_(i, j) = -dfdy_(i, j);
122  }
123 
124  a_(i, i) += 1/dx;
125  }
126 
127  LUDecompose(a_, pivotIndices_);
128 
129  scalar xnew = x0 + dx;
130  odes_.derivatives(xnew, y0, li, dy_);
131  LUBacksubstitute(a_, pivotIndices_, dy_);
132 
133  yTemp_ = y0;
134 
135  for (label nn=1; nn<nSteps; nn++)
136  {
137  yTemp_ += dy_;
138  xnew += dx;
139 
140  if (nn == 1 && k<=1)
141  {
142  scalar dy1 = 0;
143  for (label i=0; i<n_; i++)
144  {
145  dy1 += sqr(dy_[i]/scale[i]);
146  }
147  dy1 = sqrt(dy1);
148 
149  odes_.derivatives(x0 + dx, yTemp_, li, dydx_);
150  for (label i=0; i<n_; i++)
151  {
152  dy_[i] = dydx_[i] - dy_[i]/dx;
153  }
154 
155  LUBacksubstitute(a_, pivotIndices_, dy_);
156 
157  // This form from the original paper is unreliable
158  // step size underflow for some cases
159  // const scalar denom = max(1, dy1);
160 
161  // This form is reliable but limits how large the step size can be
162  const scalar denom = min(1, dy1 + small);
163 
164  scalar dy2 = 0;
165  for (label i=0; i<n_; i++)
166  {
167  // Test of dy_[i] to avoid overflow
168  if (mag(dy_[i]) > scale[i]*denom)
169  {
170  theta_ = 1;
171  return false;
172  }
173 
174  dy2 += sqr(dy_[i]/scale[i]);
175  }
176  dy2 = sqrt(dy2);
177  theta_ = dy2/denom;
178 
179  if (theta_ > 1)
180  {
181  return false;
182  }
183  }
184 
185  odes_.derivatives(xnew, yTemp_, li, dy_);
186  LUBacksubstitute(a_, pivotIndices_, dy_);
187  }
188 
189  for (label i=0; i<n_; i++)
190  {
191  y[i] = yTemp_[i] + dy_[i];
192  }
193 
194  return true;
195 }
196 
197 
198 void Foam::seulex::extrapolate
199 (
200  const label k,
202  scalarField& y
203 ) const
204 {
205  for (int j=k-1; j>0; j--)
206  {
207  for (label i=0; i<n_; i++)
208  {
209  table[j-1][i] =
210  table(j, i) + coeff_(k, j)*(table(j, i) - table[j-1][i]);
211  }
212  }
213 
214  for (int i=0; i<n_; i++)
215  {
216  y[i] = table(0, i) + coeff_(k, 0)*(table(0, i) - y[i]);
217  }
218 }
219 
220 
222 {
223  if (ODESolver::resize())
224  {
225  table_.shallowResize(kMaxx_, n_);
226  resizeField(dfdx_);
227  resizeMatrix(dfdy_);
228  resizeMatrix(a_);
229  resizeField(pivotIndices_);
230  resizeField(y0_);
231  resizeField(ySequence_);
232  resizeField(scale_);
233  resizeField(dy_);
234  resizeField(yTemp_);
235  resizeField(dydx_);
236 
237  return true;
238  }
239  else
240  {
241  return false;
242  }
243 }
244 
245 
247 (
248  scalar& x,
249  scalarField& y,
250  const label li,
251  stepState& step
252 ) const
253 {
254  temp_[0] = great;
255  scalar dx = step.dxTry;
256  y0_ = y;
257  dxOpt_[0] = mag(0.1*dx);
258 
259  if (step.first || step.prevReject)
260  {
261  theta_ = 2*jacRedo_;
262  }
263 
264  if (step.first)
265  {
266  // NOTE: the first element of relTol_ and absTol_ are used here.
267  scalar logTol = -log10(relTol_[0] + absTol_[0])*0.6 + 0.5;
268  kTarg_ = max(1, min(kMaxx_ - 1, int(logTol)));
269  }
270 
271  forAll(scale_, i)
272  {
273  scale_[i] = absTol_[i] + relTol_[i]*mag(y[i]);
274  }
275 
276  bool jacUpdated = false;
277 
278  if (theta_ > jacRedo_)
279  {
280  odes_.jacobian(x, y, li, dfdx_, dfdy_);
281  jacUpdated = true;
282  }
283 
284  int k;
285  scalar dxNew = mag(dx);
286  bool firstk = true;
287 
288  while (firstk || step.reject)
289  {
290  dx = step.forward ? dxNew : -dxNew;
291  firstk = false;
292  step.reject = false;
293 
294  if (mag(dx) <= mag(x)*sqr(small))
295  {
297  << "step size underflow :" << dx << endl;
298  }
299 
300  scalar errOld = 0;
301 
302  for (k=0; k<=kTarg_+1; k++)
303  {
304  bool success = seul(x, y0_, li, dx, k, ySequence_, scale_);
305 
306  if (!success)
307  {
308  step.reject = true;
309  dxNew = mag(dx)*stepFactor5_;
310  break;
311  }
312 
313  if (k == 0)
314  {
315  y = ySequence_;
316  }
317  else
318  {
319  forAll(ySequence_, i)
320  {
321  table_[k-1][i] = ySequence_[i];
322  }
323  }
324 
325  if (k != 0)
326  {
327  extrapolate(k, table_, y);
328  scalar err = 0;
329  forAll(scale_, i)
330  {
331  scale_[i] = absTol_[i] + relTol_[i]*mag(y0_[i]);
332  err += sqr((y[i] - table_(0, i))/scale_[i]);
333  }
334  err = sqrt(err/n_);
335  if (err > 1/small || (k > 1 && err >= errOld))
336  {
337  step.reject = true;
338  dxNew = mag(dx)*stepFactor5_;
339  break;
340  }
341  errOld = min(4*err, 1);
342  scalar expo = 1.0/(k + 1);
343  scalar facmin = pow(stepFactor3_, expo);
344  scalar fac;
345  if (err == 0)
346  {
347  fac = 1/facmin;
348  }
349  else
350  {
351  fac = stepFactor2_/pow(err/stepFactor1_, expo);
352  fac = max(facmin/stepFactor4_, min(1/facmin, fac));
353  }
354  dxOpt_[k] = mag(dx*fac);
355  temp_[k] = cpu_[k]/dxOpt_[k];
356 
357  if ((step.first || step.last) && err <= 1)
358  {
359  break;
360  }
361 
362  if
363  (
364  k == kTarg_ - 1
365  && !step.prevReject
366  && !step.first && !step.last
367  )
368  {
369  if (err <= 1)
370  {
371  break;
372  }
373  else if (err > nSeq_[kTarg_]*nSeq_[kTarg_ + 1]*4)
374  {
375  step.reject = true;
376  kTarg_ = k;
377  if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
378  {
379  kTarg_--;
380  }
381  dxNew = dxOpt_[kTarg_];
382  break;
383  }
384  }
385 
386  if (k == kTarg_)
387  {
388  if (err <= 1)
389  {
390  break;
391  }
392  else if (err > nSeq_[k + 1]*2)
393  {
394  step.reject = true;
395  if (kTarg_>1 && temp_[k-1] < kFactor1_*temp_[k])
396  {
397  kTarg_--;
398  }
399  dxNew = dxOpt_[kTarg_];
400  break;
401  }
402  }
403 
404  if (k == kTarg_+1)
405  {
406  if (err > 1)
407  {
408  step.reject = true;
409  if
410  (
411  kTarg_ > 1
412  && temp_[kTarg_-1] < kFactor1_*temp_[kTarg_]
413  )
414  {
415  kTarg_--;
416  }
417  dxNew = dxOpt_[kTarg_];
418  }
419  break;
420  }
421  }
422  }
423  if (step.reject)
424  {
425  step.prevReject = true;
426  if (!jacUpdated)
427  {
428  theta_ = 2*jacRedo_;
429 
430  if (theta_ > jacRedo_ && !jacUpdated)
431  {
432  odes_.jacobian(x, y, li, dfdx_, dfdy_);
433  jacUpdated = true;
434  }
435  }
436  }
437  }
438 
439  jacUpdated = false;
440 
441  step.dxDid = dx;
442  x += dx;
443 
444  label kopt;
445  if (k == 1)
446  {
447  kopt = 2;
448  }
449  else if (k <= kTarg_)
450  {
451  kopt=k;
452  if (temp_[k-1] < kFactor1_*temp_[k])
453  {
454  kopt = k - 1;
455  }
456  else if (temp_[k] < kFactor2_*temp_[k - 1])
457  {
458  kopt = min(k + 1, kMaxx_ - 1);
459  }
460  }
461  else
462  {
463  kopt = k - 1;
464  if (k > 2 && temp_[k-2] < kFactor1_*temp_[k - 1])
465  {
466  kopt = k - 2;
467  }
468  if (temp_[k] < kFactor2_*temp_[kopt])
469  {
470  kopt = min(k, kMaxx_ - 1);
471  }
472  }
473 
474  if (step.prevReject)
475  {
476  kTarg_ = min(kopt, k);
477  dxNew = min(mag(dx), dxOpt_[kTarg_]);
478  step.prevReject = false;
479  }
480  else
481  {
482  if (kopt <= k)
483  {
484  dxNew = dxOpt_[kopt];
485  }
486  else
487  {
488  if (k < kTarg_ && temp_[k] < kFactor2_*temp_[k - 1])
489  {
490  dxNew = dxOpt_[k]*cpu_[kopt + 1]/cpu_[k];
491  }
492  else
493  {
494  dxNew = dxOpt_[k]*cpu_[kopt]/cpu_[k];
495  }
496  }
497  kTarg_ = kopt;
498  }
499 
500  step.dxTry = step.forward ? dxNew : -dxNew;
501 }
502 
503 
504 // ************************************************************************* //
bool success
scalar y
label k
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
Macros for easy insertion into run-time selection tables.
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:51
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:89
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:47
A list of keyword definitions, which are a keyword followed by any number of values (e....
Definition: dictionary.H:162
An ODE solver for chemistry.
Definition: ode.H:53
An extrapolation-algorithm, based on the linearly implicit Euler method with step size control and or...
Definition: seulex.H:62
seulex(const ODESystem &ode, const dictionary &dict)
Construct from ODESystem.
Definition: seulex.C:50
virtual bool resize()
Resize the ODE solver.
Definition: seulex.C:221
virtual void solve(scalar &x, scalarField &y, const label li, stepState &step) const
Solve the ODE system and the update the state.
Definition: seulex.C:247
#define WarningInFunction
Report a warning using Foam::Warning.
const dimensionedScalar e
Elementary charge.
Namespace for OpenFOAM.
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
addToRunTimeSelectionTable(polyPatch, mergedCyclicPolyPatch, word)
RectangularMatrix< scalar > scalarRectangularMatrix
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:257
dimensionedScalar y0(const dimensionedScalar &ds)
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution.
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
dimensionedScalar log10(const dimensionedScalar &ds)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< scalar > mag(const dimensioned< Type > &)
defineTypeNameAndDebug(combustionModel, 0)
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
dictionary dict