Skip to content

Commit d58cd15

Browse files
friedenhefriedenhe
authored
Enabled using objFunc std as convergence criterion. (#624)
* Enabled using objFunc std as convergence criterion. * Enables using objFuncStd as convergence criterion.
1 parent 2e1b6a2 commit d58cd15

File tree

5 files changed

+173
-22
lines changed

5 files changed

+173
-22
lines changed

‎dafoam/pyDAFoam.py

+4
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,10 @@ def __init__(self):
7474
## of magnitude (default) higher than this tolerance, the primal solution will return fail=True
7575
self.primalMinResTol = 1.0e-8
7676

77+
## The convergence tolerance based on the selected objective function's standard deviation
78+
## The standard deviation is calculated based on the last N (default 200) steps of the objective function series
79+
self.primalObjStdTol = {"active": False, "objFuncName": "None", "steps": 200, "tol": 1e-5, "tolDiff": 1e2}
80+
7781
## The boundary condition for primal solution. The keys should include "variable", "patch",
7882
## and "value". For turbulence variable, one can also set "useWallFunction" [bool].
7983
## Note that setting "primalBC" will overwrite any values defined in the "0" folder.

‎src/adjoint/DASolver/DASolver.C

+133-6
Original file line numberDiff line numberDiff line change
@@ -64,6 +64,9 @@ DASolver::DASolver(
6464
primalMinResTol_ = daOptionPtr_->getOption<scalar>("primalMinResTol");
6565
primalMinIters_ = daOptionPtr_->getOption<label>("primalMinIters");
6666

67+
// initialize the objStd variables.
68+
this->initObjStd();
69+
6770
Info << "DAOpton initialized " << endl;
6871
}
6972

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

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

138-
// check exit condition
139-
if (DAUtility::primalMaxInitRes_ < primalMinResTol_ && runTime.timeIndex() > primalMinIters_)
141+
// calculate the objective function standard deviation. It will be used in determining if the primal converges
142+
this->calcObjStd(runTime);
143+
144+
// check exit condition, we need to satisfy both the residual and objFunc std condition
145+
if (DAUtility::primalMaxInitRes_ < primalMinResTol_ && runTime.timeIndex() > primalMinIters_ && primalObjStd_ < primalObjStdTol_)
140146
{
141147
Info << "Time = " << t << endl;
148+
142149
Info << "Minimal residual " << DAUtility::primalMaxInitRes_ << " satisfied the prescribed tolerance " << primalMinResTol_ << endl
143150
<< endl;
151+
152+
if (primalObjStdActive_)
153+
{
154+
Info << "ObjFunc standard deviation " << primalObjStd_ << " satisfied the prescribed tolerance " << primalObjStdTol_ << endl
155+
<< endl;
156+
}
157+
144158
this->printAllObjFuncs();
145159
runTime.writeNow();
146160
prevPrimalSolTime_ = t;
@@ -160,6 +174,99 @@ label DASolver::loop(Time& runTime)
160174
}
161175
}
162176

