From d58cd15176d3b5572c8cbda193fde2b52c853122 Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Wed, 10 Apr 2024 08:28:37 -0500
Subject: [PATCH] Enabled using objFunc std as convergence criterion. (#624)

* Enabled using objFunc std as convergence criterion.

* Enables using objFuncStd as convergence criterion.
---
 dafoam/pyDAFoam.py                            |   4 +
 src/adjoint/DASolver/DASolver.C               | 139 +++++++++++++++++-
 src/adjoint/DASolver/DASolver.H               |  21 ++-
 .../DAFoam_Test_DARhoSimpleFoamUBendRef.txt   |  28 ++--
 tests/runTests_DARhoSimpleFoamUBend.py        |   3 +-
 5 files changed, 173 insertions(+), 22 deletions(-)

diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
index aab119d1..8be65864 100755
--- a/dafoam/pyDAFoam.py
+++ b/dafoam/pyDAFoam.py
@@ -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.
diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index b4322bc8..6bff5954 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -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;
 }
 
@@ -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...
     */
 
@@ -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;
@@ -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()
 {
     /*
@@ -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];
@@ -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;
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index d85a2950..7ca2630c 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -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");
@@ -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);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/tests/refs/DAFoam_Test_DARhoSimpleFoamUBendRef.txt b/tests/refs/DAFoam_Test_DARhoSimpleFoamUBendRef.txt
index a7721973..38209269 100755
--- a/tests/refs/DAFoam_Test_DARhoSimpleFoamUBendRef.txt
+++ b/tests/refs/DAFoam_Test_DARhoSimpleFoamUBendRef.txt
@@ -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
\ No newline at end of file
diff --git a/tests/runTests_DARhoSimpleFoamUBend.py b/tests/runTests_DARhoSimpleFoamUBend.py
index 5e749e41..505f6539 100755
--- a/tests/runTests_DARhoSimpleFoamUBend.py
+++ b/tests/runTests_DARhoSimpleFoamUBend.py
@@ -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]},