Skip to content

Commit

Permalink
Enabled using objFunc std as convergence criterion. (mdolab#624)
Browse files Browse the repository at this point in the history
* Enabled using objFunc std as convergence criterion.

* Enables using objFuncStd as convergence criterion.
  • Loading branch information
friedenhe authored Apr 10, 2024
1 parent 2e1b6a2 commit d58cd15
Show file tree
Hide file tree
Showing 5 changed files with 173 additions and 22 deletions.
4 changes: 4 additions & 0 deletions dafoam/pyDAFoam.py
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,10 @@ def __init__(self):
## of magnitude (default) higher than this tolerance, the primal solution will return fail=True
self.primalMinResTol = 1.0e-8

## The convergence tolerance based on the selected objective function's standard deviation
## The standard deviation is calculated based on the last N (default 200) steps of the objective function series
self.primalObjStdTol = {"active": False, "objFuncName": "None", "steps": 200, "tol": 1e-5, "tolDiff": 1e2}

## The boundary condition for primal solution. The keys should include "variable", "patch",
## and "value". For turbulence variable, one can also set "useWallFunction" [bool].
## Note that setting "primalBC" will overwrite any values defined in the "0" folder.
Expand Down
139 changes: 133 additions & 6 deletions src/adjoint/DASolver/DASolver.C
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,9 @@ DASolver::DASolver(
primalMinResTol_ = daOptionPtr_->getOption<scalar>("primalMinResTol");
primalMinIters_ = daOptionPtr_->getOption<label>("primalMinIters");

// initialize the objStd variables.
this->initObjStd();

Info << "DAOpton initialized " << endl;
}

Expand Down Expand Up @@ -116,7 +119,7 @@ label DASolver::loop(Time& runTime)
{
/*
Description:
The loop method to increment the runtime. The reason we implent this is
The loop method to increment the runtime. The reason we implement this is
because the runTime.loop() and simple.loop() give us seg fault...
*/

Expand All @@ -135,12 +138,23 @@ label DASolver::loop(Time& runTime)
funcObj.execute();
}

// check exit condition
if (DAUtility::primalMaxInitRes_ < primalMinResTol_ && runTime.timeIndex() > primalMinIters_)
// calculate the objective function standard deviation. It will be used in determining if the primal converges
this->calcObjStd(runTime);

// check exit condition, we need to satisfy both the residual and objFunc std condition
if (DAUtility::primalMaxInitRes_ < primalMinResTol_ && runTime.timeIndex() > primalMinIters_ && primalObjStd_ < primalObjStdTol_)
{
Info << "Time = " << t << endl;

Info << "Minimal residual " << DAUtility::primalMaxInitRes_ << " satisfied the prescribed tolerance " << primalMinResTol_ << endl
<< endl;

if (primalObjStdActive_)
{
Info << "ObjFunc standard deviation " << primalObjStd_ << " satisfied the prescribed tolerance " << primalObjStdTol_ << endl
<< endl;
}

this->printAllObjFuncs();
runTime.writeNow();
prevPrimalSolTime_ = t;
Expand All @@ -160,6 +174,99 @@ label DASolver::loop(Time& runTime)
}
}

void DASolver::initObjStd()
{
/*
Description:
Initialize the objStd variables.
*/

// check if the objective function std is used in determining the primal convergence
primalObjStdActive_ = daOptionPtr_->getSubDictOption<label>("primalObjStdTol", "active");
if (primalObjStdActive_)
{
// if it is active, read in the tolerance and set a large value for the initial std
primalObjStdTol_ = daOptionPtr_->getSubDictOption<scalar>("primalObjStdTol", "tol");
primalObjStd_ = 999.0;

label steps = daOptionPtr_->getSubDictOption<label>("primalObjStdTol", "steps");
primalObjSeries_.setSize(steps, 0.0);

word objFuncNameWanted = daOptionPtr_->getSubDictOption<word>("primalObjStdTol", "objFuncName");

label objFuncNameFound = 0;
const dictionary& objFuncDict = daOptionPtr_->getAllOptions().subDict("objFunc");
forAll(objFuncDict.toc(), idxI)
{
word objFuncName = objFuncDict.toc()[idxI];
if (objFuncName == objFuncNameWanted)
{
objFuncNameFound = 1;
}
}
if (objFuncNameFound == 0)
{
FatalErrorIn("initObjStd") << "objStd->objFuncName not found! "
<< abort(FatalError);
}
}
else
{
// if it is not active, set primalObjStdTol_ > primalObjStd_, such that it will
// always pass the condition in DASolver::loop (ignore primalObjStd)
primalObjStdTol_ = 1e-5;
primalObjStd_ = 0.0;
}
}

