60 tmp<scalarField> designVarsConstraints =
designVars_().constraintValues();
61 if (designVarsConstraints)
72 for (adjointSolverManager& adjSolvManager : adjointSolvManagers_)
74 for (adjointSolver& solver : adjSolvManager.adjointSolvers())
76 if (!isA<adjointNull>(solver))
90 <<
setw(width) <<
"#Cycle" <<
" " 91 <<
setw(width) <<
"LineSearchIters" <<
" " 92 <<
setw(width) <<
"CycleCPUcost" <<
" " 93 <<
setw(width) <<
"CyclePrimalSolutions" <<
" " 94 <<
setw(width) <<
"CycleAdjointSolutions" <<
" " 95 <<
setw(width) <<
"TotalCPUcost" <<
" " 96 <<
setw(width) <<
"TotalPrimalSolutions" <<
" " 97 <<
setw(width) <<
"TotalAdjointSolutions" <<
endl;
104 label cyclePrimalSolutions(nPrimalsPerCycle_);
105 label cycleAdjointSolutions(nAdjointsPerCycle_);
106 label lineSearchIters(1);
109 lineSearchIters = lineSearch_().innerIter();
110 cyclePrimalSolutions *= lineSearchIters;
111 if (lineSearch_().computeGradient())
113 cycleAdjointSolutions *= lineSearchIters;
116 if (zeroAdjointSolns)
118 cycleAdjointSolutions = 0;
120 nPrimalSolutions_ += cyclePrimalSolutions;
121 nAdjointSolutions_ += cycleAdjointSolutions;
122 const scalar elapsedCpuTime = mesh_.time().elapsedCpuTime();
123 const scalar cycleCost = elapsedCpuTime - CPUcost_;
124 CPUcost_ = elapsedCpuTime;
127 <<
setw(width) << mesh_.time().timeName() <<
" " 128 <<
setw(width) << lineSearchIters <<
" " 129 <<
setw(width) << cycleCost <<
" " 130 <<
setw(width) << cyclePrimalSolutions <<
" " 131 <<
setw(width) << cycleAdjointSolutions <<
" " 132 <<
setw(width) << CPUcost_ <<
" " 133 <<
setw(width) << nPrimalSolutions_ <<
" " 134 <<
setw(width) << nAdjointSolutions_ <<
endl;
144 bool converged(
false);
146 if (designVarsThreshold_)
149 designVars_->activeDesignVariables();
151 const scalarField activeVars(designVars_->getVars(), activeVarIDs);
152 const scalar scaledCorrection =
156 <<
"Active vars " << activeVars <<
endl;
157 Info<<
"Max. scaled correction of the design variables = " 160 if (scaledCorrection < designVarsThreshold_())
162 Info<<
tab <<
"Design variables have converged " <<
endl;
167 if (objectiveThreshold_)
169 const scalar newObjective = computeObjectiveValue();
170 const scalar oldObjective = updateMethod_->getObjectiveValueOld();
171 const scalar relativeUpdate =
172 mag(newObjective - oldObjective)/(
mag(oldObjective) + SMALL);
173 Info<<
"Relative change of the objective value = " 176 if (relativeUpdate < objectiveThreshold_())
178 Info<<
tab <<
"Objective function has converged " <<
endl;
183 const scalarField& constraints = updateMethod_->getConstraintValues();
184 const scalar feasibility =
sum(
pos(constraints)*constraints);
185 Info<<
"Feasibility = " << feasibility <<
endl;
186 if (converged && feasibility < feasibilityThreshold_)
188 Info<<
"Stopping criteria met and all constraints satisfied." <<
nl 189 <<
"Optimisation has converged, stopping ..." <<
nl <<
nl 194 for (adjointSolverManager& am : adjointSolvManagers_)
196 for (adjointSolver& as : am.adjointSolvers())
199 as.getObjectiveManager().writeObjectives();
202 writeToCostFile(
true);
217 Foam::designVariablesUpdate::designVariablesUpdate
220 const dictionary&
dict,
222 autoPtr<designVariables>& designVars
228 designVars_(designVars),
234 dict_.subDict(
"updateMethod"),
243 dict_.subDict(
"updateMethod").subOrEmptyDict(
"lineSearch"),
248 CPUcostFile_(mesh_.time().globalPath()/
"optimisation"/
"CPUcost"),
250 nAdjointsPerCycle_(nAdjointSolvers()),
251 nPrimalSolutions_(nPrimalsPerCycle_),
252 nAdjointSolutions_(nAdjointsPerCycle_),
254 designVarsThreshold_(nullptr),
255 objectiveThreshold_(nullptr),
256 convergenceCriteriaDefined_(false),
257 feasibilityThreshold_
259 dict.subOrEmptyDict(
"convergence").
260 getOrDefault<scalar>(
"feasibilityThreshold", 1.
e-06)
264 if (convergenceDict.found(
"designVariables"))
267 (
new scalar(convergenceDict.get<scalar>(
"designVariables")));
269 if (convergenceDict.found(
"objective"))
272 (
new scalar(convergenceDict.get<scalar>(
"objective")));
279 <<
"Neither eta (updateMethod) or maxInitChange (designVariables) " 293 const auto& cnstrTable =
294 *(constrainedOptimisationMethod::dictionaryConstructorTablePtr_);
298 <<
"Found " << nConstr <<
" adjoint solvers corresponding to " 299 <<
"constraints but the optimisation method (" 301 <<
") is not a constrainedOptimisationMethod." <<
nl 302 <<
"Available constrainedOptimisationMethods:" <<
nl 303 << cnstrTable.sortedToc()
314 <<
"Did not find any adjoint solvers corresponding to " 315 <<
"constraints but the optimisation method (" 317 <<
") is a constrainedOptimisationMethod." <<
nl <<
nl 318 <<
"This can cause some constraintOptimisationMethods to misbehave." 320 <<
"Either the isConstraint bool is not set in one of the adjoint " 321 <<
"solvers or you should consider using an updateMethod " 322 <<
"that is not a constrainedOptimisationMethod" 328 Info<<
"Optimisation will run until the max. scaled correction " 334 Info<<
"Optimisation will run until the relative update of the " 340 Info<<
"No convergence criterion defined for optimsation" <<
nl 342 <<
" optimisation cycles " <<
nl <<
endl;
360 setOldObjectiveValue();
387 updateGradientsAndValues();
388 updateMethod_->computeCorrection();
392 if (!updateMethod_->initialEtaSet() || designVars_->resetEta())
394 const scalar eta(designVars_->computeEta(
correction));
395 updateMethod_->modifyStep(eta);
396 updateMethod_->initialEtaSet() =
true;
407 PtrList<scalarField> constraintSens;
408 scalar objectiveValue(
Zero);
409 DynamicList<scalar> constraintValues;
411 for (adjointSolverManager& adjSolvManager : adjointSolvManagers_)
413 const scalar opWeight = adjSolvManager.operatingPointWeight();
417 tmp<scalarField> tadjointSolverManagerSens =
418 adjSolvManager.aggregateSensitivities();
422 objectiveValue += opWeight*adjSolvManager.objectiveValue();
425 PtrList<scalarField> adjointSolverManagerConstSens =
426 adjSolvManager.constraintSensitivities();
431 designVars_->postProcessSens
433 tadjointSolverManagerSens.ref(),
434 adjointSolverManagerConstSens,
435 adjSolvManager.adjointSolversNames(),
436 adjSolvManager.isMaster()
439 if (objectiveSens.empty())
441 objectiveSens.setSize(tadjointSolverManagerSens().size(),
Zero);
445 objectiveSens += opWeight*tadjointSolverManagerSens();
446 forAll(adjointSolverManagerConstSens, sI)
449 push_back(adjointSolverManagerConstSens.set(sI,
nullptr));
451 constraintValues.push_back(adjSolvManager.constraintValues());
454 tmp<scalarField> designVarsConstValues = designVars_->constraintValues();
455 PtrList<scalarField> designVarsConstDerivs =
456 designVars_->constraintDerivatives();
457 if (designVarsConstValues && designVarsConstDerivs.size())
459 if (designVarsConstValues().size() != designVarsConstDerivs.size())
462 <<
"Size of design variables constraints and derivatives differ" 466 constraintValues.push_back(designVarsConstValues());
467 constraintSens.push_back(std::move(designVarsConstDerivs));
471 updateMethod_->setObjectiveDeriv(objectiveSens);
472 updateMethod_->setConstraintDeriv(constraintSens);
473 updateMethod_->setObjectiveValue(objectiveValue);
474 updateMethod_->setConstraintValues
481 scalar objectiveValue(
Zero);
484 const scalar opWeight = adjSolvManager.operatingPointWeight();
485 objectiveValue += opWeight*adjSolvManager.objectiveValue();
487 return objectiveValue;
493 updateMethod_->setObjectiveValueOld(computeObjectiveValue());
501 scalar objectiveValue(
Zero);
506 const scalar opWeight = adjSolvManager.operatingPointWeight();
508 objectiveValue += opWeight*adjSolvManager.objectiveValue();
509 constraintValues.
push_back(adjSolvManager.constraintValues());
513 tmp<scalarField> designVarsConstValues = designVars_->constraintValues();
514 if (designVarsConstValues)
516 constraintValues.
push_back(designVarsConstValues());
518 updateMethod_->setObjectiveValue(objectiveValue);
519 updateMethod_->setConstraintValues
522 return updateMethod_->computeMeritFunction();
528 return updateMethod_->meritFunctionDirectionalDerivative();
537 updateMethod_->updateOldCorrection(oldCorrection);
543 updateMethod_->writeAuxiliaryData();
544 designVars_->writeDesignVars();
551 updateOldCorrection(oldCorrection);
553 designVars_->evolveNumber();
556 lineSearch_->postUpdate();
566 as.getObjectiveManager().getObjectiveFunctions();
569 obj.writeInstantaneousSeparator();
574 if (convergenceCriteriaDefined_)
576 checkConvergence(oldCorrection);
PtrList< adjointSolverManager > & adjointSolvManagers_
tmp< fvMatrix< Type > > correction(const fvMatrix< Type > &)
Return the correction form of the given matrix by subtracting the matrix multiplied by the current fi...
A class encapsulating functionality neccessary to perform an optimisation loop, such as updating the ...
void writeCPUcostHeader()
Write CPU cost header.
const Type & value() const noexcept
Return const reference to value.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
void postUpdate(const scalarField &oldCorrection)
Steps to be executed after the susccessfull update of the design varibles, i.e. the last step of line...
errorManipArg< error, int > exit(error &err, const int errNo=1)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
label nAdjointSolvers() const
Get total number of adjoint solvers.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
constexpr char nl
The newline '\n' character (0x0a)
bool convergenceCriteriaDefined_
Is at least a single convergence criterion defined.
Ostream & endl(Ostream &os)
Add newline and flush stream.
static bool & parRun() noexcept
Test if this a parallel run.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tf1, const word &name, const dimensionSet &dimensions, const bool initCopy=false)
Global function forwards to reuseTmpDimensionedField::New.
static unsigned int defaultPrecision() noexcept
Return the default precision.
Base class for adjoint solvers.
constexpr char tab
The tab '\t' character(0x09)
virtual dimensionedScalar endTime() const
Return end time.
const Time & time() const
Return the top-level database.
autoPtr< designVariables > & designVars_
label nConstraints(PtrList< adjointSolverManager > &adjointSolverManagers) const
Get the number of adjoint solvers that correspond to constraints.
#define forAll(list, i)
Loop across all elements in list.
void reset(T *p=nullptr) noexcept
Delete managed object and set to new given pointer.
dimensionedScalar pos(const dimensionedScalar &ds)
autoPtr< scalar > objectiveThreshold_
The relative update of the objective value w.r.t. to its last value should be smaller than this value...
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
tmp< scalarField > computeDirection()
Compute update direction.
const dimensionedScalar e
Elementary charge.
static void exit(int errNo=1)
Shutdown (finalize) MPI as required and exit program with errNo.
scalar feasibilityThreshold_
In case of a constrained optimisation, the sum of positive constraints should be lower than this valu...
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void write()
Write useful quantities to files.
scalar computeObjectiveValue()
Sum objective values from all adjointSolverManagers.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
#define DebugInfo
Report an information message using Foam::Info.
Istream and Ostream manipulators taking arguments.
Abstract base class for objective functions. No point in making this runTime selectable since its chi...
Type gMax(const FieldField< Field, Type > &f)
defineTypeNameAndDebug(combustionModel, 0)
scalar meritFunctionDirectionalDerivative()
Derivative of the merit function.
scalar computeMeritFunction()
Compute the merit function of the optimisation problem.
void push_back(const T &val)
Copy append an element to the end of this list.
void writeToCostFile(bool zeroAdjointSolns=false)
Write to cost file.
Class for managing adjoint solvers, which may be more than one per operating point.
autoPtr< scalar > designVarsThreshold_
The maximum of the correction/designVariables values must be lower that this threshold to consider th...
#define WarningInFunction
Report a warning using Foam::Warning.
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
void updateOldCorrection(const scalarField &)
Update old correction. Needed for quasi-Newton Methods.
void update()
Update design variables.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Omanip< int > setw(const int i)
Macros to ease declaration of run-time selection tables.
autoPtr< updateMethod > updateMethod_
Method to update the design variables based on the provided sensitivity derivatives.
A class for managing temporary objects.
void updateGradientsAndValues()
Compute cumulative objective and constraint gradients.
PtrList< adjointSolverManager > & adjointSolverManagers
void setOldObjectiveValue()
Set the old objective value known by the updateMethod.
void checkConvergence(const scalarField &oldCorrection)
Check if the optimisation loop has converged based on the provided criteria.
static constexpr const zero Zero
Global zero (0)