Changed how the primalResTol is calculated (#587)

* Changed how primalResTol is computed. * Fixed tests.
mdolab · Feb 15, 2024 · 376843f · 376843f
1 parent 561f127
commit 376843f
Show file tree

Hide file tree

Showing 50 changed files with 178 additions and 203 deletions.
diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
@@ -704,8 +704,8 @@ def __init__(self):
         }
 
         ## number of minimal primal iterations. The primal has to run this many iterations, even the primal residual
-        ## has reduced below the tolerance. The default is a negative value (always satisfied).
-        self.primalMinIters = -1
+        ## has reduced below the tolerance. The default is 1: the primal has to run for at least one iteration
+        self.primalMinIters = 1
 
         ## tensorflow related functions
         self.tensorflow = {

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C
@@ -454,11 +454,8 @@ void DASpalartAllmaras::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverNuTilda = solve(nuTildaEqn);
-        if (printToScreen)
-        {
-            Info << "nuTilda Initial residual: " << solverNuTilda.initialResidual() << endl
-                 << "          Final residual: " << solverNuTilda.finalResidual() << endl;
-        }
+
+        DAUtility::primalResidualControl(solverNuTilda, printToScreen, "nuTilda");
 
         DAUtility::boundVar(allOptions_, nuTilda_, printToScreen);
         nuTilda_.correctBoundaryConditions();

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C
@@ -479,11 +479,7 @@ void DASpalartAllmarasFv3::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverNuTilda = solve(nuTildaEqn);
-        if (printToScreen)
-        {
-            Info << "nuTilda Initial residual: " << solverNuTilda.initialResidual() << endl
-                 << "          Final residual: " << solverNuTilda.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverNuTilda, printToScreen, "nuTilda");
 
         DAUtility::boundVar(allOptions_, nuTilda_, printToScreen);
         nuTilda_.correctBoundaryConditions();

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3FieldInversion.C b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3FieldInversion.C
@@ -465,11 +465,7 @@ void DASpalartAllmarasFv3FieldInversion::calcResiduals(const dictionary& options
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverNuTilda = solve(nuTildaEqn);
-        if (printToScreen)
-        {
-            Info << "nuTilda Initial residual: " << solverNuTilda.initialResidual() << endl
-                 << "          Final residual: " << solverNuTilda.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverNuTilda, printToScreen, "nuTilda");
 
         DAUtility::boundVar(allOptions_, nuTilda_, printToScreen);
         nuTilda_.correctBoundaryConditions();

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DAkEpsilon.C b/src/adjoint/DAModel/DATurbulenceModel/DAkEpsilon.C
@@ -572,11 +572,7 @@ void DAkEpsilon::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverEpsilon = solve(epsEqn);
-        if (printToScreen)
-        {
-            Info << "epsilon Initial residual: " << solverEpsilon.initialResidual() << endl
-                 << "          Final residual: " << solverEpsilon.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverEpsilon, printToScreen, "epsilon");
 
         DAUtility::boundVar(allOptions_, epsilon_, printToScreen);
     }
@@ -610,11 +606,7 @@ void DAkEpsilon::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverK = solve(kEqn);
-        if (printToScreen)
-        {
-            Info << "k Initial residual: " << solverK.initialResidual() << endl
-                 << "    Final residual: " << solverK.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverK, printToScreen, "k");
 
         DAUtility::boundVar(allOptions_, k_, printToScreen);
 

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DAkOmega.C b/src/adjoint/DAModel/DATurbulenceModel/DAkOmega.C
@@ -577,11 +577,7 @@ void DAkOmega::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverOmega = solve(omegaEqn);
-        if (printToScreen)
-        {
-            Info << "omega Initial residual: " << solverOmega.initialResidual() << endl
-                 << "        Final residual: " << solverOmega.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverOmega, printToScreen, "omega");
 
         DAUtility::boundVar(allOptions_, omega_, printToScreen);
     }
@@ -614,11 +610,7 @@ void DAkOmega::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverK = solve(kEqn);
-        if (printToScreen)
-        {
-            Info << "k Initial residual: " << solverK.initialResidual() << endl
-                 << "    Final residual: " << solverK.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverK, printToScreen, "k");
 
         DAUtility::boundVar(allOptions_, k_, printToScreen);
 

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaFieldInversionOmega.C b/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaFieldInversionOmega.C
@@ -581,11 +581,7 @@ void DAkOmegaFieldInversionOmega::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverOmega = solve(omegaEqn);
-        if (printToScreen)
-        {
-            Info << "omega Initial residual: " << solverOmega.initialResidual() << endl
-                 << "        Final residual: " << solverOmega.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverOmega, printToScreen, "omega");
 
         DAUtility::boundVar(allOptions_, omega_, printToScreen);
     }
