72 y.map(objectiveDerivatives_ - derivativesOld_, activeDesignVars_);
73 s.map(correctionOld_, activeDesignVars_);
75 scalar ys = globalSum(
s*
y);
77 if (counter_ == 1 && scaleFirstHessian_)
79 scalar scaleFactor = ys/globalSum(
y*
y);
80 Info<<
"Scaling Hessian with factor " << scaleFactor <<
endl;
81 forAll(activeDesignVars_, varI)
83 HessianInvOld_[varI][varI] *= scaleFactor;
88 if (ys < curvatureThreshold_)
91 <<
" y*s is below threshold! y*s=" << ys <<
endl;
95 <<
"Hessian curvature index " << ys <<
endl;
100 + (ys + globalSum(leftMult(
y, HessianInvOld_)*
y))/
sqr(ys)*outerProd(
s,
s)
103 outerProd(rightMult(HessianInvOld_,
y),
s)
104 + outerProd(
s, leftMult(
y, HessianInvOld_))
113 if (counter_ < nSteepestDescent_)
115 Info<<
"Using steepest descent to update design variables" <<
endl;
116 correction_ = -eta_*objectiveDerivatives_;
123 activeDerivs.
map(objectiveDerivatives_, activeDesignVars_);
126 -etaHessian_*rightMult(HessianInv_, activeDerivs)
131 forAll(activeDesignVars_, varI)
133 correction_[activeDesignVars_[varI]] = activeCorrection[varI];
138 derivativesOld_ = objectiveDerivatives_;
139 correctionOld_ = correction_;
140 HessianInvOld_ = HessianInv_;
146 if (optMethodIODict_.headerOk())
148 optMethodIODict_.readEntry(
"HessianInvOld", HessianInvOld_);
149 optMethodIODict_.readEntry(
"derivativesOld", derivativesOld_);
150 optMethodIODict_.readEntry(
"correctionOld", correctionOld_);
151 optMethodIODict_.readEntry(
"counter", counter_);
152 optMethodIODict_.readEntry(
"eta", eta_);
154 const label
n(HessianInvOld_.n());
158 if (activeDesignVars_.empty())
171 const dictionary&
dict 177 coeffsDict().getOrDefault<scalar>(
"etaHessian", 1)
181 coeffsDict().getOrDefault<label>(
"nSteepestDescent", 1)
183 activeDesignVars_(0),
186 coeffsDict().getOrDefault<bool>(
"scaleFirstHessian", false)
190 coeffsDict().getOrDefault<scalar>(
"curvatureThreshold", 1
e-10)
206 Info<<
"\t Did not find explicit definition of active design variables." 207 <<
" Treating all available ones as active" <<
endl;
236 correctionOld_ = oldCorrection;
248 optMethodIODict_.add<
scalarField>(
"derivativesOld", derivativesOld_,
true);
249 optMethodIODict_.add<
scalarField>(
"correctionOld", correctionOld_,
true);
250 optMethodIODict_.add<label>(
"counter", counter_,
true);
void size(const label n)
Older name for setAddressableSize.
SquareMatrix< scalar > HessianInvOld_
The previous Hessian inverse.
A list of keyword definitions, which are a keyword followed by a number of values (eg...
void allocateMatrices()
Allocate matrices in the first optimisation cycle.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
bool empty() const noexcept
True if List is empty (ie, size() is zero)
Ostream & endl(Ostream &os)
Add newline and flush stream.
virtual void write()
Write old info to dict.
Macros for easy insertion into run-time selection tables.
#define forAll(list, i)
Loop across all elements in list.
virtual void updateOldCorrection(const scalarField &oldCorrection)
Update old correction. useful for quasi-newton methods coupled with line search.
const dimensionedScalar e
Elementary charge.
static const Identity< scalar > I
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
void readFromDict()
Read old info from dict.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void updateHessian()
Update approximation of the inverse Hessian.
#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
defineTypeNameAndDebug(combustionModel, 0)
void update()
Update design variables.
addToRunTimeSelectionTable(decompositionMethod, kahipDecomp, dictionary)
SquareMatrix< scalar > HessianInv_
The Hessian inverse. Should have the size of the active design variables.
scalarField objectiveDerivatives_
Derivatives of the objective functions.
#define WarningInFunction
Report a warning using Foam::Warning.
meshDefDict readIfPresent("polyMeshPatches", polyPatchNames)
messageStream Info
Information stream (stdout output on master, null elsewhere)
virtual void write()
Write useful quantities to files.
labelList activeDesignVars_
Map to active design variables.
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))
virtual void updateOldCorrection(const scalarField &oldCorrection)
Update old correction. Useful for quasi-Newton methods coupled with line search.
dictionary coeffsDict()
Return optional dictionary with parameters specific to each method.
void computeCorrection()
Compute design variables correction.
static constexpr const zero Zero
Global zero (0)