diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
index 75f65b12..b1d93227 100755
--- 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
index 7ac2d0eb..8b0782c3 100755
--- 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
index 4c819f9f..afe68f68 100755
--- 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
index 5e2f545e..0bd65416 100755
--- 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
index 7731d30f..3e9d7994 100755
--- 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
index 272b5f30..ffeeef65 100755
--- 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
index b9d1cb9c..acfa006e 100755
--- 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
index 330ad573..ee9651fb 100755
--- 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
index 949c82d8..dc47e2d2 100755
--- 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
index bc7a9177..1e92966a 100755
--- 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
index 21ac75e4..43c8b87a 100755
--- 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
index 3658c349..bd6595d6 100755
--- 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
index b1c99584..4d20504a 100755
--- a/src/adjoint/DASolver/DALaplacianFoam/DALaplacianFoam.C
+++ b/src/adjoint/DASolver/DALaplacianFoam/DALaplacianFoam.C
@@ -145,7 +145,6 @@ 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;
@@ -153,6 +152,8 @@ label DALaplacianFoam::solvePrimal(
     // 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
index a6b63c88..1507f456 100755
--- 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
index bc2064a5..eb29677f 100755
--- 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
index fd212002..bfab46d8 100755
--- 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
index 76200097..1e67c9e9 100755
--- a/src/adjoint/DASolver/DAPisoFoam/DAPisoFoam.C
+++ b/src/adjoint/DASolver/DAPisoFoam/DAPisoFoam.C
@@ -169,7 +169,6 @@ 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;
@@ -177,6 +176,8 @@ label DAPisoFoam::solvePrimal(
     // 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
index a0447700..5bd73fa4 100755
--- 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
index 2f2e4560..e00e5816 100755
--- 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
index c0b8809f..75e50f5e 100755
--- 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
index b691ce77..35fce84e 100755
--- 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);
diff --git a/src/adjoint/DASolver/DARhoSimpleCFoam/UEqnRhoSimpleC.H b/src/adjoint/DASolver/DARhoSimpleCFoam/UEqnRhoSimpleC.H
index 50881137..ad3a88ba 100755
--- a/src/adjoint/DASolver/DARhoSimpleCFoam/UEqnRhoSimpleC.H
+++ b/src/adjoint/DASolver/DARhoSimpleCFoam/UEqnRhoSimpleC.H
@@ -11,7 +11,7 @@ UEqn.relax();
 // 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
 DAUtility::boundVar(allOptions, U, printToScreen);
diff --git a/src/adjoint/DASolver/DARhoSimpleCFoam/pEqnRhoSimpleC.H b/src/adjoint/DASolver/DARhoSimpleCFoam/pEqnRhoSimpleC.H
index 31445230..0e973f3e 100755
--- a/src/adjoint/DASolver/DARhoSimpleCFoam/pEqnRhoSimpleC.H
+++ b/src/adjoint/DASolver/DARhoSimpleCFoam/pEqnRhoSimpleC.H
@@ -57,7 +57,7 @@ while (simple.correctNonOrthogonal())
     // and final residuals
     SolverPerformance<scalar> solverP = pEqn.solve();
 
-    this->primalResidualControl<scalar>(solverP, printToScreen, printInterval, "p");
+    DAUtility::primalResidualControl(solverP, printToScreen, "p");
 
     if (simple.finalNonOrthogonalIter())
     {
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
index e3f10f3a..c27d761a 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
@@ -149,12 +149,12 @@ label DARhoSimpleFoam::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/DARhoSimpleFoam/EEqnRhoSimple.H b/src/adjoint/DASolver/DARhoSimpleFoam/EEqnRhoSimple.H
index e5ac7579..f17ac6de 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/EEqnRhoSimple.H
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/EEqnRhoSimple.H
@@ -21,7 +21,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);
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H b/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H
index 1f5fea8d..10c04ba7 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H
@@ -23,7 +23,7 @@ fvOptions.constrain(UEqn);
 // and final residuals
 SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-this->primalResidualControl<vector>(solverU, printToScreen, printInterval, "U");
+DAUtility::primalResidualControl(solverU, printToScreen, "U");
 
 fvOptions.correct(U);
 
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H b/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H
index 6c0870e1..e100626e 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H
@@ -47,7 +47,7 @@ while (simple.correctNonOrthogonal())
     // and final residuals
     SolverPerformance<scalar> solverP = pEqn.solve();
 
-    this->primalResidualControl<scalar>(solverP, printToScreen, printInterval, "p");
+    DAUtility::primalResidualControl(solverP, printToScreen, "p");
 
     if (simple.finalNonOrthogonalIter())
     {
diff --git a/src/adjoint/DASolver/DAScalarTransportFoam/DAScalarTransportFoam.C b/src/adjoint/DASolver/DAScalarTransportFoam/DAScalarTransportFoam.C
index d2947573..a93b23db 100755
--- a/src/adjoint/DASolver/DAScalarTransportFoam/DAScalarTransportFoam.C
+++ b/src/adjoint/DASolver/DAScalarTransportFoam/DAScalarTransportFoam.C
@@ -147,7 +147,6 @@ label DAScalarTransportFoam::solvePrimal(
     // right after mesh.movePoints() calls.
     //mesh.moving(false);
 
-    primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printIntervalUnsteady");
     label printToScreen = 0;
     label timeInstanceI = 0;
@@ -155,6 +154,8 @@ label DAScalarTransportFoam::solvePrimal(
     // 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)
@@ -172,7 +173,7 @@ label DAScalarTransportFoam::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/DASimpleFoam/DASimpleFoam.C b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
index 3528f39f..c90ae46d 100755
--- a/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
+++ b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
@@ -160,12 +160,12 @@ label DASimpleFoam::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/DASimpleFoam/UEqnSimple.H b/src/adjoint/DASolver/DASimpleFoam/UEqnSimple.H
index 9d171b17..c5449737 100755
--- a/src/adjoint/DASolver/DASimpleFoam/UEqnSimple.H
+++ b/src/adjoint/DASolver/DASimpleFoam/UEqnSimple.H
@@ -26,7 +26,7 @@ if (simple.momentumPredictor())
     // and final residuals
     SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-    this->primalResidualControl<vector>(solverU, printToScreen, printInterval, "U");
+    DAUtility::primalResidualControl(solverU, printToScreen, "U");
 
     fvOptions.correct(U);
 }
diff --git a/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H b/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H
index 69b8161f..c18a5f5c 100755
--- a/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H
+++ b/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H
@@ -49,7 +49,7 @@
         // and final residuals
         SolverPerformance<scalar> solverP = pEqn.solve();
 
-        this->primalResidualControl<scalar>(solverP, printToScreen, printInterval, "p");
+        DAUtility::primalResidualControl(solverP, printToScreen, "p");
 
         if (simple.finalNonOrthogonalIter())
         {
diff --git a/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C b/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C
index b8b2e78c..23175d90 100755
--- a/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C
+++ b/src/adjoint/DASolver/DASimpleTFoam/DASimpleTFoam.C
@@ -139,12 +139,12 @@ label DASimpleTFoam::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/DASimpleTFoam/TEqnSimpleT.H b/src/adjoint/DASolver/DASimpleTFoam/TEqnSimpleT.H
index 682eedca..493ace16 100644
--- a/src/adjoint/DASolver/DASimpleTFoam/TEqnSimpleT.H
+++ b/src/adjoint/DASolver/DASimpleTFoam/TEqnSimpleT.H
@@ -14,7 +14,7 @@
     // and final residuals
     SolverPerformance<scalar> solverT = TEqn.solve();
 
-    this->primalResidualControl<scalar>(solverT, printToScreen, printInterval, "T");
+    DAUtility::primalResidualControl(solverT, printToScreen, "T");
 
     // bound T
     DAUtility::boundVar(allOptions, T, printToScreen);
diff --git a/src/adjoint/DASolver/DASimpleTFoam/UEqnSimpleT.H b/src/adjoint/DASolver/DASimpleTFoam/UEqnSimpleT.H
index cb3d468e..6787487b 100755
--- a/src/adjoint/DASolver/DASimpleTFoam/UEqnSimpleT.H
+++ b/src/adjoint/DASolver/DASimpleTFoam/UEqnSimpleT.H
@@ -24,7 +24,7 @@ if (simple.momentumPredictor())
     // and final residuals
     SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-    this->primalResidualControl<vector>(solverU, printToScreen, printInterval, "U");
+    DAUtility::primalResidualControl(solverU, printToScreen, "U");
 
     fvOptions.correct(U);
 }
diff --git a/src/adjoint/DASolver/DASimpleTFoam/pEqnSimpleT.H b/src/adjoint/DASolver/DASimpleTFoam/pEqnSimpleT.H
index 69b8161f..c18a5f5c 100755
--- a/src/adjoint/DASolver/DASimpleTFoam/pEqnSimpleT.H
+++ b/src/adjoint/DASolver/DASimpleTFoam/pEqnSimpleT.H
@@ -49,7 +49,7 @@
         // and final residuals
         SolverPerformance<scalar> solverP = pEqn.solve();
 
-        this->primalResidualControl<scalar>(solverP, printToScreen, printInterval, "p");
+        DAUtility::primalResidualControl(solverP, printToScreen, "p");
 
         if (simple.finalNonOrthogonalIter())
         {
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C b/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C
index 21f18f3d..7cf4b7a1 100755
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C
+++ b/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C
@@ -108,11 +108,11 @@ label DASolidDisplacementFoam::solvePrimal(
         return 1;
     }
 
-    primalMinRes_ = 1e10;
     label printInterval = daOptionPtr_->getOption<label>("printInterval");
     label printToScreen = 0;
     while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
     {
+        DAUtility::primalMaxInitRes_ = -1e16;
 
         printToScreen = this->isPrintTime(runTime, printInterval);
 
@@ -135,7 +135,7 @@ label DASolidDisplacementFoam::solvePrimal(
         // and final residuals
         SolverPerformance<vector> solverD = DEqn.solve();
 
-        this->primalResidualControl<vector>(solverD, printToScreen, printInterval, "D");
+        DAUtility::primalResidualControl(solverD, printToScreen, "D");
 
         if (this->validateStates())
         {
diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index 77443e2e..cdd6d8fa 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -19,6 +19,8 @@ pyComputeInterface Foam::DAUtility::pyCalcBetaInterface = NULL;
 void* Foam::DAUtility::pyCalcBetaJacVecProd = NULL;
 pyJacVecProdInterface Foam::DAUtility::pyCalcBetaJacVecProdInterface = NULL;
 
+scalar Foam::DAUtility::primalMaxInitRes_ = -1e16;
+
 namespace Foam
 {
 
@@ -134,10 +136,10 @@ label DASolver::loop(Time& runTime)
     }
 
     // check exit condition
-    if (primalMinRes_ < primalMinResTol_ && runTime.timeIndex() > primalMinIters_)
+    if (DAUtility::primalMaxInitRes_ < primalMinResTol_ && runTime.timeIndex() > primalMinIters_)
     {
         Info << "Time = " << t << endl;
-        Info << "Minimal residual " << primalMinRes_ << " satisfied the prescribed tolerance " << primalMinResTol_ << endl
+        Info << "Minimal residual " << DAUtility::primalMaxInitRes_ << " satisfied the prescribed tolerance " << primalMinResTol_ << endl
              << endl;
         this->printAllObjFuncs();
         runTime.writeNow();
@@ -8481,10 +8483,10 @@ label DASolver::checkResidualTol()
 
     scalar tol = daOptionPtr_->getOption<scalar>("primalMinResTol");
     scalar tolMax = daOptionPtr_->getOption<scalar>("primalMinResTolDiff");
-    if (primalMinRes_ / tol > tolMax)
+    if (DAUtility::primalMaxInitRes_ / tol > tolMax)
     {
         Info << "********************************************" << endl;
-        Info << "Primal min residual " << primalMinRes_ << endl
+        Info << "Primal min residual " << DAUtility::primalMaxInitRes_ << endl
              << "did not satisfy the prescribed tolerance "
              << tol << endl;
         Info << "Primal solution failed!" << endl;
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index b3995f3e..315954d7 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -139,20 +139,9 @@ protected:
     /// the derivative value computed by the forward mode primal solution.
     HashTable<PetscScalar> forwardADDerivVal_;
 
-    /// the maximal residual for primal solution
-    scalar primalMinRes_ = 1e10;
-
     /// the solution time for the previous primal solution
     scalar prevPrimalSolTime_ = -1e10;
 
-    /// setup maximal residual control and print the residual as needed
-    template<class classType>
-    void primalResidualControl(
-        const SolverPerformance<classType>& solverP,
-        const label printToScreen,
-        const label printInterval,
-        const word varName);
-
     label isPrintTime(
         const Time& runTime,
         const label printInterval) const;
@@ -1286,28 +1275,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-template<class classType>
-void DASolver::primalResidualControl(
-    const SolverPerformance<classType>& solverP,
-    const label printToScreen,
-    const label printInterval,
-    const word varName)
-{
-    /*
-    Description:
-        Setup maximal residual control and print the residual as needed
-    */
-    if (mag(solverP.initialResidual()) < primalMinRes_)
-    {
-        primalMinRes_ = mag(solverP.initialResidual());
-    }
-    if (printToScreen)
-    {
-        Info << varName << " Initial residual: " << solverP.initialResidual() << endl
-             << varName << "   Final residual: " << solverP.finalResidual() << endl;
-    }
-}
-
 template<class classType>
 label DASolver::validateField(const classType& field)
 {
diff --git a/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C b/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C
index 3ed8d3a2..400dbf8f 100755
--- a/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C
+++ b/src/adjoint/DASolver/DATurboFoam/DATurboFoam.C
@@ -102,12 +102,12 @@ label DATurboFoam::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/DATurboFoam/EEqnTurbo.H b/src/adjoint/DASolver/DATurboFoam/EEqnTurbo.H
index 340ce801..b01b71ab 100755
--- a/src/adjoint/DASolver/DATurboFoam/EEqnTurbo.H
+++ b/src/adjoint/DASolver/DATurboFoam/EEqnTurbo.H
@@ -19,7 +19,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);
diff --git a/src/adjoint/DASolver/DATurboFoam/UEqnTurbo.H b/src/adjoint/DASolver/DATurboFoam/UEqnTurbo.H
index a4bc8394..545aebd5 100755
--- a/src/adjoint/DASolver/DATurboFoam/UEqnTurbo.H
+++ b/src/adjoint/DASolver/DATurboFoam/UEqnTurbo.H
@@ -14,7 +14,7 @@ UEqn.relax();
 // 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
 DAUtility::boundVar(allOptions, U, printToScreen);
diff --git a/src/adjoint/DASolver/DATurboFoam/pEqnTurbo.H b/src/adjoint/DASolver/DATurboFoam/pEqnTurbo.H
index ffc76af7..28694f09 100755
--- a/src/adjoint/DASolver/DATurboFoam/pEqnTurbo.H
+++ b/src/adjoint/DASolver/DATurboFoam/pEqnTurbo.H
@@ -45,7 +45,7 @@ if (simple.transonic())
         // and final residuals
         SolverPerformance<scalar> solverP = pEqn.solve();
 
-        this->primalResidualControl<scalar>(solverP, printToScreen, printInterval, "p");
+        DAUtility::primalResidualControl(solverP, printToScreen, "p");
 
         if (simple.finalNonOrthogonalIter())
         {
@@ -78,7 +78,7 @@ else
         // and final residuals
         SolverPerformance<scalar> solverP = pEqn.solve();
 
-        this->primalResidualControl<scalar>(solverP, printToScreen, printInterval, "p");
+        DAUtility::primalResidualControl(solverP, printToScreen, "p");
 
         if (simple.finalNonOrthogonalIter())
         {
diff --git a/src/adjoint/DAUtility/DAUtility.C b/src/adjoint/DAUtility/DAUtility.C
index 0294bc5f..423740e6 100755
--- a/src/adjoint/DAUtility/DAUtility.C
+++ b/src/adjoint/DAUtility/DAUtility.C
@@ -731,6 +731,67 @@ label DAUtility::isValueCloseToRef(
     }
 }
 
+void DAUtility::primalResidualControl(
+    const SolverPerformance<scalar>& solverP,
+    const label printToScreen,
+    const word varName)
+{
+    /*
+    Description:
+        Setup maximal residual control and print the residual as needed
+    */
+
+    // calculate the initial residual mag and set it to primalResidualNorms_
+
+    scalar initRes = solverP.initialResidual();
+
+    if (initRes > DAUtility::primalMaxInitRes_)
+    {
+        DAUtility::primalMaxInitRes_ = initRes;
+    }
+
+    if (printToScreen)
+    {
+        Info << varName << " Initial residual: " << solverP.initialResidual() << endl
+             << varName << "   Final residual: " << solverP.finalResidual() << endl;
+    }
+}
+
+void DAUtility::primalResidualControl(
+    const SolverPerformance<vector>& solverP,
+    const label printToScreen,
+    const word varName)
+{
+    /*
+    Description:
+        Setup maximal residual control and print the residual as needed
+    */
+
+    // calculate the initial residual mag and set it to primalResidualNorms_
+
+    // for vectors, we need to use the median value for the residual
+    // this is because we often need to run 2D simulations with symmetry 
+    // BC, so one component of the residual vector, which is related to the symmetry BC,
+    // may be high while the other two components' residuals are low.
+    // In this case, we can use the median value for the residual vector, which better
+    // represents the convergence.
+
+    vector initRes = solverP.initialResidual();
+    scalarList initResList = {initRes[0], initRes[1], initRes[2]};
+    sort(initResList);
+
+    if (initResList[1] > DAUtility::primalMaxInitRes_)
+    {
+        DAUtility::primalMaxInitRes_ = initResList[1];
+    }
+
+    if (printToScreen)
+    {
+        Info << varName << " Initial residual: " << solverP.initialResidual() << endl
+             << varName << "   Final residual: " << solverP.finalResidual() << endl;
+    }
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/adjoint/DAUtility/DAUtility.H b/src/adjoint/DAUtility/DAUtility.H
index 6f5e6963..f1463756 100755
--- a/src/adjoint/DAUtility/DAUtility.H
+++ b/src/adjoint/DAUtility/DAUtility.H
@@ -131,6 +131,21 @@ public:
 
     static void* pyCalcBetaJacVecProd;
     static pyJacVecProdInterface pyCalcBetaJacVecProdInterface;
+
+    /// the max initial residual norms for the primal solution for each state
+    static scalar primalMaxInitRes_;
+
+    /// control when to print the residual and also compute the maxInitRes
+    static void primalResidualControl(
+        const SolverPerformance<scalar>& solverP,
+        const label printToScreen,
+        const word varName);
+    
+    /// control when to print the residual and also compute the maxInitRes
+    static void primalResidualControl(
+        const SolverPerformance<vector>& solverP,
+        const label printToScreen,
+        const word varName);
 };
 
 template<class classType>
diff --git a/tests/refs/DAFoam_Test_DARhoSimpleCFoamRef.txt b/tests/refs/DAFoam_Test_DARhoSimpleCFoamRef.txt
index 64e237ed..dd4a61a8 100755
--- a/tests/refs/DAFoam_Test_DARhoSimpleCFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DARhoSimpleCFoamRef.txt
@@ -1,22 +1,22 @@
 Dictionary Key: CD
-@value  0.01016659179082796 1e-08 1e-10
+@value  0.01016659179618512 1e-08 1e-10
 Dictionary Key: CL
-@value   0.4536846351797039 1e-08 1e-10
+@value   0.4536846350009314 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
 Dictionary Key: CD
 Dictionary Key: alpha
-@value 0.004336901757235039 0.0001 1e-06
+@value 0.004419288619714606 0.0001 1e-06
 Dictionary Key: shapey
-@value   0.0267621536016065 0.0001 1e-06
-@value  0.04958022286719843 0.0001 1e-06
-@value   0.1319330548362239 0.0001 1e-06
+@value  0.02619137460608506 0.0001 1e-06
+@value  0.05002489637978214 0.0001 1e-06
+@value   0.1314687490333257 0.0001 1e-06
 Dictionary Key: CL
 Dictionary Key: alpha
-@value   0.1540391990713826 0.0001 1e-06
+@value   0.1586052418041611 0.0001 1e-06
 Dictionary Key: shapey
-@value    3.621966829337841 0.0001 1e-06
-@value     2.44543902410904 0.0001 1e-06
-@value    3.016228732003395 0.0001 1e-06
+@value     3.58362402336514 0.0001 1e-06
+@value    2.456281547160383 0.0001 1e-06
+@value     2.99368621643047 0.0001 1e-06
 Dictionary Key: fail
 @value                    0 0.0001 1e-06
\ No newline at end of file
diff --git a/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFIMLRef.txt b/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFIMLRef.txt
index eec87ccf..23c9553f 100755
--- a/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFIMLRef.txt
+++ b/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFIMLRef.txt
@@ -1,14 +1,14 @@
 Dictionary Key: CD
-@value 0.006024441240321594 1e-06 1e-10
+@value 0.006026347733544404 1e-06 1e-10
 Dictionary Key: CL
-@value   0.0879110908067339 1e-06 1e-10
+@value  0.08789446345115981 1e-06 1e-10
 Dictionary Key: fail
 @value                    0 1e-06 1e-10
 Dictionary Key: CD
 Dictionary Key: alpha
-@value 0.000427152282947212 0.001 1e-06
+@value 0.0004276629226289779 0.001 1e-06
 Dictionary Key: CL
 Dictionary Key: alpha
-@value    0.102641313442945 0.001 1e-06
+@value   0.1026067641371463 0.001 1e-06
 Dictionary Key: fail
 @value                    0 0.001 1e-06
\ No newline at end of file
diff --git a/tests/refs/DAFoam_Test_DATurboFoamSubsonicRef.txt b/tests/refs/DAFoam_Test_DATurboFoamSubsonicRef.txt
index cbbaad4e..d81218bb 100755
--- a/tests/refs/DAFoam_Test_DATurboFoamSubsonicRef.txt
+++ b/tests/refs/DAFoam_Test_DATurboFoamSubsonicRef.txt
@@ -1,30 +1,30 @@
 Dictionary Key: CMZ
-@value   0.1430241062895427 1e-08 1e-10
+@value   0.1430241062697062 1e-08 1e-10
 Dictionary Key: TPR
-@value    1.027591736890468 1e-08 1e-10
+@value    1.027591736884627 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
 Dictionary Key: CMZ
 Dictionary Key: shapey
-@value   0.2847715919611993 1e-04 1e-06
-@value   0.5020122841985006 1e-04 1e-06
-@value    0.339989430432045 1e-04 1e-06
-@value   0.4471545417329905 1e-04 1e-06
+@value    0.284774235720788 0.0001 1e-06
+@value   0.5019992820686592 0.0001 1e-06
+@value   0.3399894733398309 0.0001 1e-06
+@value   0.4471495481357648 0.0001 1e-06
 Dictionary Key: shapez
-@value -0.03754549250758193 1e-04 1e-06
-@value  -0.1668372689661463 1e-04 1e-06
-@value  -0.1865409331445353 1e-04 1e-06
-@value  -0.2092364378494078 1e-04 1e-06
+@value -0.03754471534814914 0.0001 1e-06
+@value  -0.1668261693510449 0.0001 1e-06
+@value  -0.1865385648670614 0.0001 1e-06
+@value  -0.2092321326526177 0.0001 1e-06
 Dictionary Key: TPR
 Dictionary Key: shapey
-@value  0.08984549038801276 0.0001 1e-06
-@value  -0.1026120377600116 0.0001 1e-06
-@value  0.07889370056517067 0.0001 1e-06
-@value  -0.0846173656569997 0.0001 1e-06
+@value  0.08984231476276661 0.0001 1e-06
+@value  -0.1026282832287896 0.0001 1e-06
+@value  0.07889381125344438 0.0001 1e-06
+@value -0.08462041661143266 0.0001 1e-06
 Dictionary Key: shapez
-@value -0.03253030126355309 0.0001 1e-06
-@value  0.09548883025920442 0.0001 1e-06
-@value  -0.0452150125277994 0.0001 1e-06
-@value   0.1037784570108749 0.0001 1e-06
+@value -0.03252865833279311 0.0001 1e-06
+@value  0.09550057053737057 0.0001 1e-06
+@value -0.04521482413177378 0.0001 1e-06
+@value    0.103781649612844 0.0001 1e-06
 Dictionary Key: fail
-@value                    0 1e-04 1e-06
+@value                    0 0.0001 1e-06
diff --git a/tests/runTests_DARhoSimpleCFoam.py b/tests/runTests_DARhoSimpleCFoam.py
index f81dd960..af569702 100755
--- a/tests/runTests_DARhoSimpleCFoam.py
+++ b/tests/runTests_DARhoSimpleCFoam.py
@@ -43,7 +43,7 @@
     "solverName": "DARhoSimpleCFoam",
     "useAD": {"mode": "fd"},
     "designSurfaces": ["wing"],
-    "primalMinResTol": 1e-12,
+    "primalMinResTol": 1e-11,
     "writeJacobians": ["all"],
     "primalBC": {
         "UIn": {"variable": "U", "patches": ["inout"], "value": [U0, 0.0, 0.0]},
diff --git a/tests/runTests_DARhoSimpleCFoamAD.py b/tests/runTests_DARhoSimpleCFoamAD.py
index a5732734..9b5cae03 100755
--- a/tests/runTests_DARhoSimpleCFoamAD.py
+++ b/tests/runTests_DARhoSimpleCFoamAD.py
@@ -44,7 +44,7 @@
     "designSurfaces": ["wing"],
     "useAD": {"mode": "reverse"},
     "writeJacobians": ["all"],
-    "primalMinResTol": 1e-12,
+    "primalMinResTol": 3e-12,
     "primalBC": {
         "UIn": {"variable": "U", "patches": ["inout"], "value": [U0, 0.0, 0.0]},
         "p0": {"variable": "p", "patches": ["inout"], "value": [p0]},
diff --git a/tests/runTests_DARhoSimpleFoamUBend.py b/tests/runTests_DARhoSimpleFoamUBend.py
index 197725c6..a9047782 100755
--- a/tests/runTests_DARhoSimpleFoamUBend.py
+++ b/tests/runTests_DARhoSimpleFoamUBend.py
@@ -37,6 +37,7 @@
     "useAD": {"mode": "fd"},
     "designSurfaces": ["ubend"],
     "primalMinResTol": 1e-5,
+    "primalMinResTolDiff": 1e5,
     "primalBC": {
         "U0": {"variable": "U", "patches": ["inlet"], "value": [U0, 0.0, 0.0]},
         "p0": {"variable": "p", "patches": ["outlet"], "value": [p0]},