@@ -618,11 +614,7 @@ void DAkOmegaFieldInversionOmega::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverK = solve(kEqn);
-        if (printToScreen)
-        {
-            Info << "k Initial residual: " << solverK.initialResidual() << endl
-                 << "    Final residual: " << solverK.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverK, printToScreen, "k");
 
         DAUtility::boundVar(allOptions_, k_, printToScreen);
 

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSST.C b/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSST.C
@@ -741,11 +741,7 @@ void DAkOmegaSST::calcResiduals(const dictionary& options)
             // get the solver performance info such as initial
             // and final residuals
             SolverPerformance<scalar> solverOmega = solve(omegaEqn);
-            if (printToScreen)
-            {
-                Info << "omega Initial residual: " << solverOmega.initialResidual() << endl
-                     << "        Final residual: " << solverOmega.finalResidual() << endl;
-            }
+            DAUtility::primalResidualControl(solverOmega, printToScreen, "omega");
 
             DAUtility::boundVar(allOptions_, omega_, printToScreen);
         }
@@ -782,11 +778,7 @@ void DAkOmegaSST::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverK = solve(kEqn);
-        if (printToScreen)
-        {
-            Info << "k Initial residual: " << solverK.initialResidual() << endl
-                 << "    Final residual: " << solverK.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverK, printToScreen, "k");
 
         DAUtility::boundVar(allOptions_, k_, printToScreen);
 

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSSTFIML.C b/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSSTFIML.C
@@ -1047,11 +1047,7 @@ void DAkOmegaSSTFIML::calcResiduals(const dictionary& options)
             // get the solver performance info such as initial
             // and final residuals
             SolverPerformance<scalar> solverOmega = solve(omegaEqn);
-            if (printToScreen)
-            {
-                Info << "omega Initial residual: " << solverOmega.initialResidual() << endl
-                     << "        Final residual: " << solverOmega.finalResidual() << endl;
-            }
+            DAUtility::primalResidualControl(solverOmega, printToScreen, "omega");
 
             DAUtility::boundVar(allOptions_, omega_, printToScreen);
         }
@@ -1088,11 +1084,7 @@ void DAkOmegaSSTFIML::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverK = solve(kEqn);
-        if (printToScreen)
-        {
-            Info << "k Initial residual: " << solverK.initialResidual() << endl
-                 << "    Final residual: " << solverK.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverK, printToScreen, "k");
 
         DAUtility::boundVar(allOptions_, k_, printToScreen);
 

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSSTFieldInversion.C b/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSSTFieldInversion.C
@@ -792,11 +792,7 @@ void DAkOmegaSSTFieldInversion::calcResiduals(const dictionary& options)
             // get the solver performance info such as initial
             // and final residuals
             SolverPerformance<scalar> solverOmega = solve(omegaEqn);
-            if (printToScreen)
-            {
-                Info << "omega Initial residual: " << solverOmega.initialResidual() << endl
-                     << "        Final residual: " << solverOmega.finalResidual() << endl;
-            }
+            DAUtility::primalResidualControl(solverOmega, printToScreen, "omega");
 
             DAUtility::boundVar(allOptions_, omega_, printToScreen);
         }
@@ -833,11 +829,7 @@ void DAkOmegaSSTFieldInversion::calcResiduals(const dictionary& options)
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverK = solve(kEqn);
-        if (printToScreen)
-        {
-            Info << "k Initial residual: " << solverK.initialResidual() << endl
-                 << "    Final residual: " << solverK.finalResidual() << endl;
-        }
+        DAUtility::primalResidualControl(solverK, printToScreen, "k");
 
         DAUtility::boundVar(allOptions_, k_, printToScreen);
 

diff --git a/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSSTLM.C b/src/adjoint/DAModel/DATurbulenceModel/DAkOmegaSSTLM.C
@@ -1105,11 +1105,7 @@ void DAkOmegaSSTLM::calcResiduals(const dictionary& options)
                 // get the solver performance info such as initial
                 // and final residuals
                 SolverPerformance<scalar> solverOmega = solve(omegaEqn);
