BFGS.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-2023 PCOpt/NTUA
9  Copyright (C) 2013-2023 FOSS GP
10  Copyright (C) 2019-2020 OpenCFD Ltd.
11 -------------------------------------------------------------------------------
12 License
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 "BFGS.H"
32 
33 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37  defineTypeNameAndDebug(BFGS, 0);
39  (
40  updateMethod,
41  BFGS,
43  );
44 }
45 
46 
47 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
48 
50 {
51  // Vectors needed to construct the inverse HessianInv matrix
56 
57  scalar ys = globalSum(s*y);
58  if (counter_ == 1 && scaleFirstHessian_)
59  {
60  if (ys > scalar(0))
61  {
62  scalar scaleFactor = ys/globalSum(y*y);
63  Info<< "Scaling Hessian with factor " << scaleFactor << endl;
65  {
66  Hessian_()[varI][varI] *= scaleFactor;
67  }
68  }
69  else
70  {
72  << "y*s is negative. Skipping the scaling of the first Hessian"
73  << endl;
74  }
75  }
76 
77  // Check curvature condition
78  if (ys < curvatureThreshold_)
79  {
81  << " y*s is below threshold! y*s=" << ys << endl;
82  }
83 
84  DebugInfo
85  << "Hessian curvature index " << ys << endl;
86 
87  // Construct the inverse HessianInv
88  Hessian_() +=
89  (ys + globalSum(leftMult(y, Hessian_())*y))/sqr(ys)*outerProd(s, s)
90  - (scalar(1)/ys)*
91  (
93  + outerProd(s, leftMult(y, Hessian_()))
94  );
95 }
96 
97 
98 void Foam::BFGS::update()
99 {
100  // In the first few iterations, use steepest descent but update the Hessian
101  // matrix
102  if (counter_ < nSteepestDescent_)
103  {
104  Info<< "Using steepest descent to update design variables" << endl;
105  for (const label varI : activeDesignVars_)
106  {
107  correction_[varI] = -eta_*objectiveDerivatives_[varI];
108  }
109  }
110  // Else use BFGS formula to update the design variables
111  else
112  {
113  // Compute correction for active design variables
114  scalarField activeDerivs(activeDesignVars_.size(), Zero);
115  activeDerivs.map(objectiveDerivatives_, activeDesignVars_);
116  scalarField activeCorrection
117  (
118  -etaHessian_*rightMult(Hessian_(), activeDerivs)
119  );
120 
121  // Transfer correction to the global list
122  correction_ = Zero;
123  forAll(activeDesignVars_, varI)
124  {
125  correction_[activeDesignVars_[varI]] = activeCorrection[varI];
126  }
127  }
128 
129  // Store fields for the next iteration
130  derivativesOld_ = objectiveDerivatives_;
131  correctionOld_ = correction_;
132 }
133 
134 
135 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
136 
137 Foam::BFGS::BFGS
138 (
139  const fvMesh& mesh,
140  const dictionary& dict,
141  autoPtr<designVariables>& designVars,
142  const label nConstraints,
143  const word& type
144 )
145 :
146  quasiNewton(mesh, dict, designVars, nConstraints, type),
147  curvatureThreshold_
148  (
149  coeffsDict(type).getOrDefault<scalar>("curvatureThreshold", 1e-10)
150  )
151 {
152  allocateHessian();
153 }
154 
155 
156 // ************************************************************************* //
dictionary dict
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:116
A list of keyword definitions, which are a keyword followed by a number of values (eg...
Definition: dictionary.H:129
scalar curvatureThreshold_
Curvature threshold.
Definition: BFGS.H:76
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:531
bool scaleFirstHessian_
Scale the initial unitary Hessian approximation.
Definition: quasiNewton.H:70
const scalarField rightMult(const SquareMatrix< scalar > &, const scalarField &)
Definition: updateMethod.C:65
const scalarField leftMult(const scalarField &, const SquareMatrix< scalar > &)
Definition: updateMethod.C:38
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:421
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
Definition: POSIX.C:799
scalar y
const labelList & activeDesignVars_
Map to active design variables.
Definition: updateMethod.H:75
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
dynamicFvMesh & mesh
Base class for quasi-Newton methods.
Definition: quasiNewton.H:49
A class for handling words, derived from Foam::string.
Definition: word.H:63
autoPtr< SquareMatrix< scalar > > Hessian_
The Hessian or its inverse, depending on the deriving class.
Definition: quasiNewton.H:78
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
label counter_
Optimisation cycle count.
Definition: updateMethod.H:123
scalarField derivativesOld_
The previous derivatives.
Definition: quasiNewton.H:83
virtual void updateHessian()
Update approximation of the inverse Hessian.
Definition: BFGS.C:42
#define DebugInfo
Report an information message using Foam::Info.
void map(const UList< Type > &mapF, const labelUList &mapAddressing)
1 to 1 map from the given field
Definition: Field.C:303
scalar globalSum(const scalarField &field)
Compute either global or local sum, based on globalSum flag.
Definition: updateMethod.C:168
void allocateHessian()
Allocate the Hessian matrix.
Definition: quasiNewton.C:35
defineTypeNameAndDebug(combustionModel, 0)
virtual void update()
Update design variables.
Definition: BFGS.C:91
scalarField objectiveDerivatives_
Derivatives of the objective functions.
Definition: updateMethod.H:80
#define WarningInFunction
Report a warning using Foam::Warning.
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
SquareMatrix< scalar > outerProd(const scalarField &, const scalarField &)
Definition: updateMethod.C:92
messageStream Info
Information stream (stdout output on master, null elsewhere)
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
Definition: HashPtrTable.H:48
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))
Namespace for OpenFOAM.
addToRunTimeSelectionTable(functionObject, pointHistory, dictionary)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:127
scalarField correctionOld_
The previous correction.
Definition: quasiNewton.H:88