177+
void DASolver::initObjStd()
178+
{
179+
/*
180+
Description:
181+
Initialize the objStd variables.
182+
*/
183+
184+
// check if the objective function std is used in determining the primal convergence
185+
primalObjStdActive_ = daOptionPtr_->getSubDictOption<label>("primalObjStdTol", "active");
186+
if (primalObjStdActive_)
187+
{
188+
// if it is active, read in the tolerance and set a large value for the initial std
189+
primalObjStdTol_ = daOptionPtr_->getSubDictOption<scalar>("primalObjStdTol", "tol");
190+
primalObjStd_ = 999.0;
191+
192+
label steps = daOptionPtr_->getSubDictOption<label>("primalObjStdTol", "steps");
193+
primalObjSeries_.setSize(steps, 0.0);
194+
195+
word objFuncNameWanted = daOptionPtr_->getSubDictOption<word>("primalObjStdTol", "objFuncName");
196+
197+
label objFuncNameFound = 0;
198+
const dictionary& objFuncDict = daOptionPtr_->getAllOptions().subDict("objFunc");
199+
forAll(objFuncDict.toc(), idxI)
200+
{
201+
word objFuncName = objFuncDict.toc()[idxI];
202+
if (objFuncName == objFuncNameWanted)
203+
{
204+
objFuncNameFound = 1;
205+
}
206+
}
207+
if (objFuncNameFound == 0)
208+
{
209+
FatalErrorIn("initObjStd") << "objStd->objFuncName not found! "
210+
<< abort(FatalError);
211+
}
212+
}
213+
else
214+
{
215+
// if it is not active, set primalObjStdTol_ > primalObjStd_, such that it will
216+
// always pass the condition in DASolver::loop (ignore primalObjStd)
217+
primalObjStdTol_ = 1e-5;
218+
primalObjStd_ = 0.0;
219+
}
220+
}
221+
222+
void DASolver::calcObjStd(Time& runTime)
223+
{
224+
/*
225+
Description:
226+
calculate the objective function's std, this will be used to stop the primal simulation and also
227+
evaluate whether the primal converges. We will start calculating the objStd when primalObjSeries_
228+
is filled at least once, i.e., runTime.timeIndex() >= steps
229+
*/
230+
231+
if (!primalObjStdActive_ || runTime.timeIndex() < 1)
232+
{
233+
return;
234+
}
235+
236+
label steps = daOptionPtr_->getSubDictOption<label>("primalObjStdTol", "steps");
237+
238+
word objFuncNameWanted = daOptionPtr_->getSubDictOption<word>("primalObjStdTol", "objFuncName");
239+
240+
scalar objFunPartSum = 0.0;
241+
forAll(daObjFuncPtrList_, idxI)
242+
{
243+
DAObjFunc& daObjFunc = daObjFuncPtrList_[idxI];
244+
word objFuncName = daObjFunc.getObjFuncName();
245+
if (objFuncName == objFuncNameWanted)
246+
{
247+
objFunPartSum += daObjFunc.getObjFuncValue();
248+
}
249+
}
250+
label seriesI = (runTime.timeIndex() - 1) % steps;
251+
primalObjSeries_[seriesI] = objFunPartSum;
252+
253+
if (runTime.timeIndex() >= steps)
254+
{
255+
scalar mean = 0;
256+
forAll(primalObjSeries_, idxI)
257+
{
258+
mean += primalObjSeries_[idxI];
259+
}
260+
mean /= steps;
261+
primalObjStd_ = 0.0;
262+
forAll(primalObjSeries_, idxI)
263+
{
264+
primalObjStd_ += (primalObjSeries_[idxI] - mean) * (primalObjSeries_[idxI] - mean);
265+
}
266+
primalObjStd_ = sqrt(primalObjStd_ / steps);
267+
}
268+
}
269+
163270
void DASolver::calcUnsteadyObjFuncs()
164271
{
165272
/*
@@ -248,6 +355,14 @@ void DASolver::printAllObjFuncs()
248355
<< "-" << objFuncPart
249356
<< "-" << daObjFunc.getObjFuncType()
250357
<< ": " << objFuncVal;
358+
if (primalObjStdActive_)
359+
{
360+
word objFuncNameWanted = daOptionPtr_->getSubDictOption<word>("primalObjStdTol", "objFuncName");
361+
if (objFuncNameWanted == objFuncName)
362+
{
363+
Info << " Std " << primalObjStd_;
364+
}
365+
}
251366
if (timeOperator == "average" || timeOperator == "sum")
252367
{
253368
Info << " Unsteady " << timeOperator << " " << unsteadyObjFuncs_[uKey];
@@ -8481,14 +8596,26 @@ label DASolver::checkResidualTol()
84818596
If yes, return 0 else return 1
84828597
*/
84838598

8484-
scalar tol = daOptionPtr_->getOption<scalar>("primalMinResTol");
8599+
// when checking the tolerance, we relax the criteria by tolMax
8600+
84858601
scalar tolMax = daOptionPtr_->getOption<scalar>("primalMinResTolDiff");
8486-
if (DAUtility::primalMaxInitRes_ / tol > tolMax)
8602+
scalar stdTolMax = daOptionPtr_->getSubDictOption<scalar>("primalObjStdTol", "tolDiff");
8603+
if (DAUtility::primalMaxInitRes_ / primalMinResTol_ > tolMax)
84878604
{
84888605
Info << "********************************************" << endl;
84898606
Info << "Primal min residual " << DAUtility::primalMaxInitRes_ << endl
84908607
<< "did not satisfy the prescribed tolerance "
8491-
<< tol << endl;
8608+
<< primalMinResTol_ << endl;
8609+
Info << "Primal solution failed!" << endl;
8610+
Info << "********************************************" << endl;
8611+
return 1;
8612+
}
8613+
else if (primalObjStd_ / primalObjStdTol_ > stdTolMax)
8614+
{
8615+
Info << "********************************************" << endl;
8616+
Info << "Primal objFunc standard deviation " << primalObjStd_ << endl
8617+
<< "did not satisfy the prescribed tolerance "
8618+
<< primalObjStdTol_ << endl;
84928619
Info << "Primal solution failed!" << endl;
84938620
Info << "********************************************" << endl;
84948621
return 1;

‎src/adjoint/DASolver/DASolver.H

+20-1
Original file line numberDiff line numberDiff line change
@@ -224,6 +224,19 @@ protected:
224224
/// step-averaged surfaceScalar states
225225
PtrList<surfaceScalarField> meanSurfaceScalarStates_;
226226

227+
/// whether to use the primal's std as the convergence criterion
228+
label primalObjStdActive_;
229+
230+
/// the step series of the primal objective, used if primalObjStdTol->active is True
231+
scalarList primalObjSeries_;
232+
233+
/// the standard deviation of the primal objective, used if primalObjStdTol->active is True
234+
scalar primalObjStd_;
235+
236+
/// the tolerance of the primal objective std, used if primalObjStdTol->active is True
237+
scalar primalObjStdTol_;
238+
239+
227240
public:
228241
/// Runtime type information
229242
TypeName("DASolver");
@@ -1294,13 +1307,19 @@ public:
12941307
const word fieldType,
12951308
const double timeName);
12961309

1310+
/// get the latest time solution from the case folder.
12971311
scalar getLatestTime()
12981312
{
1299-
// get the latest time solution from the case folder.
13001313
instantList timeDirs = runTimePtr_->findTimes(runTimePtr_->path(), runTimePtr_->constant());
13011314
scalar latestTime = timeDirs.last().value();
13021315
return latestTime;
13031316
}
1317+
1318+
/// initialize the objStd vars
1319+
void initObjStd();
1320+
1321+
/// calculate the objective function's std
1322+
void calcObjStd(Time& runTime);
13041323
};
13051324

