50 if (
s.size() != m.
n())
53 <<
"scalar derivative and HessianInv matrix do not have the " 63 res[i] +=
s[j]*m[j][i];
73 const SquareMatrix<scalar>& m,
77 if (
s.size() != m.n())
80 <<
"scalar derivative and HessianInv matrix do not have the " 90 res[i] += m[i][j]*
s[j];
104 if (a.size() !=
b.size())
107 <<
"operands of outerProduct do not have the same dimension" 111 SquareMatrix<scalar> res(a.size(),
Zero);
116 res[i][j] = a[i]*
b[j];
134 <<
"LU decomposed A " <<
A <<
endl;
137 for (label j = 0; j <
n; j++)
143 for (label i = 0; i <
n; i++)
187 scalar value = globalSum(tfield());
210 Foam::updateMethod::updateMethod
215 const label nConstraints,
234 designVars_(designVars),
235 nConstraints_(nConstraints),
236 activeDesignVars_(designVars().activeDesignVariables()),
237 objectiveDerivatives_(designVars().getVars().size(),
Zero),
238 constraintDerivatives_(0),
240 objectiveValueOld_(0),
242 correction_(readOrZeroField(
"correction", designVars().getVars().size())),
243 cumulativeCorrection_(0),
245 counter_(getOrDefault<label>(
"counter",
Zero)),
246 initialEtaSet_(false),
247 correctionFolder_(mesh_.time().globalPath()/
"optimisation"/
"correction"),
248 globalSum_(designVars_->globalSum())
253 if (
dict.readIfPresent(
"eta",
eta_))
281 const label nConstraints
286 Info<<
"updateMethod type : " << modelType <<
endl;
288 auto* ctorPtr = dictionaryConstructorTable(modelType);
297 *dictionaryConstructorTablePtr_
302 (ctorPtr(
mesh,
dict, designVars, nConstraints, modelType));
316 objectiveDerivatives_ = derivs;
325 constraintDerivatives_ = derivs;
331 objectiveValue_ = value;
337 objectiveValueOld_ = value;
349 return objectiveValue_;
355 return objectiveValueOld_;
385 globalSum_ = useGlobalSum;
391 nConstraints_ = nConstraints;
397 return nConstraints_;
412 if (cumulativeCorrection_.empty())
414 cumulativeCorrection_.setSize(correction_.size(),
Zero);
417 cumulativeCorrection_ += correction_;
419 fileName correctionFile
421 correctionFolder_/
"correction" + mesh_.time().timeName()
423 fileName cumulativeCorrectionFile
425 correctionFolder_/
"cumulativeCorrection" + mesh_.time().timeName()
428 OFstream corFile(correctionFile);
429 OFstream cumulCorFile(cumulativeCorrectionFile);
433 << cI <<
" " << correction_[cI] <<
endl;
435 << cI <<
" " << cumulativeCorrection_[cI] <<
endl;
443 return objectiveValue_;
449 return globalSum(objectiveDerivatives_*correction_);
455 return initialEtaSet_;
464 correction_ = oldCorrection;
void setConstraintDeriv(const PtrList< scalarField > &derivs)
Set constraints derivative.
virtual scalar computeMeritFunction()
Compute merit function. Could be different than the objective in the presence of constraints.
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &f1)
void setConstaintsNumber(const label nConstraints)
Set the number of constraints.
virtual scalar meritFunctionDirectionalDerivative()
Directional derivative of the merit function, in the direction of the correction. Could be different ...
errorManipArg< error, int > exit(error &err, const int errNo=1)
label nConstraints() const
Get the number of constraints.
error FatalError
Error stream (stdout output on all processes), with additional 'FOAM FATAL ERROR' header text and sta...
A list of keyword definitions, which are a keyword followed by a number of values (eg...
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void LUDecompose(scalarSquareMatrix &matrix, labelList &pivotIndices)
LU decompose the matrix with pivoting.
static autoPtr< updateMethod > New(const fvMesh &mesh, const dictionary &dict, autoPtr< designVariables > &designVars, const label nConstraints)
Return a reference to the selected turbulence model.
Abstract base class for optimisation methods.
volSymmTensorField invA(inv(I *UEqn.A()+drag->Dcu()))
Ostream & endl(Ostream &os)
Add newline and flush stream.
Ostream & writeEntry(const keyType &key, const T &value)
Write a keyword/value entry.
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a T. FatalIOError if not found, or if the number of tokens is incorrect.
scalar getObjectiveValueOld() const
Get old objective value.
const scalarField rightMult(const SquareMatrix< scalar > &, const scalarField &)
bool & initialEtaSet()
Return whether initial eta was set.
T returnReduce(const T &value, const BinaryOp &bop, const int tag=UPstream::msgType(), const label comm=UPstream::worldComm)
Perform reduction on a copy, using specified binary operation.
scalar getObjectiveValue() const
Get objective value.
const scalarField leftMult(const scalarField &, const SquareMatrix< scalar > &)
void setConstraintValues(const scalarField &values)
Set values of constraints.
scalar eta_
Step multiplying the correction.
#define forAll(list, i)
Loop across all elements in list.
localIOdictionary is derived from IOdictionary but excludes parallel master reading.
label n() const noexcept
The number of columns.
List< T > values(const HashTable< T, Key, Hash > &tbl, const bool doSort=false)
List of values from HashTable, optionally sorted.
const dictionary & optionalSubDict(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
Find and return a sub-dictionary, otherwise return this dictionary.
fileName::Type type(const fileName &name, const bool followLink=true)
Return the file type: DIRECTORY or FILE, normally following symbolic links.
virtual void updateOldCorrection(const scalarField &oldCorrection)
Update old correction. Useful for quasi-newton methods coupled with line search.
Type gSum(const FieldField< Field, Type > &f)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for expressions::valueTypeCode::INVALID.
void setStep(const scalar eta)
Set step for optimisation methods.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
void modifyStep(const scalar multiplier)
Multiply step.
virtual bool writeAuxiliaryData()
Write non-continuation data, usually under the optimisation folder.
void setObjectiveValueOld(const scalar value)
Set old objective value.
A class for handling words, derived from Foam::string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static tmp< T > New(Args &&... args)
Construct tmp with forwarding arguments.
errorManip< error > abort(error &err)
#define DebugInfo
Report an information message using Foam::Info.
SquareMatrix< scalar > inv(SquareMatrix< scalar > A)
An Ostream is an abstract base class for all output systems (streams, files, token lists...
void setObjectiveValue(const scalar value)
Set objective value.
void setObjectiveDeriv(const scalarField &derivs)
Set objective derivative.
scalar globalSum(const scalarField &field)
Compute either global or local sum, based on globalSum flag.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
bool readIfPresent(const word &keyword, T &val, enum keyType::option matchOpt=keyType::REGEX) const
Find an entry if present, and assign to T val. FatalIOError if it is found and the number of tokens i...
OBJstream os(runTime.globalPath()/outputName)
defineTypeNameAndDebug(combustionModel, 0)
dictionary coeffsDict(const word &type) const
Return optional dictionary with parameters specific to each method.
void setGlobalSum(const bool useGlobalSum)
Set globalSum variable.
bool initialEtaSet_
Is initially set?
A list of pointers to objects of type <T>, with allocation/deallocation management of the pointers...
const scalarField & getConstraintValues() const
Get values of constraints.
Mesh data needed to do the Finite Volume discretisation.
SquareMatrix< scalar > outerProd(const scalarField &, const scalarField &)
static bool master(const label communicator=worldComm)
True if process corresponds to the master rank in the communicator.
label getCycle() const
Get optimisation cycle.
messageStream Info
Information stream (stdout output on master, null elsewhere)
Reading is optional [identical to READ_IF_PRESENT].
Pointer management similar to std::unique_ptr, with some additional methods and type checking...
List< label > labelList
A List of labels.
A class for managing temporary objects.
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))
static const Foam::dimensionedScalar A("", Foam::dimPressure, 611.21)
scalarField & returnCorrection()
Return the correction of the design variables.
Defines the attributes of an object for which implicit objectRegistry management is supported...
#define FatalIOErrorInLookup(ios, lookupTag, lookupName, lookupTable)
Report an error message using Foam::FatalIOError.
tmp< scalarField > readOrZeroField(const word &name, const label size)
Helper function to either read a scalarField of certain size from a dictionary, or construct a zero f...
virtual bool writeData(Ostream &os) const
Write continuation data under uniform.
void LUBacksubstitute(const scalarSquareMatrix &luMmatrix, const labelList &pivotIndices, List< Type > &source)
LU back-substitution with given source, returning the solution in the source.
IOerror FatalIOError
Error stream (stdout output on all processes), with additional 'FOAM FATAL IO ERROR' header text and ...
static constexpr const zero Zero
Global zero (0)