void DASolver::calcObjStd(Time& runTime)
{
/*
Description:
calculate the objective function's std, this will be used to stop the primal simulation and also
evaluate whether the primal converges. We will start calculating the objStd when primalObjSeries_
is filled at least once, i.e., runTime.timeIndex() >= steps
*/

if (!primalObjStdActive_ || runTime.timeIndex() < 1)
{
return;
}

label steps = daOptionPtr_->getSubDictOption<label>("primalObjStdTol", "steps");

word objFuncNameWanted = daOptionPtr_->getSubDictOption<word>("primalObjStdTol", "objFuncName");

scalar objFunPartSum = 0.0;
forAll(daObjFuncPtrList_, idxI)
{
DAObjFunc& daObjFunc = daObjFuncPtrList_[idxI];
word objFuncName = daObjFunc.getObjFuncName();
if (objFuncName == objFuncNameWanted)
{
objFunPartSum += daObjFunc.getObjFuncValue();
}
}
label seriesI = (runTime.timeIndex() - 1) % steps;
primalObjSeries_[seriesI] = objFunPartSum;

if (runTime.timeIndex() >= steps)
{
scalar mean = 0;
forAll(primalObjSeries_, idxI)
{
mean += primalObjSeries_[idxI];
}
mean /= steps;
primalObjStd_ = 0.0;
forAll(primalObjSeries_, idxI)
{
primalObjStd_ += (primalObjSeries_[idxI] - mean) * (primalObjSeries_[idxI] - mean);
}
primalObjStd_ = sqrt(primalObjStd_ / steps);
}
}

