Rosenbrock34.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 | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2013-2016 OpenFOAM Foundation
9  Copyright (C) 2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "Rosenbrock34.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36  defineTypeNameAndDebug(Rosenbrock34, 0);
37  addToRunTimeSelectionTable(ODESolver, Rosenbrock34, dictionary);
38 
39 const scalar
40  // Constants by Shampine
41  // More accurate than the L-Stable coefficients for small step-size
42  // but less stable for large step-size
43  Rosenbrock34::a21 = 2,
44  Rosenbrock34::a31 = 48.0/25.0,
45  Rosenbrock34::a32 = 6.0/25.0,
46 
47  Rosenbrock34::c21 = -8,
48  Rosenbrock34::c31 = 372.0/25.0,
49  Rosenbrock34::c32 = 12.0/5.0,
50 
51  Rosenbrock34::c41 = -112.0/125.0,
52  Rosenbrock34::c42 = -54.0/125.0,
53  Rosenbrock34::c43 = -2.0/5.0,
54 
55  Rosenbrock34::b1 = 19.0/9.0,
56  Rosenbrock34::b2 = 1.0/2.0,
57  Rosenbrock34::b3 = 25.0/108.0,
58  Rosenbrock34::b4 = 125.0/108.0,
59 
60  Rosenbrock34::e1 = 34.0/108.0,
61  Rosenbrock34::e2 = 7.0/36.0,
62  Rosenbrock34::e3 = 0,
63  Rosenbrock34::e4 = 125.0/108.0,
64 
65  Rosenbrock34::gamma = 1.0/2.0,
66  Rosenbrock34::c2 = 1,
67  Rosenbrock34::c3 = 3.0/5.0,
68 
69  Rosenbrock34::d1 = 1.0/2.0,
70  Rosenbrock34::d2 = -3.0/2.0,
71  Rosenbrock34::d3 = 605.0/250.0,
72  Rosenbrock34::d4 = 29.0/250.0;
73 
74  /*
75  // L-Stable constants from Hairer et. al.
76  Rosenbrock34::a21 = 2,
77  Rosenbrock34::a31 = 1.867943637803922,
78  Rosenbrock34::a32 = 0.2344449711399156,
79 
80  Rosenbrock34::c21 = -7.137615036412310,
81  Rosenbrock34::c31 = 2.580708087951457,
82  Rosenbrock34::c32 = 0.6515950076447975,
83  Rosenbrock34::c41 = -2.137148994382534,
84  Rosenbrock34::c42 = -0.3214669691237626,
85  Rosenbrock34::c43 = -0.6949742501781779,
86 
87  Rosenbrock34::b1 = 2.255570073418735,
88  Rosenbrock34::b2 = 0.2870493262186792,
89  Rosenbrock34::b3 = 0.435317943184018,
90  Rosenbrock34::b4 = 1.093502252409163,
91 
92  Rosenbrock34::e1 = -0.2815431932141155,
93  Rosenbrock34::e2 = -0.0727619912493892,
94  Rosenbrock34::e3 = -0.1082196201495311,
95  Rosenbrock34::e4 = -1.093502252409163,
96 
97  Rosenbrock34::gamma = 0.57282,
98  Rosenbrock34::c2 = 1.14564,
99  Rosenbrock34::c3 = 0.65521686381559,
100 
101  Rosenbrock34::d1 = 0.57282,
102  Rosenbrock34::d2 = -1.769193891319233,
103  Rosenbrock34::d3 = 0.7592633437920482,
104  Rosenbrock34::d4 = -0.1049021087100450;
105  */
106 }
107 
108 
109 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
110 
112 :
113  ODESolver(ode, dict),
115  k1_(n_),
116  k2_(n_),
117  k3_(n_),
118  k4_(n_),
119  err_(n_),
120  dydx_(n_),
121  dfdx_(n_),
122  dfdy_(n_, n_),
123  a_(n_, n_),
124  pivotIndices_(n_)
125 {}
126 
127 
128 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
129 
131 {
132  if (ODESolver::resize())
133  {
135 
136  resizeField(k1_);
137  resizeField(k2_);
138  resizeField(k3_);
139  resizeField(k4_);
140  resizeField(err_);
141  resizeField(dydx_);
142  resizeField(dfdx_);
143  resizeMatrix(dfdy_);
144  resizeMatrix(a_);
145  resizeField(pivotIndices_);
146 
147  return true;
148  }
149 
150  return false;
151 }
152 
153 
154 Foam::scalar Foam::Rosenbrock34::solve
155 (
156  const scalar x0,
157  const scalarField& y0,
158  const scalarField& dydx0,
159  const scalar dx,
160  scalarField& y
161 ) const
162 {
163  odes_.jacobian(x0, y0, dfdx_, dfdy_);
164 
165  for (label i=0; i<n_; i++)
166  {
167  for (label j=0; j<n_; j++)
168  {
169  a_(i, j) = -dfdy_(i, j);
170  }
171 
172  a_(i, i) += 1.0/(gamma*dx);
173  }
174 
175  LUDecompose(a_, pivotIndices_);
176 
177  // Calculate k1:
178  forAll(k1_, i)
179  {
180  k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
181  }
182 
183  LUBacksubstitute(a_, pivotIndices_, k1_);
184 
185  // Calculate k2:
186  forAll(y, i)
187  {
188  y[i] = y0[i] + a21*k1_[i];
189  }
190 
191  odes_.derivatives(x0 + c2*dx, y, dydx_);
192 
193  forAll(k2_, i)
194  {
195  k2_[i] = dydx_[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
196  }
197 
198  LUBacksubstitute(a_, pivotIndices_, k2_);
199 
200  // Calculate k3:
201  forAll(y, i)
202  {
203  y[i] = y0[i] + a31*k1_[i] + a32*k2_[i];
204  }
205 
206  odes_.derivatives(x0 + c3*dx, y, dydx_);
207 
208  forAll(k3_, i)
209  {
210  k3_[i] = dydx_[i] + dx*d3*dfdx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
211  }
212 
213  LUBacksubstitute(a_, pivotIndices_, k3_);
214 
215  // Calculate k4:
216  forAll(k4_, i)
217  {
218  k4_[i] = dydx_[i] + dx*d4*dfdx_[i]
219  + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
220  }
221 
222  LUBacksubstitute(a_, pivotIndices_, k4_);
223 
224  // Calculate error and update state:
225  forAll(y, i)
226  {
227  y[i] = y0[i] + b1*k1_[i] + b2*k2_[i] + b3*k3_[i] + b4*k4_[i];
228  err_[i] = e1*k1_[i] + e2*k2_[i] + e4*k4_[i];
229  }
230 
231  return normalizeError(y0, y, err_);
232 }
233 
234 
236 (
237  scalar& x,
238  scalarField& y,
239  scalar& dxTry
240 ) const
241 {
242  adaptiveSolver::solve(odes_, x, y, dxTry);
243 }
244 
245 
246 // ************************************************************************* //
Foam::addToRunTimeSelectionTable
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
Foam::ODESolver
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:56
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::adaptiveSolver::solve
virtual scalar solve(const scalar x0, const scalarField &y0, const scalarField &dydx0, const scalar dx, scalarField &y) const =0
Solve a single step dx and return the error.
Foam::Field< scalar >
Foam::Rosenbrock34::resize
virtual bool resize()
Resize the ODE solver.
Definition: Rosenbrock34.C:130
Foam::y0
dimensionedScalar y0(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:281
Foam::Rosenbrock34::Rosenbrock34
Rosenbrock34(const ODESystem &ode, const dictionary &dict)
Construct from ODESystem.
Definition: Rosenbrock34.C:111
Foam::ode
An ODE solver for chemistry.
Definition: ode.H:52
Foam::Rosenbrock34::solve
virtual scalar solve(const scalar x0, const scalarField &y0, const scalarField &dydx0, const scalar dx, scalarField &y) const
Solve a single step dx and return the error.
Definition: Rosenbrock34.C:155
dict
dictionary dict
Definition: searchingEngine.H:14
Foam::dictionary
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:123
addToRunTimeSelectionTable.H
Macros for easy insertion into run-time selection tables.
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
Foam::constant::physicoChemical::c2
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Foam::LUBacksubstitute
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
Definition: scalarMatricesTemplates.C:120
Foam::LUDecompose
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
Definition: scalarMatrices.C:34
Foam::ODESystem
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:49
Foam::adaptiveSolver::resize
bool resize(const label n)
Resize the ODE solver.
Definition: adaptiveSolver.C:52
gamma
const scalar gamma
Definition: EEqn.H:9
x
x
Definition: LISASMDCalcMethod2.H:52
Rosenbrock34.H
Foam::defineTypeNameAndDebug
defineTypeNameAndDebug(combustionModel, 0)
Foam::ODESolver::resize
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:92
Foam::adaptiveSolver
Definition: adaptiveSolver.H:53
y
scalar y
Definition: LISASMDCalcMethod1.H:14