rodas23.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-------------------------------------------------------------------------------
11License
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 "rodas23.H"
31
32// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33
34namespace Foam
35{
38
39const scalar
40 rodas23::c3 = 1,
41 rodas23::d1 = 1.0/2.0,
42 rodas23::d2 = 3.0/2.0,
43 rodas23::a31 = 2,
44 rodas23::a41 = 2,
45 rodas23::c21 = 4,
46 rodas23::c31 = 1,
47 rodas23::c32 = -1,
48 rodas23::c41 = 1,
49 rodas23::c42 = -1,
50 rodas23::c43 = -8.0/3.0,
51 rodas23::gamma = 1.0/2.0;
52}
53
54
55// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
56
58:
61 k1_(n_),
62 k2_(n_),
63 k3_(n_),
64 dy_(n_),
65 err_(n_),
66 dydx_(n_),
67 dfdx_(n_),
68 dfdy_(n_, n_),
69 a_(n_, n_),
70 pivotIndices_(n_)
71{}
72
73
74// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
75
77{
79 {
81
82 resizeField(k1_);
83 resizeField(k2_);
84 resizeField(k3_);
85 resizeField(dy_);
86 resizeField(err_);
87 resizeField(dydx_);
88 resizeField(dfdx_);
89 resizeMatrix(dfdy_);
90 resizeMatrix(a_);
91 resizeField(pivotIndices_);
92
93 return true;
94 }
95
96 return false;
97}
98
99
101(
102 const scalar x0,
103 const scalarField& y0,
104 const scalarField& dydx0,
105 const scalar dx,
107) const
108{
109 odes_.jacobian(x0, y0, dfdx_, dfdy_);
110
111 for (label i=0; i<n_; i++)
112 {
113 for (label j=0; j<n_; j++)
114 {
115 a_(i, j) = -dfdy_(i, j);
116 }
117
118 a_(i, i) += 1.0/(gamma*dx);
119 }
120
121 LUDecompose(a_, pivotIndices_);
122
123 // Calculate k1:
124 forAll(k1_, i)
125 {
126 k1_[i] = dydx0[i] + dx*d1*dfdx_[i];
127 }
128
129 LUBacksubstitute(a_, pivotIndices_, k1_);
130
131 // Calculate k2:
132 forAll(k2_, i)
133 {
134 k2_[i] = dydx0[i] + dx*d2*dfdx_[i] + c21*k1_[i]/dx;
135 }
136
137 LUBacksubstitute(a_, pivotIndices_, k2_);
138
139 // Calculate k3:
140 forAll(y, i)
141 {
142 dy_[i] = a31*k1_[i];
143 y[i] = y0[i] + dy_[i];
144 }
145
146 odes_.derivatives(x0 + dx, y, dydx_);
147
148 forAll(k3_, i)
149 {
150 k3_[i] = dydx_[i] + (c31*k1_[i] + c32*k2_[i])/dx;
151 }
152
153 LUBacksubstitute(a_, pivotIndices_, k3_);
154
155 // Calculate new state and error
156 forAll(y, i)
157 {
158 dy_[i] += k3_[i];
159 y[i] = y0[i] + dy_[i];
160 }
161
162 odes_.derivatives(x0 + dx, y, dydx_);
163
164 forAll(err_, i)
165 {
166 err_[i] = dydx_[i] + (c41*k1_[i] + c42*k2_[i] + c43*k3_[i])/dx;
167 }
168
169 LUBacksubstitute(a_, pivotIndices_, err_);
170
171 forAll(y, i)
172 {
173 y[i] = y0[i] + dy_[i] + err_[i];
174 }
175
176 return normalizeError(y0, y, err_);
177}
178
179
181(
182 scalar& x,
183 scalarField& y,
184 scalar& dxTry
185) const
186{
187 adaptiveSolver::solve(odes_, x, y, dxTry);
188}
189
190
191// ************************************************************************* //
scalar y
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
virtual bool resize()
Resize the ODE solver.
Definition: Euler.C:53
Abstract base-class for ODE system solvers.
Definition: ODESolver.H:57
virtual bool resize()=0
Resize the ODE solver.
Definition: ODESolver.C:92
Abstract base class for the systems of ordinary differential equations.
Definition: ODESystem.H:50
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
An ODE solver for chemistry.
Definition: ode.H:55
L-stable, stiffly-accurate embedded Rosenbrock ODE solver of order (2)3.
Definition: rodas23.H:68
virtual bool resize()
Resize the ODE solver.
Definition: rodas23.C:76
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const scalar gamma
Definition: EEqn.H:9
Namespace for OpenFOAM.
dimensionedScalar y0(const dimensionedScalar &ds)
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
dictionary dict
CEqn solve()
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333