void DASolver::calcUnsteadyObjFuncs()
{
/*
Expand Down Expand Up @@ -248,6 +355,14 @@ void DASolver::printAllObjFuncs()
<< "-" << objFuncPart
<< "-" << daObjFunc.getObjFuncType()
<< ": " << objFuncVal;
if (primalObjStdActive_)
{
word objFuncNameWanted = daOptionPtr_->getSubDictOption<word>("primalObjStdTol", "objFuncName");
if (objFuncNameWanted == objFuncName)
{
Info << " Std " << primalObjStd_;
}
}
if (timeOperator == "average" || timeOperator == "sum")
{
Info << " Unsteady " << timeOperator << " " << unsteadyObjFuncs_[uKey];
Expand Down Expand Up @@ -8481,14 +8596,26 @@ label DASolver::checkResidualTol()
If yes, return 0 else return 1
*/

scalar tol = daOptionPtr_->getOption<scalar>("primalMinResTol");
// when checking the tolerance, we relax the criteria by tolMax

scalar tolMax = daOptionPtr_->getOption<scalar>("primalMinResTolDiff");
if (DAUtility::primalMaxInitRes_ / tol > tolMax)
scalar stdTolMax = daOptionPtr_->getSubDictOption<scalar>("primalObjStdTol", "tolDiff");
if (DAUtility::primalMaxInitRes_ / primalMinResTol_ > tolMax)
{
Info << "********************************************" << endl;
Info << "Primal min residual " << DAUtility::primalMaxInitRes_ << endl
<< "did not satisfy the prescribed tolerance "
<< tol << endl;
<< primalMinResTol_ << endl;
Info << "Primal solution failed!" << endl;
Info << "********************************************" << endl;
return 1;
}
else if (primalObjStd_ / primalObjStdTol_ > stdTolMax)
{
Info << "********************************************" << endl;
Info << "Primal objFunc standard deviation " << primalObjStd_ << endl
<< "did not satisfy the prescribed tolerance "
<< primalObjStdTol_ << endl;
Info << "Primal solution failed!" << endl;
Info << "********************************************" << endl;
return 1;
Expand Down
21 changes: 20 additions & 1 deletion src/adjoint/DASolver/DASolver.H
Original file line number Diff line number Diff line change
Expand Up @@ -224,6 +224,19 @@ protected:
/// step-averaged surfaceScalar states
PtrList<surfaceScalarField> meanSurfaceScalarStates_;

/// whether to use the primal's std as the convergence criterion
label primalObjStdActive_;

/// the step series of the primal objective, used if primalObjStdTol->active is True
scalarList primalObjSeries_;

/// the standard deviation of the primal objective, used if primalObjStdTol->active is True
scalar primalObjStd_;

/// the tolerance of the primal objective std, used if primalObjStdTol->active is True
scalar primalObjStdTol_;


public:
/// Runtime type information
TypeName("DASolver");
Expand Down Expand Up @@ -1294,13 +1307,19 @@ public:
const word fieldType,
const double timeName);

/// get the latest time solution from the case folder.
scalar getLatestTime()
{
// get the latest time solution from the case folder.
instantList timeDirs = runTimePtr_->findTimes(runTimePtr_->path(), runTimePtr_->constant());
scalar latestTime = timeDirs.last().value();
return latestTime;
}

/// initialize the objStd vars
void initObjStd();

/// calculate the objective function's std
void calcObjStd(Time& runTime);
};

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Expand Down
28 changes: 14 additions & 14 deletions tests/refs/DAFoam_Test_DARhoSimpleFoamUBendRef.txt
Original file line number Diff line number Diff line change
@@ -1,26 +1,26 @@
Dictionary Key: NU
@value 322.2353921318797 1e-08 1e-10
@value 322.5176001120205 1e-08 1e-10
Dictionary Key: PL
@value 2.310286713734058 1e-08 1e-10
@value 2.31161488891712 1e-08 1e-10
Dictionary Key: fail
@value 0 1e-08 1e-10
Dictionary Key: NU
Dictionary Key: shapey
@value -140.1709592966036 0.0001 1e-06
@value 60.78016818466478 0.0001 1e-06
@value -68.73500447578726 0.0001 1e-06
@value 152.1541024728527 0.0001 1e-06
@value 138.432194899197 0.0001 1e-06
@value 29.84397412606532 0.0001 1e-06
Dictionary Key: shapez
@value -4.362814513888678 0.0001 1e-06
@value 6.619514670089002 0.0001 1e-06
@value 8.125727969252937 0.0001 1e-06
@value -20.5162746556844 0.0001 1e-06
@value -168.3587625937871 0.0001 1e-06
@value -86.23596162719087 0.0001 1e-06
Dictionary Key: PL
Dictionary Key: shapey
@value 3.545442933817327 0.0001 1e-06
@value 0.6419043697747193 0.0001 1e-06
@value -0.8440552457793308 0.0001 1e-06
@value 2.677963956813089 0.0001 1e-06
@value 7.13696478920283 0.0001 1e-06
@value 4.288518523483048 0.0001 1e-06
Dictionary Key: shapez
@value -0.0006037858402438943 0.0001 1e-06
@value 0.2348244172148388 0.0001 1e-06
@value 0.1332789140632325 0.0001 1e-06
@value -1.12856634205818 0.0001 1e-06
@value -13.40132223266928 0.0001 1e-06
@value -7.605506838786253 0.0001 1e-06
Dictionary Key: fail
@value 0 0.0001 1e-06
3 changes: 2 additions & 1 deletion tests/runTests_DARhoSimpleFoamUBend.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,8 +40,9 @@
"useAD": {"mode": "reverse"},
"useMeanStates": {"active": True, "start": 0.5},
"designSurfaces": ["ubend"],
"primalMinResTol": 1e-5,
"primalMinResTol": 1e0,
"primalMinResTolDiff": 1e5,
"primalObjStdTol": {"active": True, "objFuncName": "NU", "steps": 50, "tol": 0.27, "tolDiff": 1e2},
"primalBC": {
"U0": {"variable": "U", "patches": ["inlet"], "value": [U0, 0.0, 0.0]},
"p0": {"variable": "p", "patches": ["outlet"], "value": [p0]},
Expand Down

0 comments on commit d58cd15

Please sign in to comment.