-                if (printToScreen)
-                {
-                    Info << "omega Initial residual: " << solverOmega.initialResidual() << endl
-                         << "        Final residual: " << solverOmega.finalResidual() << endl;
-                }
+                DAUtility::primalResidualControl(solverOmega, printToScreen, "omega");
 
                 DAUtility::boundVar(allOptions_, omega_, printToScreen);
             }
@@ -1146,11 +1142,7 @@ void DAkOmegaSSTLM::calcResiduals(const dictionary& options)
             // get the solver performance info such as initial
             // and final residuals
             SolverPerformance<scalar> solverK = solve(kEqn);
-            if (printToScreen)
-            {
-                Info << "k Initial residual: " << solverK.initialResidual() << endl
-                     << "    Final residual: " << solverK.finalResidual() << endl;
-            }
+            DAUtility::primalResidualControl(solverK, printToScreen, "k");
 
             DAUtility::boundVar(allOptions_, k_, printToScreen);
 
@@ -1207,11 +1199,7 @@ void DAkOmegaSSTLM::calcResiduals(const dictionary& options)
                 // get the solver performance info such as initial
                 // and final residuals
                 SolverPerformance<scalar> solverReThetat = solve(ReThetatEqn);
-                if (printToScreen)
-                {
-                    Info << "ReThetat Initial residual: " << solverReThetat.initialResidual() << endl
-                         << "        Final residual: " << solverReThetat.finalResidual() << endl;
-                }
+                DAUtility::primalResidualControl(solverReThetat, printToScreen, "ReThetat");
 
                 DAUtility::boundVar(allOptions_, ReThetat_, printToScreen);
             }
@@ -1259,11 +1247,7 @@ void DAkOmegaSSTLM::calcResiduals(const dictionary& options)
                 // get the solver performance info such as initial
                 // and final residuals
                 SolverPerformance<scalar> solverGammaInt = solve(gammaIntEqn);
-                if (printToScreen)
-                {
-                    Info << "gammaInt Initial residual: " << solverGammaInt.initialResidual() << endl
-                         << "        Final residual: " << solverGammaInt.finalResidual() << endl;
-                }
+                DAUtility::primalResidualControl(solverGammaInt, printToScreen, "gammaInt");
 
                 DAUtility::boundVar(allOptions_, gammaInt_, printToScreen);
 

diff --git a/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C b/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C
@@ -82,12 +82,13 @@ label DAHeatTransferFoam::solvePrimal(
         return 1;
     }
 
