mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
ENH: treating issues with the convergence of ISQP
The solution of the QP subproblem can become quite expensive, especially for cases with many design variables (e.g. topology optimisation). A (potentially dense) matrix with the size of the design variables is solved using a matrix-free CG solver. The convergence speed greatly depends on the used preconditioner. This commit adds preconditioner-vector products based on the L-BFGS inverse Hessian and, more importantly, a preconditioner computed using the Sherman-Morrison formula. The latter is applicable here since the LHS of the QP problem is computed as the sum of rank-2 L-BFGS updates, a sum of rank-1 updates (as many as the flow-related constraints) and a diagonal matrix depending on the bound constraints. Additionally, the QP subproblem could have no feasible points. To relax this, constraints can be applied gradually through the targetConstraintReduction enty (typical value of 0.1 for topology optimisation).
This commit is contained in:
Vaggelis Papoutsis
committed by
Andrew Heather
Andrew Heather
parent
0120d7720b
commit
0af46becdd
@ -5,8 +5,8 @@
|
||||
\\ / A nd | www.openfoam.com
|
||||
\\/ M anipulation |
|
||||
-------------------------------------------------------------------------------
|
||||
Copyright (C) 2020-2021 PCOpt/NTUA
|
||||
Copyright (C) 2020-2021 FOSS GP
|
||||
Copyright (C) 2020-2023 PCOpt/NTUA
|
||||
Copyright (C) 2020-2023 FOSS GP
|
||||
-------------------------------------------------------------------------------
|
||||
License
|
||||
This file is part of OpenFOAM.
|
||||
@ -34,7 +34,7 @@ License
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
defineTypeNameAndDebug(ISQP, 0);
|
||||
defineTypeNameAndDebug(ISQP, 2);
|
||||
addToRunTimeSelectionTable
|
||||
(
|
||||
updateMethod,
|
||||
@ -50,6 +50,15 @@ namespace Foam
|
||||
}
|
||||
|
||||
|
||||
const Foam::Enum<Foam::ISQP::preconditioners>
|
||||
Foam::ISQP::preconditionerNames
|
||||
({
|
||||
{ preconditioners::diag, "diag" },
|
||||
{ preconditioners::invHessian, "invHessian" },
|
||||
{ preconditioners::ShermanMorrison, "ShermanMorrison" }
|
||||
});
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
void Foam::ISQP::updateSizes()
|
||||
@ -208,6 +217,7 @@ void Foam::ISQP::zeroUpdates()
|
||||
|
||||
void Foam::ISQP::solveDeltaPEqn()
|
||||
{
|
||||
addProfiling(ISQP, "ISQP::solveDeltaPEqn");
|
||||
// Explicit part of the right hand side of the deltaX equation
|
||||
scalarField FDx(-resFL());
|
||||
if (includeBoundConstraints_)
|
||||
@ -233,29 +243,6 @@ void Foam::ISQP::solveDeltaPEqn()
|
||||
}
|
||||
}
|
||||
CGforDeltaP(FDx);
|
||||
/*
|
||||
//Info<< "FDx::" << FDx << endl;
|
||||
//Info<< "invHFL::" << invHFL() << endl;
|
||||
|
||||
// Loop to obtain deltaX. Part of the LHS is treated explicitly, to
|
||||
// avoid solving a (potentially) dense system with the size of the design
|
||||
// variables
|
||||
scalarField rhs(computeRHSForDeltaX(FDx));
|
||||
scalar res(sqrt(globalSum(magSqr(deltaP_ - rhs))));
|
||||
scalar resInit(res);
|
||||
label iter(0);
|
||||
do
|
||||
{
|
||||
deltaP_ = rhs;
|
||||
rhs = computeRHSForDeltaX(FDx);
|
||||
res = sqrt(globalSum(magSqr(deltaP_ - rhs)));
|
||||
DebugInfo
|
||||
<< "Solving for deltaX, Initial Residual " << resInit
|
||||
<< ", Final Residual " << res << endl;
|
||||
resInit = res;
|
||||
} while (iter++ < maxDxIters_ && res > 1.e-07);
|
||||
// Info<< "deltaX solution " << deltaP_ << endl;
|
||||
// */
|
||||
}
|
||||
|
||||
|
||||
@ -302,9 +289,10 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::computeRHSForDeltaX
|
||||
|
||||
void Foam::ISQP::CGforDeltaP(const scalarField& FDx)
|
||||
{
|
||||
scalarField precond(deltaPDiagPreconditioner());
|
||||
scalarField r(FDx - DeltaPMatrixVectorProduct(deltaP_));
|
||||
scalarField z(precond*r);
|
||||
addProfiling(ISQP, "ISQP::CGforDeltaP");
|
||||
autoPtr<scalarField> precon(nullptr);
|
||||
scalarField r(FDx - matrixVectorProduct(deltaP_));
|
||||
scalarField z(preconVectorProduct(r, precon));
|
||||
scalarField p(z);
|
||||
scalar res(sqrt(globalSum(r*r)));
|
||||
scalar resInit(res);
|
||||
@ -313,17 +301,17 @@ void Foam::ISQP::CGforDeltaP(const scalarField& FDx)
|
||||
label iter(0);
|
||||
do
|
||||
{
|
||||
scalarField Ap(DeltaPMatrixVectorProduct(p));
|
||||
scalarField Ap(matrixVectorProduct(p));
|
||||
scalar a = rz/globalSum(p*Ap);
|
||||
deltaP_ += a*p;
|
||||
r -= a*Ap;
|
||||
res = sqrt(globalSum(r*r));
|
||||
z = precond*r;
|
||||
z = preconVectorProduct(r, precon);
|
||||
rz = globalSum(r*z);
|
||||
scalar beta = rz/rzOld;
|
||||
p = z + beta*p;
|
||||
rzOld = rz;
|
||||
} while (iter++ < maxDxIters_ && res > 1.e-09);
|
||||
} while (iter++ < maxDxIters_ && res > relTol_*resInit);
|
||||
DebugInfo
|
||||
<< "CG, Solving for deltaP, Initial Residual " << resInit
|
||||
<< ", Final Residual " << res
|
||||
@ -331,42 +319,13 @@ void Foam::ISQP::CGforDeltaP(const scalarField& FDx)
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::deltaPDiagPreconditioner()
|
||||
{
|
||||
tmp<scalarField> tpreconditioner(HessianDiag());
|
||||
//tmp<scalarField> tpreconditioner(SR1HessianDiag());
|
||||
scalarField& preconditioner = tpreconditioner.ref();
|
||||
|
||||
// Part related to the constraints
|
||||
forAll(constraintDerivatives_, cI)
|
||||
{
|
||||
scalarField cDerivs(constraintDerivatives_[cI], activeDesignVars_);
|
||||
scalar mult(gs_[cI]/lamdas_[cI]);
|
||||
if (includeExtraVars_)
|
||||
{
|
||||
mult += extraVars_()[cI]/z_()[cI];
|
||||
}
|
||||
preconditioner += sqr(cDerivs)/mult;
|
||||
}
|
||||
|
||||
if (includeBoundConstraints_)
|
||||
{
|
||||
preconditioner += lTilda_()/ls_() + uTilda_()/us_();
|
||||
}
|
||||
|
||||
preconditioner = 1./preconditioner;
|
||||
|
||||
return tpreconditioner;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::DeltaPMatrixVectorProduct
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::matrixVectorProduct
|
||||
(
|
||||
const scalarField& vector
|
||||
)
|
||||
{
|
||||
addProfiling(ISQP, "ISQP::MatrixVectorProduct");
|
||||
tmp<scalarField> tAp(HessianVectorProduct(vector));
|
||||
//tmp<scalarField> tAp(SR1HessianVectorProduct(vector));
|
||||
scalarField& Ap = tAp.ref();
|
||||
scalarField GsLAv(cValues_.size(), Zero);
|
||||
forAll(constraintDerivatives_, cI)
|
||||
@ -403,6 +362,134 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::DeltaPMatrixVectorProduct
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::preconVectorProduct
|
||||
(
|
||||
const scalarField& vector,
|
||||
autoPtr<scalarField>& precon
|
||||
)
|
||||
{
|
||||
addProfiling(ISQP, "ISQP::preconVectorProduct");
|
||||
if (preconType_ == preconditioners::diag)
|
||||
{
|
||||
if (!precon)
|
||||
{
|
||||
precon.reset(diagPreconditioner().ptr());
|
||||
}
|
||||
return precon()*vector;
|
||||
}
|
||||
else if (preconType_ == preconditioners::invHessian)
|
||||
{
|
||||
return invHessianVectorProduct(vector);
|
||||
}
|
||||
else if (preconType_ == preconditioners::ShermanMorrison)
|
||||
{
|
||||
return ShermanMorrisonPrecon(vector);
|
||||
}
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::diagPreconditioner()
|
||||
{
|
||||
addProfiling(ISQP, "ISQP::deltaPDiagPreconditioner");
|
||||
tmp<scalarField> tpreconditioner(HessianDiag());
|
||||
scalarField& preconditioner = tpreconditioner.ref();
|
||||
|
||||
// Part related to the constraints
|
||||
forAll(constraintDerivatives_, cI)
|
||||
{
|
||||
scalarField cDerivs(constraintDerivatives_[cI], activeDesignVars_);
|
||||
scalar mult(gs_[cI]/lamdas_[cI]);
|
||||
if (includeExtraVars_)
|
||||
{
|
||||
mult += extraVars_()[cI]/z_()[cI];
|
||||
}
|
||||
preconditioner += sqr(cDerivs)/mult;
|
||||
}
|
||||
|
||||
if (includeBoundConstraints_)
|
||||
{
|
||||
preconditioner += lTilda_()/ls_() + uTilda_()/us_();
|
||||
}
|
||||
|
||||
preconditioner = 1./preconditioner;
|
||||
|
||||
return tpreconditioner;
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::ShermanMorrisonPrecon
|
||||
(
|
||||
const scalarField& vector
|
||||
)
|
||||
{
|
||||
// Recursively update the inv(LHS)*vector since the LHS consists of the
|
||||
// L-BFGS-based Hessian, computed with rank-2 updates, and the part related
|
||||
// to flow constraints, computed as rank-1 updates. In the inversion
|
||||
// process, the diagonal matrix related to bound constraints is treated as
|
||||
// the initial matrix of the L-BFGS update.
|
||||
|
||||
// Constribution from bound constraints, treated as the seed of the
|
||||
// L-BFGS inverse
|
||||
tmp<scalarField> tdiag(nullptr);
|
||||
if (includeBoundConstraints_)
|
||||
{
|
||||
tdiag.reset(lTilda_()/ls_() + uTilda_()/us_());
|
||||
}
|
||||
|
||||
// List of vectors to be used in the rank-1 updates related to the flow
|
||||
// constraitns
|
||||
PtrList<scalarField> r1Updates(cValues_.size());
|
||||
|
||||
forAll(constraintDerivatives_, cI)
|
||||
{
|
||||
const scalarField& cDerivsI = constraintDerivatives_[cI];
|
||||
r1Updates.set(cI, new scalarField(cDerivsI, activeDesignVars_));
|
||||
}
|
||||
// Multipliers of the rank-1 updates
|
||||
scalarField mult(gs_/lamdas_);
|
||||
if (includeExtraVars_)
|
||||
{
|
||||
mult += extraVars_()/z_();
|
||||
}
|
||||
|
||||
return
|
||||
ShermanMorrisonRank1Update(r1Updates, vector, tdiag, mult, mult.size());
|
||||
}
|
||||
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::ShermanMorrisonRank1Update
|
||||
(
|
||||
const PtrList<scalarField>& r1Updates,
|
||||
const scalarField& p,
|
||||
const tmp<scalarField>& diag,
|
||||
const scalarField& mult,
|
||||
label n
|
||||
)
|
||||
{
|
||||
auto tAp(tmp<scalarField>::New(activeDesignVars_.size(), Zero));
|
||||
scalarField& Ap = tAp.ref();
|
||||
|
||||
if (n == 0)
|
||||
{
|
||||
Ap = invHessianVectorProduct(p, counter_, diag);
|
||||
return tAp;
|
||||
}
|
||||
|
||||
do
|
||||
{
|
||||
--n;
|
||||
Ap = ShermanMorrisonRank1Update(r1Updates, p, diag, mult, n);
|
||||
tmp<scalarField> tAv =
|
||||
ShermanMorrisonRank1Update(r1Updates, r1Updates[n], diag, mult, n);
|
||||
scalarField& Av = tAv.ref();
|
||||
scalar yHs = globalSum(r1Updates[n]*Av)/mult[n];
|
||||
Ap -= Av*globalSum(r1Updates[n]*Ap)/(1 + yHs)/mult[n];
|
||||
} while (n > 0);
|
||||
return tAp;
|
||||
}
|
||||
|
||||
|
||||
void Foam::ISQP::computeNewtonDirection()
|
||||
{
|
||||
// Zero the updates computed in the previous optimisation cycle
|
||||
@ -410,6 +497,7 @@ void Foam::ISQP::computeNewtonDirection()
|
||||
|
||||
// Solve equation for deltaP_. The expensive part of the step. Everything
|
||||
// else can be computed based on this
|
||||
addProfiling(ISQP, "ISQP::computeNewtonDirection");
|
||||
solveDeltaPEqn();
|
||||
|
||||
// deltaLamda
|
||||
@ -698,7 +786,10 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::invHFL()
|
||||
|
||||
Foam::tmp<Foam::scalarField> Foam::ISQP::resFGs()
|
||||
{
|
||||
tmp<scalarField> tFGs(tmp<scalarField>::New(gs_ + cValues_));
|
||||
tmp<scalarField> tFGs
|
||||
(
|
||||
tmp<scalarField>::New(gs_ + cValues_ - max((1 - cRed_)*cValues_, Zero))
|
||||
);
|
||||
scalarField& FGs = tFGs.ref();
|
||||
|
||||
forAll(constraintDerivatives_, cI)
|
||||
@ -795,6 +886,7 @@ Foam::tmp<Foam::scalarField> Foam::ISQP::resFz()
|
||||
|
||||
void Foam::ISQP::solveSubproblem()
|
||||
{
|
||||
addProfiling(ISQP, "ISQP::solveSubproblem");
|
||||
zeroUpdates();
|
||||
if (includeBoundConstraints_ || !cValues_.empty())
|
||||
{
|
||||
@ -825,13 +917,14 @@ void Foam::ISQP::solveSubproblem()
|
||||
(
|
||||
iter++ < maxNewtonIters_ && (eps_ > epsMin_ || resMax > 0.9*eps_)
|
||||
);
|
||||
Info<< "Finished solving the QP subproblem" << nl << endl;
|
||||
if (iter == maxNewtonIters_)
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Iterative solution of the QP problem did not converge"
|
||||
<< endl;
|
||||
}
|
||||
if (debug)
|
||||
if (debug > 2)
|
||||
{
|
||||
scalarField vars(designVars_().getVars(), activeDesignVars_);
|
||||
scalarField newVars(vars + p_);
|
||||
@ -961,8 +1054,21 @@ Foam::ISQP::ISQP
|
||||
coeffsDict(type).getOrDefault<label>("maxLineSearchIters", 10)
|
||||
),
|
||||
maxDxIters_(coeffsDict(type).getOrDefault<label>("maxDpIters", 1000)),
|
||||
relTol_(coeffsDict(type).getOrDefault<scalar>("relTol", 0.01)),
|
||||
preconType_
|
||||
(
|
||||
preconditionerNames.getOrDefault
|
||||
(
|
||||
"preconditioner", coeffsDict(type), preconditioners::diag
|
||||
)
|
||||
),
|
||||
cRed_
|
||||
(coeffsDict(type).getOrDefault<scalar>("targetConstraintReduction", 1)),
|
||||
meritFunctionFile_(nullptr)
|
||||
{
|
||||
Info<< "Preconditioner type of the SQP subproblem is ::"
|
||||
<< preconditionerNames.names()[preconType_]
|
||||
<< endl;
|
||||
// Always apply damping of s in ISQP
|
||||
useYDamping_ = true;
|
||||
useSDamping_ = false;
|
||||
|
||||
@ -34,14 +34,17 @@ Description
|
||||
solved using the interior point method (hence the I in the classe name,
|
||||
which is not standard in the terminology used in the literature). The
|
||||
(potentially dense) linear system formulated by the interior point method
|
||||
is solved using Conjugate Gradient with a diagonal preconditioner, using
|
||||
matrix-vector products to avoid storing the LHS matrix.
|
||||
is solved using Conjugate Gradient and a choice from three preconditioners
|
||||
(diagonal, inverse Hessian, Sherman-Morrison), using matrix-vector products
|
||||
to avoid storing the LHS matrix.
|
||||
|
||||
Bound constraints on the design variables will also be included, if set by
|
||||
the designVariables known by the updateMethod. If the QP problem is
|
||||
infeasible, the algorithm can still be used by setting includeExtraVars_
|
||||
to true, to allow a computation of an update, despite not being able to
|
||||
satisfy all the constraints of the QP problem.
|
||||
satisfy all the constraints of the QP problem. Alternatively,
|
||||
targetConstraintReduction can be set to a number lower than 1 to help with
|
||||
the feasibility of the dual problem.
|
||||
|
||||
|
||||
SourceFiles
|
||||
@ -69,6 +72,18 @@ class ISQP
|
||||
public LBFGS,
|
||||
public SQPBase
|
||||
{
|
||||
public:
|
||||
|
||||
//- Enumeration of preconditioners for the subproblem
|
||||
enum preconditioners
|
||||
{
|
||||
diag, invHessian, ShermanMorrison
|
||||
};
|
||||
|
||||
//- Names of preconditioners for the subproblem
|
||||
static const Enum<preconditioners> preconditionerNames;
|
||||
|
||||
|
||||
protected:
|
||||
|
||||
// Protected data
|
||||
@ -92,10 +107,6 @@ protected:
|
||||
|
||||
// Lagrange multipliers and slack variables
|
||||
|
||||
//- Lagrange multipliers for the inequality constraints
|
||||
// Inheritated from SQPBase
|
||||
//scalarField lamdas_;
|
||||
|
||||
//- Slack variables for the inequality constraints
|
||||
scalarField gs_;
|
||||
|
||||
@ -154,6 +165,16 @@ protected:
|
||||
//- Maxmimum number of iterations for the deltaX equation
|
||||
label maxDxIters_;
|
||||
|
||||
//- Relative tolerance of the CG solver, solving for deltaP_
|
||||
scalar relTol_;
|
||||
|
||||
//- Which preconditioner to use for the solution of the subproblem
|
||||
label preconType_;
|
||||
|
||||
//- Percentage reduction for the constraints.
|
||||
// Can be used to relax QP problems with no feasible points
|
||||
scalar cRed_;
|
||||
|
||||
//- File including the l1 merit function
|
||||
autoPtr<OFstream> meritFunctionFile_;
|
||||
|
||||
@ -193,20 +214,45 @@ protected:
|
||||
void solveDeltaPEqn();
|
||||
|
||||
//- Compute the RHS for the deltaX equation
|
||||
// Currently not in use
|
||||
tmp<scalarField> computeRHSForDeltaX(const scalarField& FDx);
|
||||
|
||||
//- CG algorithm for the solution of the deltaP eqn
|
||||
void CGforDeltaP(const scalarField& FDx);
|
||||
|
||||
//- Diagonal preconditioner of the deltaP eqn
|
||||
tmp<scalarField> deltaPDiagPreconditioner();
|
||||
|
||||
//- Procudt of the LHS of the deltaP eqn with a vector
|
||||
tmp<scalarField> DeltaPMatrixVectorProduct
|
||||
tmp<scalarField> matrixVectorProduct
|
||||
(
|
||||
const scalarField& vector
|
||||
);
|
||||
|
||||
//- Preconditioner-vector product for CG
|
||||
// In case a diagonal preconditioner is used, it is stored in
|
||||
// precon for all CG iterations
|
||||
tmp<scalarField> preconVectorProduct
|
||||
(
|
||||
const scalarField& vector,
|
||||
autoPtr<scalarField>& precon
|
||||
);
|
||||
|
||||
//- Diagonal preconditioner of the deltaP eqn
|
||||
tmp<scalarField> diagPreconditioner();
|
||||
|
||||
//- Use the Sherman-Morrison formula to compute the matrix-free
|
||||
//- product of the approximate inverse of the LHS with a vector
|
||||
tmp<scalarField> ShermanMorrisonPrecon(const scalarField& vector);
|
||||
|
||||
//- Recursive (naive) implementation of the rank-1 update
|
||||
// WIP, efficient for up to 2 flow-related constraints
|
||||
tmp<scalarField> ShermanMorrisonRank1Update
|
||||
(
|
||||
const PtrList<scalarField>& r1Updates,
|
||||
const scalarField& p,
|
||||
const tmp<scalarField>& diag,
|
||||
const scalarField& mult,
|
||||
label n
|
||||
);
|
||||
|
||||
//- Perform line search and return max residual corresponding to
|
||||
//- the updated solution
|
||||
scalar lineSearch();
|
||||
|
||||
@ -174,29 +174,36 @@ Foam::tmp<Foam::scalarField>
|
||||
Foam::LBFGS::invHessianVectorProduct
|
||||
(
|
||||
const scalarField& vector,
|
||||
const label counter
|
||||
const label counter,
|
||||
tmp<scalarField> diag
|
||||
)
|
||||
{
|
||||
// Sanity checks
|
||||
tmp<scalarField> tq(tmp<scalarField>::New(activeDesignVars_.size(), Zero));
|
||||
scalarField& q = tq.ref();
|
||||
if (vector.size() == designVars_().getVars().size())
|
||||
label nv = designVars_().getVars().size();
|
||||
label nav = activeDesignVars_.size();
|
||||
if (vector.size() == nv)
|
||||
{
|
||||
q.map(vector, activeDesignVars_);
|
||||
}
|
||||
else if (vector.size() == activeDesignVars_.size())
|
||||
else if (vector.size() == nav)
|
||||
{
|
||||
q = vector;
|
||||
}
|
||||
else
|
||||
{
|
||||
FatalErrorInFunction
|
||||
<< "Size of input vector is equal to neither the number of "
|
||||
<< " design variabes nor that of the active design variables"
|
||||
<< "Size of input vector "
|
||||
<< "(" << vector.size() << ") "
|
||||
<< "is equal to neither the number of design variabes "
|
||||
<< "(" << nv << ")"
|
||||
<< " nor that of the active design variables "
|
||||
<< "(" << nav << ")"
|
||||
<< exit(FatalError);
|
||||
}
|
||||
|
||||
if (counter != 0)
|
||||
if (counter)
|
||||
{
|
||||
// L-BFGS two loop recursion
|
||||
//~~~~~~~~~~~~~~~~~~~~~~~~~~
|
||||
@ -206,16 +213,21 @@ Foam::LBFGS::invHessianVectorProduct
|
||||
scalarField r(nSteps, 0.);
|
||||
for (label i = nLast; i > -1; --i)
|
||||
{
|
||||
//Info << "Y " << y_[i] << endl;
|
||||
//Info << "S " << s_[i] << endl;
|
||||
r[i] = 1./globalSum(y_[i]*s_[i]);
|
||||
a[i] = r[i]*globalSum(s_[i]*q);
|
||||
q -= a[i]*y_[i];
|
||||
}
|
||||
|
||||
scalar gamma =
|
||||
globalSum(y_[nLast]*s_[nLast])/globalSum(y_[nLast]*y_[nLast]);
|
||||
q *= gamma;
|
||||
globalSum(y_[nLast]*y_[nLast])/globalSum(y_[nLast]*s_[nLast]);
|
||||
if (diag)
|
||||
{
|
||||
q /= (gamma + diag());
|
||||
}
|
||||
else
|
||||
{
|
||||
q /= gamma;
|
||||
}
|
||||
|
||||
scalarField b(activeDesignVars_.size(), Zero);
|
||||
for (label i = 0; i < nSteps; ++i)
|
||||
@ -224,6 +236,10 @@ Foam::LBFGS::invHessianVectorProduct
|
||||
q += s_[i]*(a[i] - b);
|
||||
}
|
||||
}
|
||||
else if (diag)
|
||||
{
|
||||
q /= (1 + diag());
|
||||
}
|
||||
|
||||
return tq;
|
||||
}
|
||||
@ -243,6 +259,7 @@ Foam::LBFGS::HessianVectorProduct
|
||||
const label counter
|
||||
)
|
||||
{
|
||||
addProfiling(LBFGS, "LBFGS::HessianVectorProduct");
|
||||
// Sanity checks
|
||||
tmp<scalarField> tq(tmp<scalarField>::New(activeDesignVars_.size(), Zero));
|
||||
scalarField& q = tq.ref();
|
||||
@ -273,7 +290,7 @@ Foam::LBFGS::HessianVectorProduct
|
||||
|
||||
// Product of the last matrix on the right with the input vector
|
||||
scalarField SKsource(2*nSteps, Zero);
|
||||
for(label i = 0; i < nSteps; ++i)
|
||||
for (label i = 0; i < nSteps; ++i)
|
||||
{
|
||||
SKsource[i] = delta*globalSum(s_[i]*source);
|
||||
SKsource[i + nSteps] = globalSum(y_[i]*source);
|
||||
@ -374,11 +391,11 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::HessianDiag()
|
||||
|
||||
// Product of the inverse of the middle matrix with the right vector
|
||||
List<scalarField> MR(2*nSteps, scalarField(n, Zero));
|
||||
for(label k = 0; k < n; ++k)
|
||||
for (label k = 0; k < n; ++k)
|
||||
{
|
||||
for(label i = 0; i < 2*nSteps; ++i)
|
||||
for (label i = 0; i < 2*nSteps; ++i)
|
||||
{
|
||||
for(label j = 0; j < nSteps; ++j)
|
||||
for (label j = 0; j < nSteps; ++j)
|
||||
{
|
||||
MR[i][k] +=
|
||||
invM[i][j]*delta*s_[j][k]
|
||||
@ -390,9 +407,9 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::HessianDiag()
|
||||
// Part of the Hessian diagonal computed by the multiplication
|
||||
// of the above matrix with the left matrix of the recursive Hessian
|
||||
// reconstruction
|
||||
for(label k = 0; k < n; ++k)
|
||||
for (label k = 0; k < n; ++k)
|
||||
{
|
||||
for(label j = 0; j < nSteps; ++j)
|
||||
for (label j = 0; j < nSteps; ++j)
|
||||
{
|
||||
diag[k] -=
|
||||
delta*s_[j][k]*MR[j][k] + y_[j][k]*MR[j + nSteps][k];
|
||||
@ -448,7 +465,7 @@ Foam::LBFGS::SR1HessianVectorProduct
|
||||
|
||||
// Product of the last matrix on the right with the input vector
|
||||
scalarField YBSsource(nSteps, Zero);
|
||||
for(label i = 0; i < nSteps; ++i)
|
||||
for (label i = 0; i < nSteps; ++i)
|
||||
{
|
||||
YBSsource[i] = globalSum((y_[i] - delta*s_[i])*source);
|
||||
}
|
||||
@ -544,11 +561,11 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::SR1HessianDiag()
|
||||
|
||||
// Product of the inverse of the middle matrix with the right vector
|
||||
List<scalarField> MR(nSteps, scalarField(n, Zero));
|
||||
for(label k = 0; k < n; ++k)
|
||||
for (label k = 0; k < n; ++k)
|
||||
{
|
||||
for(label i = 0; i < nSteps; ++i)
|
||||
for (label i = 0; i < nSteps; ++i)
|
||||
{
|
||||
for(label j = 0; j < nSteps; ++j)
|
||||
for (label j = 0; j < nSteps; ++j)
|
||||
{
|
||||
MR[i][k] += invM[i][j]*(y_[j][k] - delta*s_[j][k]);
|
||||
}
|
||||
@ -558,9 +575,9 @@ Foam::tmp<Foam::scalarField> Foam::LBFGS::SR1HessianDiag()
|
||||
// Part of the Hessian diagonal computed by the multiplication
|
||||
// of the above matrix with the left matrix of the recursive Hessian
|
||||
// reconstruction
|
||||
for(label k = 0; k < n; ++k)
|
||||
for (label k = 0; k < n; ++k)
|
||||
{
|
||||
for(label j = 0; j < nSteps; ++j)
|
||||
for (label j = 0; j < nSteps; ++j)
|
||||
{
|
||||
diag[k] += (y_[j][k] - delta*s_[j][k])*MR[j][k];
|
||||
}
|
||||
|
||||
@ -117,9 +117,10 @@ protected:
|
||||
//- Update design variables
|
||||
virtual void update();
|
||||
|
||||
//- Compute the inverse Hessian - vector product
|
||||
//- Compute the inverse Hessian - vector product.
|
||||
// Input should have the size of all design variables or the active
|
||||
// ones, output is the size of the active design variables
|
||||
// ones, output is the size of the active design variables.
|
||||
// Can optionally received a seed diagonal matrix
|
||||
virtual tmp<scalarField> invHessianVectorProduct
|
||||
(
|
||||
const scalarField& vector
|
||||
@ -129,7 +130,8 @@ protected:
|
||||
tmp<scalarField> invHessianVectorProduct
|
||||
(
|
||||
const scalarField& vector,
|
||||
const label counter
|
||||
const label counter,
|
||||
tmp<scalarField> diag = nullptr
|
||||
);
|
||||
|
||||
//- Compute the Hessian - vector product
|
||||
|
||||
@ -174,16 +174,11 @@ Foam::updateMethod::inv(SquareMatrix<scalar> A)
|
||||
|
||||
Foam::scalar Foam::updateMethod::globalSum(const scalarField& field)
|
||||
{
|
||||
scalar value(0);
|
||||
if (globalSum_)
|
||||
{
|
||||
value = gSum(field);
|
||||
return gSum(field);
|
||||
}
|
||||
else
|
||||
{
|
||||
value = sum(field);
|
||||
}
|
||||
return value;
|
||||
return sum(field);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -1479,7 +1479,7 @@ scalar NURBS3DSurface::lengthU
|
||||
}
|
||||
|
||||
// Integrate.
|
||||
for(label uI = 0; uI<(uLenSize - 1); uI++)
|
||||
for (label uI = 0; uI<(uLenSize - 1); uI++)
|
||||
{
|
||||
const label ptI((uIStart+uI)*nVPts_ + vIConst);
|
||||
|
||||
@ -1512,7 +1512,7 @@ scalar NURBS3DSurface::lengthU
|
||||
}
|
||||
|
||||
// Integrate.
|
||||
for(label uI = 0; uI<(nPts - 1); uI++)
|
||||
for (label uI = 0; uI<(nPts - 1); uI++)
|
||||
{
|
||||
uLength +=
|
||||
0.5*(mag(dxdu[uI + 1]) + mag(dxdu[uI]))*(localU[uI + 1]-localU[uI]);
|
||||
@ -1556,7 +1556,7 @@ scalar NURBS3DSurface::lengthV
|
||||
}
|
||||
|
||||
// Integrate.
|
||||
for(label vI = 0; vI<(vLenSize - 1); vI++)
|
||||
for (label vI = 0; vI<(vLenSize - 1); vI++)
|
||||
{
|
||||
const label ptI((uIConst)*nVPts_ + (vIStart + vI));
|
||||
|
||||
@ -1589,7 +1589,7 @@ scalar NURBS3DSurface::lengthV
|
||||
}
|
||||
|
||||
// Integrate.
|
||||
for(label vI = 0; vI<(nPts - 1); vI++)
|
||||
for (label vI = 0; vI<(nPts - 1); vI++)
|
||||
{
|
||||
vLength +=
|
||||
0.5*(mag(dxdv[vI + 1]) + mag(dxdv[vI]))*(localV[vI + 1]-localV[vI]);
|
||||
@ -1979,7 +1979,7 @@ scalar NURBS3DSurface::lengthDerivativeU
|
||||
}
|
||||
|
||||
// Integrate.
|
||||
for(label uI=0; uI<(nPts-1); uI++)
|
||||
for (label uI=0; uI<(nPts-1); uI++)
|
||||
{
|
||||
ulDerivative +=
|
||||
0.5
|
||||
@ -2017,7 +2017,7 @@ scalar NURBS3DSurface::lengthDerivativeV
|
||||
}
|
||||
|
||||
// Integrate.
|
||||
for(label vI=0; vI<(nPts-1); vI++)
|
||||
for (label vI=0; vI<(nPts-1); vI++)
|
||||
{
|
||||
vlDerivative +=
|
||||
0.5
|
||||
@ -2046,9 +2046,9 @@ const labelList& NURBS3DSurface::getBoundaryCPIDs()
|
||||
|
||||
// v-constant cps
|
||||
label bID(0);
|
||||
for(label vI=0; vI<vNCPs; vI+=vNCPs-1)
|
||||
for (label vI=0; vI<vNCPs; vI+=vNCPs-1)
|
||||
{
|
||||
for(label uI=0; uI<uNCPs; uI++)
|
||||
for (label uI=0; uI<uNCPs; uI++)
|
||||
{
|
||||
const label CPI(vI*uNCPs + uI);
|
||||
whichBoundaryCPID_()[CPI] = bID;
|
||||
@ -2056,10 +2056,10 @@ const labelList& NURBS3DSurface::getBoundaryCPIDs()
|
||||
}
|
||||
}
|
||||
// u-constant cps
|
||||
for(label uI=0; uI<uNCPs; uI+=uNCPs-1)
|
||||
for (label uI=0; uI<uNCPs; uI+=uNCPs-1)
|
||||
{
|
||||
// corner CPS already accounted for
|
||||
for(label vI=1; vI<vNCPs-1; vI++)
|
||||
for (label vI=1; vI<vNCPs-1; vI++)
|
||||
{
|
||||
const label CPI(vI*uNCPs + uI);
|
||||
whichBoundaryCPID_()[CPI] = bID;
|
||||
|
||||
Reference in New Issue
Block a user