ODESystem.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) 2021 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 "ODESystem.H"
27 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
31 {}
32 
33 
34 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
35 
37 {}
38 
39 
40 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
41 
43 (
44  const scalar x,
45  const scalarField& y,
46  const scalarField& dy,
47  const label li
48 ) const
49 {
50  // Evaluate the derivatives using the derivatives method
51  scalarField dfdx0(nEqns());
52  derivatives(x, y, li, dfdx0);
53 
54  // Evaluate the derivatives and the Jacobian using the Jacobian method
55  scalarField dfdx1(nEqns());
56  scalarSquareMatrix d2fdxdyAnalytic(nEqns());
57  jacobian(x, y, li, dfdx1, d2fdxdyAnalytic);
58 
59  // Compare derivatives
60  Info<< "[derivatives] dfdx = ( ";
61  forAll(dfdx0, i) { Info<< dfdx0[i] << ' '; }
62  Info<< ")" << nl;
63  Info<< "[ jacobian] dfdx = ( ";
64  forAll(dfdx1, i) { Info<< dfdx1[i] << ' '; }
65  Info<< ")" << nl;
66  Info<< "[ ratio] dfdx = ( ";
67  forAll(dfdx1, i) { Info<< dfdx1[i]/stabilise(dfdx0[i], rootVSmall) << ' '; }
68  Info<< ")" << nl;
69 
70  // Construct a Jacobian using the finite differences and the derivatives
71  // method
72  scalarSquareMatrix d2fdxdyFiniteDifference(nEqns());
73  scalarField y0(y), y1(y);
74  for (label i = 0; i < nEqns(); ++ i)
75  {
76  y0[i] = y[i] - dy[i];
77  y1[i] = y[i] + dy[i];
78 
79  derivatives(x, y0, li, dfdx0);
80  derivatives(x, y1, li, dfdx1);
81 
82  for (label j = 0; j < nEqns(); j++)
83  {
84  d2fdxdyFiniteDifference(j, i) = (dfdx1[j] - dfdx0[j])/(2*dy[i]);
85  }
86 
87  y0[i] = y[i];
88  y1[i] = y[i];
89  }
90 
91  for (label i = 0; i < nEqns(); ++ i)
92  {
93  UList<scalar> FD(d2fdxdyFiniteDifference[i], nEqns());
94  Info<< "[derivatives] d2fdxdy[" << i << "] = ( ";
95  forAll(FD, i) { Info<< FD[i] << ' '; }
96  Info<< ")" << nl;
97  }
98  for (label i = 0; i < nEqns(); ++ i)
99  {
100  UList<scalar> A(d2fdxdyAnalytic[i], nEqns());
101  Info<< "[ jacobian] d2fdxdy[" << i << "] = ( ";
102  forAll(A, i) { Info<< A[i] << ' '; }
103  Info<< ")" << nl;
104  }
105  for (label i = 0; i < nEqns(); ++ i)
106  {
107  UList<scalar> FD(d2fdxdyFiniteDifference[i], nEqns());
108  UList<scalar> A(d2fdxdyAnalytic[i], nEqns());
109  Info<< "[ ratio] d2fdxdy[" << i << "] = ( ";
110  forAll(A, i) { Info<< A[i]/stabilise(FD[i], rootVSmall) << ' '; }
111  Info<< ")" << nl;
112  }
113 }
114 
115 
116 // ************************************************************************* //
static const Foam::dimensionedScalar A("A", Foam::dimPressure, 611.21)
scalar y
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
ODESystem()
Construct null.
Definition: ODESystem.C:30
virtual ~ODESystem()
Destructor.
Definition: ODESystem.C:36
void check(const scalar x, const scalarField &y, const scalarField &dy, const label li) const
Print correspondence between derivatives and Jacobian methods.
Definition: ODESystem.C:43
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:74
tmp< DimensionedField< scalar, GeoMesh > > stabilise(const DimensionedField< scalar, GeoMesh > &dsf, const dimensioned< scalar > &ds)
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
dimensionedScalar y0(const dimensionedScalar &ds)
messageStream Info
dimensionedScalar y1(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:266