-    primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printIntervalUnsteady");
     label printToScreen = 0;
     // main loop
     while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
     {
+        DAUtility::primalMaxInitRes_ = -1e16;
+
         printToScreen = this->isPrintTime(runTime, printInterval);
 
         if (printToScreen)
@@ -102,7 +103,7 @@ label DAHeatTransferFoam::solvePrimal(
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverT = TEqn.solve();
-        this->primalResidualControl<scalar>(solverT, printToScreen, printInterval, "T");
+        DAUtility::primalResidualControl(solverT, printToScreen, "T");
 
         if (this->validateStates())
         {

diff --git a/src/adjoint/DASolver/DALaplacianFoam/DALaplacianFoam.C b/src/adjoint/DASolver/DALaplacianFoam/DALaplacianFoam.C
@@ -145,14 +145,15 @@ label DALaplacianFoam::solvePrimal(
     // right after mesh.movePoints() calls.
     //mesh.moving(false);
 
-    primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printIntervalUnsteady");
     label printToScreen = 0;
     label timeInstanceI = 0;
 
     // main loop
     while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
     {
+        DAUtility::primalMaxInitRes_ = -1e16;
+
         printToScreen = this->isPrintTime(runTime, printInterval);
 
         if (printToScreen)
@@ -167,7 +168,7 @@ label DALaplacianFoam::solvePrimal(
         // get the solver performance info such as initial
         // and final residuals
         SolverPerformance<scalar> solverT = TEqn.solve();
-        this->primalResidualControl<scalar>(solverT, printToScreen, printInterval, "T");
+        DAUtility::primalResidualControl(solverT, printToScreen, "T");
 
         if (this->validateStates())
         {

diff --git a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C
@@ -154,7 +154,6 @@ label DAPimpleFoam::solvePrimal(
     // right after mesh.movePoints() calls.
     //mesh.moving(false);
 
-    primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printIntervalUnsteady");
     label printToScreen = 0;
     label timeInstanceI = 0;

diff --git a/src/adjoint/DASolver/DAPimpleFoam/UEqnPimple.H b/src/adjoint/DASolver/DAPimpleFoam/UEqnPimple.H
@@ -22,7 +22,7 @@ if (pimple.momentumPredictor())
     // and final residuals
     SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-    this->primalResidualControl<vector>(solverU, pimplePrintToScreen, printInterval, "U");
+    DAUtility::primalResidualControl(solverU, pimplePrintToScreen, "U");
 }
 
 // bound U

diff --git a/src/adjoint/DASolver/DAPimpleFoam/pEqnPimple.H b/src/adjoint/DASolver/DAPimpleFoam/pEqnPimple.H
@@ -52,7 +52,7 @@ while (pimple.correctNonOrthogonal())
     // and final residuals
     SolverPerformance<scalar> solverP = pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-    this->primalResidualControl<scalar>(solverP, pimplePrintToScreen, printInterval, "p");
+    DAUtility::primalResidualControl(solverP, pimplePrintToScreen, "p");
 
     if (pimple.finalNonOrthogonalIter())
     {

diff --git a/src/adjoint/DASolver/DAPisoFoam/DAPisoFoam.C b/src/adjoint/DASolver/DAPisoFoam/DAPisoFoam.C
@@ -169,14 +169,15 @@ label DAPisoFoam::solvePrimal(
     // right after mesh.movePoints() calls.
     //mesh.moving(false);
 
-    primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printIntervalUnsteady");
     label printToScreen = 0;
     label timeInstanceI = 0;
 
     // main loop
     while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
     {
+        DAUtility::primalMaxInitRes_ = -1e16;
+
         printToScreen = this->isPrintTime(runTime, printInterval);
 
         if (printToScreen)

diff --git a/src/adjoint/DASolver/DAPisoFoam/UEqnPiso.H b/src/adjoint/DASolver/DAPisoFoam/UEqnPiso.H
@@ -20,7 +20,7 @@ if (piso.momentumPredictor())
     // and final residuals
     SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-    this->primalResidualControl<vector>(solverU, printToScreen, printInterval, "U");
+    DAUtility::primalResidualControl(solverU, printToScreen, "U");
 }
 
 // bound U

diff --git a/src/adjoint/DASolver/DAPisoFoam/pEqnPiso.H b/src/adjoint/DASolver/DAPisoFoam/pEqnPiso.H
@@ -21,7 +21,7 @@ while (piso.correctNonOrthogonal())
     // and final residuals
     SolverPerformance<scalar> solverP = pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-    this->primalResidualControl<scalar>(solverP, printToScreen, printInterval, "p");
+    DAUtility::primalResidualControl(solverP, printToScreen, "p");
 
     if (piso.finalNonOrthogonalIter())
     {

diff --git a/src/adjoint/DASolver/DARhoSimpleCFoam/DARhoSimpleCFoam.C b/src/adjoint/DASolver/DARhoSimpleCFoam/DARhoSimpleCFoam.C
@@ -132,12 +132,13 @@ label DARhoSimpleCFoam::solvePrimal(
     // check if the parameters are set in the Python layer
     daRegressionPtr_->validate();
 
-    primalMinRes_ = 1e10;
+
     label printInterval = daOptionPtr_->getOption<label>("printInterval");
     label printToScreen = 0;
     label regModelFail = 0;
     while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
     {
+        DAUtility::primalMaxInitRes_ = -1e16;
 
         printToScreen = this->isPrintTime(runTime, printInterval);
 

diff --git a/src/adjoint/DASolver/DARhoSimpleCFoam/EEqnRhoSimpleC.H b/src/adjoint/DASolver/DARhoSimpleCFoam/EEqnRhoSimpleC.H
@@ -14,7 +14,7 @@
     // and final residuals
     SolverPerformance<scalar> solverE = EEqn.solve();
 
-    this->primalResidualControl<scalar>(solverE, printToScreen, printInterval, "he");
+    DAUtility::primalResidualControl(solverE, printToScreen, "he");
 
     // bound he
     DAUtility::boundVar(allOptions, he, printToScreen);