DBFGS.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) 2007-2019 PCOpt/NTUA
9 Copyright (C) 2013-2019 FOSS GP
10 Copyright (C) 2019-2020 OpenCFD Ltd.
11-------------------------------------------------------------------------------
12License
13 This file is part of OpenFOAM.
14
15 OpenFOAM is free software: you can redistribute it and/or modify it
16 under the terms of the GNU General Public License as published by
17 the Free Software Foundation, either version 3 of the License, or
18 (at your option) any later version.
19
20 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
23 for more details.
24
25 You should have received a copy of the GNU General Public License
26 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
27
28\*---------------------------------------------------------------------------*/
29
30#include "DBFGS.H"
32
33// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34
35namespace Foam
36{
39 (
41 DBFGS,
43 );
44}
45
46
47// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
48
50{
51 // Set active design variables, if necessary
53 {
55 }
56
57 // Set previous Hessian to be a diagonal matrix
59
60 // Allocate correct size and content to Hessian matrices
61 // has a max. capability of approximately 34000 variables.
62 HessianOld_ = temp;
63 Hessian_ = temp;
64}
65
66
68{
69 // Vectors needed to construct the inverse Hessian matrix
70 scalarField y(activeDesignVars_.size(), Zero);
71 scalarField s(activeDesignVars_.size(), Zero);
72 y.map(objectiveDerivatives_ - derivativesOld_, activeDesignVars_);
73 s.map(correctionOld_, activeDesignVars_);
74
75 scalar ys = globalSum(s*y);
76
77 if (counter_ == 1 && scaleFirstHessian_)
78 {
79 scalar scaleFactor = ys/globalSum(y*y);
80 Info<< "Scaling Hessian with factor " << scaleFactor << endl;
81 forAll(activeDesignVars_, varI)
82 {
83 HessianOld_[varI][varI] /= scaleFactor;
84 }
85 }
86
87 scalar sBs = globalSum(leftMult(s, HessianOld_)*s);
88
89 // Check curvature condition and apply dampening is necessary
90 scalar theta(1);
91 if (ys < gamma_*sBs)
92 {
94 << " y*s is below threshold. Using damped form" << endl;
95 theta = (scalar(1)-gamma_)*sBs/(sBs - ys);
96 }
97
99 << "Hessian curvature index " << ys << endl;
100
101 scalarField r(theta*y + (scalar(1)-theta)*rightMult(HessianOld_, s));
102
103 // Construct the inverse Hessian
104 Hessian_ =
105 HessianOld_
106 - outerProd(rightMult(HessianOld_, s), leftMult(s/sBs, HessianOld_))
107 + outerProd(r, r/globalSum(s*r));
108}
109
110
112{
113 SquareMatrix<scalar> HessianInv = inv(Hessian_);
114
115 // In the first few iterations, use steepest descent but update the Hessian
116 // matrix
117 if (counter_ < nSteepestDescent_)
118 {
119 Info<< "Using steepest descent to update design variables" << endl;
120 correction_ = -eta_*objectiveDerivatives_;
121 }
122 // Else use DBFGS formula to update design variables
123 else
124 {
125 // Compute correction for active design variables
126 scalarField activeDerivs(activeDesignVars_.size(), Zero);
127 activeDerivs.map(objectiveDerivatives_, activeDesignVars_);
128 scalarField activeCorrection
129 (
130 -etaHessian_*rightMult(HessianInv, activeDerivs)
131 );
132
133 // Transfer correction to the global list
134 correction_ = Zero;
135 forAll(activeDesignVars_, varI)
136 {
137 correction_[activeDesignVars_[varI]] = activeCorrection[varI];
138 }
139 }
140
141 // Store fields for the next iteration
142 derivativesOld_ = objectiveDerivatives_;
143 correctionOld_ = correction_;
144 HessianOld_ = Hessian_;
145}
146
147
149{
150 if (optMethodIODict_.headerOk())
151 {
152 optMethodIODict_.readEntry("HessianOld", HessianOld_);
153 optMethodIODict_.readEntry("derivativesOld", derivativesOld_);
154 optMethodIODict_.readEntry("correctionOld", correctionOld_);
155 optMethodIODict_.readEntry("counter", counter_);
156 optMethodIODict_.readEntry("eta", eta_);
157
158 label n = HessianOld_.n();
159 Hessian_ = SquareMatrix<scalar>(n, Zero);
160 correction_ = scalarField(correctionOld_.size(), Zero);
161
162 if (activeDesignVars_.empty())
163 {
164 activeDesignVars_ = identity(n);
165 }
166 }
167}
168
169
170// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
171
173(
174 const fvMesh& mesh,
175 const dictionary& dict
176)
177:
179
180 // Construct null matrix since we dont know the dimension yet
181 etaHessian_
182 (
183 coeffsDict().getOrDefault<scalar>("etaHessian", 1)
184 ),
185 nSteepestDescent_
186 (
187 coeffsDict().getOrDefault<label>("nSteepestDescent", 1)
188 ),
189 activeDesignVars_(0),
190 scaleFirstHessian_
191 (
192 coeffsDict().getOrDefault<bool>("scaleFirstHessian", false)
193 ),
194 curvatureThreshold_
195 (
196 coeffsDict().getOrDefault<scalar>("curvatureThreshold", 1e-10)
197 ),
198 Hessian_(),
199 HessianOld_(),
200 derivativesOld_(0),
201 correctionOld_(0),
202 counter_(0),
203 gamma_(coeffsDict().getOrDefault<scalar>("gamma", 0.2))
204
205{
206 if
207 (
208 !coeffsDict().readIfPresent("activeDesignVariables", activeDesignVars_)
209 )
210 {
211 // If not, all available design variables will be used. Number is not
212 // know at the moment
213 Info<< "\t Did not find explicit definition of active design variables."
214 " Treating all available ones as active " << endl;
215 }
216
217 // read old hessian, correction and derivatives, if present
218 readFromDict();
219}
220
221
222// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
223
225{
226 if (counter_ == 0)
227 {
228 allocateMatrices();
229 }
230 else
231 {
232 updateHessian();
233 }
234
235 update();
236 ++counter_;
237}
238
239
241{
243 correctionOld_ = oldCorrection;
244}
245
246
248{
249 optMethodIODict_.add<SquareMatrix<scalar>>("HessianOld", HessianOld_, true);
250 optMethodIODict_.add<scalarField>("derivativesOld", derivativesOld_, true);
251 optMethodIODict_.add<scalarField>("correctionOld", correctionOld_, true);
252 optMethodIODict_.add<label>("counter", counter_, true);
253
255}
256
257
258// ************************************************************************* //
scalar y
label n
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
The quasi-Newton BFGS formula with the dampening proposed by Powell.
Definition: DBFGS.H:57
void computeCorrection()
Compute design variables correction.
Definition: DBFGS.C:224
SquareMatrix< scalar > HessianOld_
The previous Hessian.
Definition: DBFGS.H:81
void readFromDict()
Read old info from dict.
Definition: DBFGS.C:148
labelList activeDesignVars_
Map to active design variables.
Definition: DBFGS.H:69
void allocateMatrices()
Allocate matrices in the first optimisation cycle.
Definition: DBFGS.C:49
SquareMatrix< scalar > Hessian_
The Hessian. Should have the size of the active design variables.
Definition: DBFGS.H:78
virtual void write()
Write old info to dict.
Definition: DBFGS.C:247
virtual void updateOldCorrection(const scalarField &oldCorrection)
Definition: DBFGS.C:240
void update()
Update design variables.
Definition: DBFGS.C:111
void updateHessian()
Update approximation of the inverse Hessian.
Definition: DBFGS.C:67
void map(const UList< Type > &mapF, const labelUList &mapAddressing)
1 to 1 map from the given field
Definition: Field.C:240
A templated (N x N) square matrix of objects of <Type>, containing N*N elements, derived from Matrix.
Definition: SquareMatrix.H:66
bool empty() const noexcept
True if the UList is empty (ie, size() is zero)
Definition: UListI.H:427
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Abstract base class for optimisation methods.
Definition: updateMethod.H:55
scalarField objectiveDerivatives_
Derivatives of the objective functions.
Definition: updateMethod.H:68
dictionary coeffsDict()
Return optional dictionary with parameters specific to each method.
Definition: updateMethod.C:254
virtual void write()
Write useful quantities to files.
Definition: updateMethod.C:398
virtual void updateOldCorrection(const scalarField &oldCorrection)
Definition: updateMethod.C:390
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
bool
Definition: EEqn.H:20
mesh update()
dynamicFvMesh & mesh
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define DebugInfo
Report an information message using Foam::Info.
#define WarningInFunction
Report a warning using Foam::Warning.
Namespace for OpenFOAM.
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
messageStream Info
Information stream (stdout output on master, null elsewhere)
static const Identity< scalar > I
Definition: Identity.H:94
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
dictionary dict
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
propsDict readIfPresent("fields", acceptFields)