SR1.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 "SR1.H"
32
33// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34
35namespace Foam
36{
39 (
41 SR1,
43 );
44}
45
46
47// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
48
50{
51 // Set active design variables, if necessary
53 {
55 }
56
57 // Set previous HessianInv to be a diagonal matrix
60 {
61 temp[i][i] = scalar(1);
62 }
63
64 // Allocate correct size and content to HessianInv matrices
65 // has a max. capability of approximately 34000 variables.
66 HessianInvOld_ = temp;
67 HessianInv_ = temp;
68}
69
70
72{
73 // Vectors needed to construct the inverse HessianInv matrix
74 scalarField y(activeDesignVars_.size(), Zero);
75 scalarField s(activeDesignVars_.size(), Zero);
76 y.map(objectiveDerivatives_ - derivativesOld_, activeDesignVars_);
77 s.map(correctionOld_, activeDesignVars_);
78
79 scalarField temp(s - rightMult(HessianInvOld_, y));
80
81 // Construct the inverse HessianInv
82 scalar tempMag = sqrt(globalSum(sqr(temp)));
83 scalar yMag = sqrt(globalSum(sqr(y)));
84 scalar HessYMag = sqrt(globalSum(sqr(rightMult(HessianInvOld_, y))));
85
86 // Stability check
87 if (tempMag > ratioThreshold_ * yMag * HessYMag)
88 {
89 HessianInv_ =
90 HessianInvOld_
91 + (scalar(1)/(globalSum(temp*y)))*outerProd(temp, temp);
92 }
93 else
94 {
96 << "Denominator of update too small. Keeping old Hessian" << endl;
97 HessianInv_ = HessianInvOld_;
98 }
99}
100
101
103{
104 // In the first few iterations, use steepest descent but update the Hessian
105 // matrix
106 if (counter_ < nSteepestDescent_)
107 {
108 Info<< "Using steepest descent to update design variables ... " << endl;
109 correction_ = -eta_*objectiveDerivatives_;
110 }
111 else
112 {
113 scalarField activeDerivs(activeDesignVars_.size(), Zero);
114 activeDerivs.map(objectiveDerivatives_, activeDesignVars_);
115 scalarField activeCorrection
116 (
117 -etaHessian_*rightMult(HessianInv_, activeDerivs)
118 );
119
120 // Transfer correction to the global list
121 correction_ = Zero;
122 forAll(activeDesignVars_, varI)
123 {
124 correction_[activeDesignVars_[varI]] = activeCorrection[varI];
125 }
126 }
127
128 // Store fields for the next iteration
129 derivativesOld_ = objectiveDerivatives_;
130 correctionOld_ = correction_;
131 HessianInvOld_ = HessianInv_;
132}
133
134
136{
137 if (optMethodIODict_.headerOk())
138 {
139 optMethodIODict_.readEntry("HessianInvOld", HessianInvOld_);
140 optMethodIODict_.readEntry("derivativesOld", derivativesOld_);
141 optMethodIODict_.readEntry("correctionOld", correctionOld_);
142 optMethodIODict_.readEntry("counter", counter_);
143 optMethodIODict_.readEntry("eta", eta_);
144
145 const label n(HessianInvOld_.n());
146 HessianInv_ = SquareMatrix<scalar>(n, Zero);
147 correction_ = scalarField(correctionOld_.size(), Zero);
148
149 if (activeDesignVars_.empty())
150 {
151 activeDesignVars_ = identity(n);
152 }
153 }
154}
155
156
157// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
158
160:
162
163 // Construct null matrix since we dont know the dimension yet
164 etaHessian_
165 (
166 coeffsDict().getOrDefault<scalar>("etaHessian", 1)
167 ),
168 nSteepestDescent_
169 (
170 coeffsDict().getOrDefault<label>("nSteepestDescent", 1)
171 ),
172 ratioThreshold_
173 (
174 coeffsDict().getOrDefault<scalar>("ratioThreshold", 1e-08)
175 ),
176 activeDesignVars_(0),
177 HessianInv_(),
178 HessianInvOld_(),
179 derivativesOld_(0),
180 correctionOld_(0),
181 counter_(0)
182{
183 if
184 (
185 !coeffsDict().readIfPresent("activeDesignVariables", activeDesignVars_)
186 )
187 {
188 // If not, all available design variables will be used. Number is not
189 // know at the moment
190 Info<< "\t Didn't find explicit definition of active design variables. "
191 << "Treating all available ones as active " << endl;
192 }
193
194 // Read old hessian, correction and derivatives, if present
195 readFromDict();
196}
197
198
199// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
200
202{
203 if (counter_ == 0)
204 {
205 allocateMatrices();
206 }
207 else
208 {
209 updateHessian();
210 }
211
212 update();
213 ++counter_;
214}
215
216
218{
220 correctionOld_ = oldCorrection;
221}
222
223
225{
226 optMethodIODict_.add<SquareMatrix<scalar>>
227 (
228 "HessianInvOld",
229 HessianInvOld_,
230 true
231 );
232 optMethodIODict_.add<scalarField>("derivativesOld", derivativesOld_, true);
233 optMethodIODict_.add<scalarField>("correctionOld", correctionOld_, true);
234 optMethodIODict_.add<label>("counter", counter_, true);
235
237}
238
239
240// ************************************************************************* //
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.
void map(const UList< Type > &mapF, const labelUList &mapAddressing)
1 to 1 map from the given field
Definition: Field.C:240
The quasi-Newton Symmetric Rank One formula.
Definition: SR1.H:57
void computeCorrection()
Compute design variables correction.
Definition: SR1.C:201
SquareMatrix< scalar > HessianInvOld_
The previous Hessian inverse.
Definition: SR1.H:79
void readFromDict()
Read old info from dict.
Definition: SR1.C:135
labelList activeDesignVars_
Map to active design variables.
Definition: SR1.H:72
void allocateMatrices()
Allocate matrices in the first optimisation cycle.
Definition: SR1.C:49
virtual void write()
Write old info to dict.
Definition: SR1.C:224
SquareMatrix< scalar > HessianInv_
Definition: SR1.H:76
virtual void updateOldCorrection(const scalarField &oldCorrection)
Definition: SR1.C:217
void update()
Update design variables.
Definition: SR1.C:102
void updateHessian()
Update approximation of the inverse Hessian.
Definition: SR1.C:71
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
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 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
dimensionedSymmTensor sqr(const dimensionedVector &dv)
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
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)