13061325
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

‎tests/refs/DAFoam_Test_DARhoSimpleFoamUBendRef.txt

+14-14
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,26 @@
11
Dictionary Key: NU
2-
@value 322.2353921318797 1e-08 1e-10
2+
@value 322.5176001120205 1e-08 1e-10
33
Dictionary Key: PL
4-
@value 2.310286713734058 1e-08 1e-10
4+
@value 2.31161488891712 1e-08 1e-10
55
Dictionary Key: fail
66
@value 0 1e-08 1e-10
77
Dictionary Key: NU
88
Dictionary Key: shapey
9-
@value -140.1709592966036 0.0001 1e-06
10-
@value 60.78016818466478 0.0001 1e-06
11-
@value -68.73500447578726 0.0001 1e-06
9+
@value 152.1541024728527 0.0001 1e-06
10+
@value 138.432194899197 0.0001 1e-06
11+
@value 29.84397412606532 0.0001 1e-06
1212
Dictionary Key: shapez
13-
@value -4.362814513888678 0.0001 1e-06
14-
@value 6.619514670089002 0.0001 1e-06
15-
@value 8.125727969252937 0.0001 1e-06
13+
@value -20.5162746556844 0.0001 1e-06
14+
@value -168.3587625937871 0.0001 1e-06
15+
@value -86.23596162719087 0.0001 1e-06
1616
Dictionary Key: PL
1717
Dictionary Key: shapey
18-
@value 3.545442933817327 0.0001 1e-06
19-
@value 0.6419043697747193 0.0001 1e-06
20-
@value -0.8440552457793308 0.0001 1e-06
18+
@value 2.677963956813089 0.0001 1e-06
19+
@value 7.13696478920283 0.0001 1e-06
20+
@value 4.288518523483048 0.0001 1e-06
2121
Dictionary Key: shapez
22-
@value -0.0006037858402438943 0.0001 1e-06
23-
@value 0.2348244172148388 0.0001 1e-06
24-
@value 0.1332789140632325 0.0001 1e-06
22+
@value -1.12856634205818 0.0001 1e-06
23+
@value -13.40132223266928 0.0001 1e-06
24+
@value -7.605506838786253 0.0001 1e-06
2525
Dictionary Key: fail
2626
@value 0 0.0001 1e-06

‎tests/runTests_DARhoSimpleFoamUBend.py

+2-1
Original file line numberDiff line numberDiff line change
@@ -40,8 +40,9 @@
4040
"useAD": {"mode": "reverse"},
4141
"useMeanStates": {"active": True, "start": 0.5},
4242
"designSurfaces": ["ubend"],
43-
"primalMinResTol": 1e-5,
43+
"primalMinResTol": 1e0,
4444
"primalMinResTolDiff": 1e5,
45+
"primalObjStdTol": {"active": True, "objFuncName": "NU", "steps": 50, "tol": 0.27, "tolDiff": 1e2},
4546
"primalBC": {
4647
"U0": {"variable": "U", "patches": ["inlet"], "value": [U0, 0.0, 0.0]},
4748
"p0": {"variable": "p", "patches": ["outlet"], "value": [p0]},

0 commit comments

Comments
 (0)