From f33075ad67aab8b100fb7efd29f7ab38a9b1918a Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Thu, 20 Jul 2023 16:07:38 -0500
Subject: [PATCH] Added the deactivate calls for all AD funcs and updated
 tests.

---
 src/adjoint/DASolver/DASolver.C               | 954 ++++++++++++------
 src/adjoint/DASolver/DASolver.H               |  22 +-
 ...SimpleFoamkOmegaFieldInversionOmegaRef.txt |   2 +-
 ...DASimpleFoamkOmegaSSTFieldInversionRef.txt |   2 +-
 tests/refs/DAFoam_Test_MphysAeroFieldRef.txt  |  29 +-
 tests/refs/DAFoam_Test_MphysAeroOptRef.txt    |   4 +-
 .../DAFoam_Test_MphysAeroPropCoupledRef.txt   | 148 +--
 tests/refs/DAFoam_Test_MphysAeroPropRef.txt   | 176 ++--
 tests/runTests_MphysAeroField.py              |  27 +-
 tests/runTests_MphysAeroOpt.py                |   4 +-
 10 files changed, 864 insertions(+), 504 deletions(-)

diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index c121b4d1..2c4807a0 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -475,13 +475,22 @@ void DASolver::calcCouplingFaceCoordsAD(
         product[i] = volCoordsArray[i].getGradient();
     }
 
-    delete[] volCoordsArray;
-    delete[] surfCoordsArray;
-
     // clean up AD
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
 
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+    for (label i = 0; i < daIndexPtr_->nLocalXv; i++)
+    {
+        this->globalADTape_.deactivateValue(volCoordsArray[i]);
+    }
+    this->calcCouplingFaceCoords(volCoordsArray, surfCoordsArray);
+
+    delete[] volCoordsArray;
+    delete[] surfCoordsArray;
+
 #endif
 }
 
@@ -1005,26 +1014,21 @@ void DASolver::getThermalAD(
                                      << abort(FatalError);
     }
 
-    // ******************************* NOTE! ******************************
-    // This step is tricky. We need to clean up the AD seeds for all variables before the next AD
-    // This is usually done by theses two calls: updateOFField(wVec) and updateOFMesh(xvVec)
-    // The above two calls will assign zeros to the gradient() part of all scalar variables because
-    // the wVec and xvVec's gradient() part is zero (they are double instead of scalar).
-    // HOWEVER, the above two calls will NOT clean the seeds for field variables's boundary values
-    // This will not cause problems for most of the functions because their boundary values always
-    // have zero seeds. However, the getThermal function actually propagates non-zeros seeds for boundary
-    // values, so here we need to manually clean up the BC's gradient part. NOT doing this will
-    // mess up the AD computation for the following function.
-    // The way we clean it is that we set a zero-seed thermal var and call getThermal to propagate
-    // the zero seeds to all intermediate vars.
-    // ******************************* NOTE! ******************************
+    // clean up AD
+    this->globalADTape_.clearAdjoints();
+    this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
     for (label i = 0; i < daIndexPtr_->nLocalXv; i++)
     {
-        volCoordsArray[i].setGradient(0.0);
+        this->globalADTape_.deactivateValue(volCoordsArray[i]);
     }
     for (label i = 0; i < daIndexPtr_->nLocalAdjointStates; i++)
     {
-        statesArray[i].setGradient(0.0);
+        this->globalADTape_.deactivateValue(statesArray[i]);
     }
     this->getThermal(volCoordsArray, statesArray, thermalArray);
 
@@ -1032,10 +1036,6 @@ void DASolver::getThermalAD(
     delete[] statesArray;
     delete[] thermalArray;
 
-    // clean up AD
-    this->globalADTape_.clearAdjoints();
-    this->globalADTape_.reset();
-
 #endif
 }
 
@@ -1869,10 +1869,7 @@ void DASolver::calcdForceProfiledXvWAD(
         FatalErrorIn("calcdFvSourcedInputsTPsiAD") << "inputMode not valid"
                                                    << abort(FatalError);
     }
-    daResidualPtr_->correctBoundaryConditions();
-    daResidualPtr_->updateIntermediateVariables();
-    daModelPtr_->correctBoundaryConditions();
-    daModelPtr_->updateIntermediateVariables();
+    this->updateStateBoundaryConditions();
 
     // Step 3
     this->calcForceProfileInternal(propName, aForce, tForce, rDist, integralForce);
@@ -1972,6 +1969,29 @@ void DASolver::calcdForceProfiledXvWAD(
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
 
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    if (inputMode == "mesh")
+    {
+        forAll(meshPoints, i)
+        {
+            for (label j = 0; j < 3; j++)
+            {
+                this->globalADTape_.deactivateValue(meshPoints[i][j]);
+            }
+        }
+        meshPtr_->movePoints(meshPoints);
+        meshPtr_->moving(false);
+    }
+    else if (inputMode == "state")
+    {
+        this->deactivateStateVariableInput4AD();
+    }
+    this->updateStateBoundaryConditions();
+    this->calcForceProfileInternal(propName, aForce, tForce, rDist, integralForce);
+
 #endif
 }
 
@@ -2640,6 +2660,48 @@ void DASolver::calcdFvSourcedInputsTPsiAD(
 
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    if (mode == "aForce")
+    {
+        for (label i = 0; i < nPoints; i++)
+        {
+            this->globalADTape_.deactivateValue(aForceField[i]);
+        }
+    }
+    else if (mode == "tForce")
+    {
+        for (label i = 0; i < nPoints; i++)
+        {
+            this->globalADTape_.deactivateValue(tForceField[i]);
+        }
+    }
+    else if (mode == "rDist")
+    {
+        for (label i = 0; i < nPoints; i++)
+        {
+            this->globalADTape_.deactivateValue(rDistField[i]);
+        }
+    }
+    else if (mode == "targetForce")
+    {
+        for (label i = 0; i < 2; i++)
+        {
+            this->globalADTape_.deactivateValue(targetForceList[i]);
+        }
+    }
+    else if (mode == "center")
+    {
+        for (label i = 0; i < 3; i++)
+        {
+            this->globalADTape_.deactivateValue(centerVector[i]);
+        }
+    }
+
+    this->calcFvSourceInternal(propName, aForceField, tForceField, rDistField, targetForceList, centerVector, fvSourceVField);
 #endif
 }
 
@@ -3537,6 +3599,159 @@ void DASolver::calcdRdAOA(
     }
 }
 
+void DASolver::setBCToOFVars(
+    const dictionary& dvSubDict,
+    const scalar& BC)
+{
+    /*
+    Description:
+        Assign the boundary condition (BC) value to OF variables
+    
+    Input:
+
+        dvSubDict: the design variable subDict that contains the inform for setting BCs
+
+        BC: the BC value
+    */
+
+    // get info from dvSubDict. This needs to be defined in the pyDAFoam
+    // name of the variable for changing the boundary condition
+    word varName = dvSubDict.getWord("variable");
+    // name of the boundary patch
+    wordList patches;
+    dvSubDict.readEntry<wordList>("patches", patches);
+    // the component of a vector variable, ignore when it is a scalar
+    label comp = dvSubDict.getLabel("comp");
+
+    // set the BC
+    forAll(patches, idxI)
+    {
+        word patchName = patches[idxI];
+        label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
+        if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
+        {
+            volVectorField& state(const_cast<volVectorField&>(
+                meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
+            // for decomposed domain, don't set BC if the patch is empty
+            if (meshPtr_->boundaryMesh()[patchI].size() > 0)
+            {
+                if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
+                {
+                    forAll(state.boundaryFieldRef()[patchI], faceI)
+                    {
+                        state.boundaryFieldRef()[patchI][faceI][comp] = BC;
+                    }
+                }
+                else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
+                         || state.boundaryFieldRef()[patchI].type() == "outletInlet")
+                {
+                    mixedFvPatchField<vector>& inletOutletPatch =
+                        refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
+                    vector val = inletOutletPatch.refValue()[0];
+                    val[comp] = BC;
+                    inletOutletPatch.refValue() = val;
+                }
+            }
+        }
+        else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
+        {
+            volScalarField& state(const_cast<volScalarField&>(
+                meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
+            // for decomposed domain, don't set BC if the patch is empty
+            if (meshPtr_->boundaryMesh()[patchI].size() > 0)
+            {
+                if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
+                {
+                    forAll(state.boundaryFieldRef()[patchI], faceI)
+                    {
+                        state.boundaryFieldRef()[patchI][faceI] = BC;
+                    }
+                }
+                else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
+                         || state.boundaryFieldRef()[patchI].type() == "outletInlet")
+                {
+                    mixedFvPatchField<scalar>& inletOutletPatch =
+                        refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
+                    inletOutletPatch.refValue() = BC;
+                }
+            }
+        }
+    }
+}
+
+void DASolver::getBCFromOFVars(
+    const dictionary& dvSubDict,
+    scalar& BC)
+{
+    /*
+    Description:
+        Get the boundary condition (BC) value from OF variables
+    
+    Input:
+
+        dvSubDict: the design variable subDict that contains the inform for setting BCs
+
+    Output:
+        BC: the BC value
+    */
+
+    // get info from dvSubDict. This needs to be defined in the pyDAFoam
+    // name of the variable for changing the boundary condition
+    word varName = dvSubDict.getWord("variable");
+    // name of the boundary patch
+    wordList patches;
+    dvSubDict.readEntry<wordList>("patches", patches);
+    // the component of a vector variable, ignore when it is a scalar
+    label comp = dvSubDict.getLabel("comp");
+
+    // Now get the BC value
+    forAll(patches, idxI)
+    {
+        word patchName = patches[idxI];
+        label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
+        if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
+        {
+            volVectorField& state(const_cast<volVectorField&>(
+                meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
+            // for decomposed domain, don't set BC if the patch is empty
+            if (meshPtr_->boundaryMesh()[patchI].size() > 0)
+            {
+                if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
+                {
+                    BC = state.boundaryFieldRef()[patchI][0][comp];
+                }
+                else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
+                         || state.boundaryFieldRef()[patchI].type() == "outletInlet")
+                {
+                    mixedFvPatchField<vector>& inletOutletPatch =
+                        refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
+                    BC = inletOutletPatch.refValue()[0][comp];
+                }
+            }
+        }
+        else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
+        {
+            volScalarField& state(const_cast<volScalarField&>(
+                meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
+            // for decomposed domain, don't set BC if the patch is empty
+            if (meshPtr_->boundaryMesh()[patchI].size() > 0)
+            {
+                if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
+                {
+                    BC = state.boundaryFieldRef()[patchI][0];
+                }
+                else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
+                         || state.boundaryFieldRef()[patchI].type() == "outletInlet")
+                {
+                    mixedFvPatchField<scalar>& inletOutletPatch =
+                        refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
+                    BC = inletOutletPatch.refValue()[0];
+                }
+            }
+        }
+    }
+}
+
 void DASolver::calcdRdBCTPsiAD(
     const Vec xvVec,
     const Vec wVec,
@@ -3615,61 +3830,8 @@ void DASolver::calcdRdBCTPsiAD(
     else
     {
 
-        // get info from dvSubDict. This needs to be defined in the pyDAFoam
-        // name of the variable for changing the boundary condition
-        word varName = dvSubDict.getWord("variable");
-        // name of the boundary patch
-        wordList patches;
-        dvSubDict.readEntry<wordList>("patches", patches);
-        // the component of a vector variable, ignore when it is a scalar
-        label comp = dvSubDict.getLabel("comp");
-
-        // Now get the BC value
-        forAll(patches, idxI)
-        {
-            word patchName = patches[idxI];
-            label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
-            if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
-            {
-                volVectorField& state(const_cast<volVectorField&>(
-                    meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
-                // for decomposed domain, don't set BC if the patch is empty
-                if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                {
-                    if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                    {
-                        BC = state.boundaryFieldRef()[patchI][0][comp];
-                    }
-                    else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                             || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                    {
-                        mixedFvPatchField<vector>& inletOutletPatch =
-                            refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
-                        BC = inletOutletPatch.refValue()[0][comp];
-                    }
-                }
-            }
-            else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
-            {
-                volScalarField& state(const_cast<volScalarField&>(
-                    meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
-                // for decomposed domain, don't set BC if the patch is empty
-                if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                {
-                    if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                    {
-                        BC = state.boundaryFieldRef()[patchI][0];
-                    }
-                    else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                             || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                    {
-                        mixedFvPatchField<scalar>& inletOutletPatch =
-                            refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
-                        BC = inletOutletPatch.refValue()[0];
-                    }
-                }
-            }
-        }
+        // get the BC value from the OF vars
+        this->getBCFromOFVars(dvSubDict, BC);
         // need to reduce the BC value across all processors, this is because some of
         // the processors might not own the prescribed patches so their BC value will be still -1e16, but
         // when calling the following reduce function, they will get the correct BC from other processors
@@ -3680,59 +3842,7 @@ void DASolver::calcdRdBCTPsiAD(
         // register BC as the input
         this->globalADTape_.registerInput(BC);
         // ******* now set BC ******
-        forAll(patches, idxI)
-        {
-            word patchName = patches[idxI];
-            label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
-            if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
-            {
-                volVectorField& state(const_cast<volVectorField&>(
-                    meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
-                // for decomposed domain, don't set BC if the patch is empty
-                if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                {
-                    if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                    {
-                        forAll(state.boundaryFieldRef()[patchI], faceI)
-                        {
-                            state.boundaryFieldRef()[patchI][faceI][comp] = BC;
-                        }
-                    }
-                    else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                             || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                    {
-                        mixedFvPatchField<vector>& inletOutletPatch =
-                            refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
-                        vector val = inletOutletPatch.refValue()[0];
-                        val[comp] = BC;
-                        inletOutletPatch.refValue() = val;
-                    }
-                }
-            }
-            else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
-            {
-                volScalarField& state(const_cast<volScalarField&>(
-                    meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
-                // for decomposed domain, don't set BC if the patch is empty
-                if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                {
-                    if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                    {
-                        forAll(state.boundaryFieldRef()[patchI], faceI)
-                        {
-                            state.boundaryFieldRef()[patchI][faceI] = BC;
-                        }
-                    }
-                    else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                             || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                    {
-                        mixedFvPatchField<scalar>& inletOutletPatch =
-                            refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
-                        inletOutletPatch.refValue() = BC;
-                    }
-                }
-            }
-        }
+        this->setBCToOFVars(dvSubDict, BC);
     }
     // ******* now set BC done******
     // compute residuals
@@ -3755,6 +3865,34 @@ void DASolver::calcdRdBCTPsiAD(
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
 
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    this->globalADTape_.deactivateValue(BC);
+    if (designVarName == "MRF")
+    {
+        const IOMRFZoneListDF& MRF = meshPtr_->thisDb().lookupObject<IOMRFZoneListDF>("MRFProperties");
+        scalar& omega = const_cast<scalar&>(MRF.getOmegaRef());
+        omega = BC;
+    }
+    else if (designVarName == "fvSource")
+    {
+        volVectorField& fvSource = const_cast<volVectorField&>(
+            meshPtr_->thisDb().lookupObject<volVectorField>("fvSource"));
+        label comp = dvSubDict.getLabel("comp");
+        forAll(fvSource, cellI)
+        {
+            fvSource[cellI][comp] = BC;
+        }
+    }
+    else
+    {
+        this->setBCToOFVars(dvSubDict, BC);
+    }
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
+
     wordList writeJacobians;
     daOptionPtr_->getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
     if (writeJacobians.found("dRdBCTPsi") || writeJacobians.found("all"))
@@ -3826,61 +3964,8 @@ void DASolver::calcdFdBCAD(
     }
     else
     {
-        // get info from dvSubDict. This needs to be defined in the pyDAFoam
-        // name of the variable for changing the boundary condition
-        word varName = dvSubDict.getWord("variable");
-        // name of the boundary patch
-        wordList patches;
-        dvSubDict.readEntry<wordList>("patches", patches);
-        // the compoent of a vector variable, ignore when it is a scalar
-        label comp = dvSubDict.getLabel("comp");
-
-        // Now get the BC value
-        forAll(patches, idxI)
-        {
-            word patchName = patches[idxI];
-            label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
-            if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
-            {
-                volVectorField& state(const_cast<volVectorField&>(
-                    meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
-                // for decomposed domain, don't set BC if the patch is empty
-                if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                {
-                    if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                    {
-                        BC = state.boundaryFieldRef()[patchI][0][comp];
-                    }
-                    else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                             || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                    {
-                        mixedFvPatchField<vector>& inletOutletPatch =
-                            refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
-                        BC = inletOutletPatch.refValue()[0][comp];
-                    }
-                }
-            }
-            else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
-            {
-                volScalarField& state(const_cast<volScalarField&>(
-                    meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
-                // for decomposed domain, don't set BC if the patch is empty
-                if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                {
-                    if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                    {
-                        BC = state.boundaryFieldRef()[patchI][0];
-                    }
-                    else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                             || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                    {
-                        mixedFvPatchField<scalar>& inletOutletPatch =
-                            refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
-                        BC = inletOutletPatch.refValue()[0];
-                    }
-                }
-            }
-        }
+        // get the BC value from the OF vars
+        this->getBCFromOFVars(dvSubDict, BC);
         // need to reduce the BC value across all processors, this is because some of
         // the processors might not own the prescribed patches so their BC value will be still -1e16, but
         // when calling the following reduce function, they will get the correct BC from other processors
@@ -3936,67 +4021,8 @@ void DASolver::calcdFdBCAD(
         }
         else
         {
-            // get info from dvSubDict. This needs to be defined in the pyDAFoam
-            // name of the variable for changing the boundary condition
-            word varName = dvSubDict.getWord("variable");
-            // name of the boundary patch
-            wordList patches;
-            dvSubDict.readEntry<wordList>("patches", patches);
-            // the compoent of a vector variable, ignore when it is a scalar
-            label comp = dvSubDict.getLabel("comp");
-            forAll(patches, idxI)
-            {
-                word patchName = patches[idxI];
-                label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
-                if (meshPtr_->thisDb().foundObject<volVectorField>(varName))
-                {
-                    volVectorField& state(const_cast<volVectorField&>(
-                        meshPtr_->thisDb().lookupObject<volVectorField>(varName)));
-                    // for decomposed domain, don't set BC if the patch is empty
-                    if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                    {
-                        if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                        {
-                            forAll(state.boundaryFieldRef()[patchI], faceI)
-                            {
-                                state.boundaryFieldRef()[patchI][faceI][comp] = BC;
-                            }
-                        }
-                        else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                                 || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                        {
-                            mixedFvPatchField<vector>& inletOutletPatch =
-                                refCast<mixedFvPatchField<vector>>(state.boundaryFieldRef()[patchI]);
-                            vector val = inletOutletPatch.refValue()[0];
-                            val[comp] = BC;
-                            inletOutletPatch.refValue() = val;
-                        }
-                    }
-                }
-                else if (meshPtr_->thisDb().foundObject<volScalarField>(varName))
-                {
-                    volScalarField& state(const_cast<volScalarField&>(
-                        meshPtr_->thisDb().lookupObject<volScalarField>(varName)));
-                    // for decomposed domain, don't set BC if the patch is empty
-                    if (meshPtr_->boundaryMesh()[patchI].size() > 0)
-                    {
-                        if (state.boundaryFieldRef()[patchI].type() == "fixedValue")
-                        {
-                            forAll(state.boundaryFieldRef()[patchI], faceI)
-                            {
-                                state.boundaryFieldRef()[patchI][faceI] = BC;
-                            }
-                        }
-                        else if (state.boundaryFieldRef()[patchI].type() == "inletOutlet"
-                                 || state.boundaryFieldRef()[patchI].type() == "outletInlet")
-                        {
-                            mixedFvPatchField<scalar>& inletOutletPatch =
-                                refCast<mixedFvPatchField<scalar>>(state.boundaryFieldRef()[patchI]);
-                            inletOutletPatch.refValue() = BC;
-                        }
-                    }
-                }
-            }
+            // ******* now set BC ******
+            this->setBCToOFVars(dvSubDict, BC);
         }
         // update all intermediate variables and boundary conditions
         this->updateStateBoundaryConditions();
@@ -4025,22 +4051,52 @@ void DASolver::calcdFdBCAD(
         VecAssemblyBegin(dFdBCPart);
         VecAssemblyEnd(dFdBCPart);
 
-        // need to clear adjoint and tape after the computation is done!
-        this->globalADTape_.clearAdjoints();
-        this->globalADTape_.reset();
-
         // we need to add dFdBCPart to dFdBC because we want to sum
         // all dFdBCPart for all parts of this objFuncName.
         VecAXPY(dFdBC, 1.0, dFdBCPart);
 
-        if (daOptionPtr_->getOption<label>("debug"))
+        VecDestroy(&dFdBCPart);
+
+        // need to clear adjoint and tape after the computation is done!
+        this->globalADTape_.clearAdjoints();
+        this->globalADTape_.reset();
+
+        // **********************************************************************************************
+        // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+        // **********************************************************************************************
+
+        this->globalADTape_.deactivateValue(BC);
+        if (designVarName == "MRF")
+        {
+            const IOMRFZoneListDF& MRF = meshPtr_->thisDb().lookupObject<IOMRFZoneListDF>("MRFProperties");
+            scalar& omega = const_cast<scalar&>(MRF.getOmegaRef());
+            omega = BC;
+        }
+        else if (designVarName == "fvSource")
+        {
+            volVectorField& fvSource = const_cast<volVectorField&>(
+                meshPtr_->thisDb().lookupObject<volVectorField>("fvSource"));
+
+            label comp = dvSubDict.getLabel("comp");
+
+            forAll(fvSource, cellI)
+            {
+                fvSource[cellI][comp] = BC;
+            }
+        }
+        else
+        {
+            this->setBCToOFVars(dvSubDict, BC);
+        }
+        this->updateStateBoundaryConditions();
+        fRef = daObjFunc->getObjFuncValue();
+
+        if (daOptionPtr_->getOption<label>("debug"))
         {
             Info << "In calcdFdBCAD" << endl;
             this->calcPrimalResidualStatistics("print");
             Info << objFuncName << ": " << fRef << endl;
         }
-
-        VecDestroy(&dFdBCPart);
     }
 
     wordList writeJacobians;
@@ -4334,6 +4390,48 @@ void DASolver::calcdRdAOATPsiAD(
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
 
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    this->globalADTape_.deactivateValue(aoa);
+
+    forAll(patches, idxI)
+    {
+        word patchName = patches[idxI];
+        label patchI = meshPtr_->boundaryMesh().findPatchID(patchName);
+
+        if (meshPtr_->boundaryMesh()[patchI].size() > 0)
+        {
+            scalar UMag = sqrt(Ux0 * Ux0 + Uy0 * Uy0);
+            scalar UxNew = UMag * cos(aoa);
+            scalar UyNew = UMag * sin(aoa);
+
+            if (U.boundaryField()[patchI].type() == "fixedValue")
+            {
+                forAll(U.boundaryField()[patchI], faceI)
+                {
+                    U.boundaryFieldRef()[patchI][faceI][flowAxisIndex] = UxNew;
+                    U.boundaryFieldRef()[patchI][faceI][normalAxisIndex] = UyNew;
+                }
+            }
+            else if (U.boundaryField()[patchI].type() == "inletOutlet")
+            {
+                mixedFvPatchField<vector>& inletOutletPatch =
+                    refCast<mixedFvPatchField<vector>>(U.boundaryFieldRef()[patchI]);
+
+                forAll(U.boundaryField()[patchI], faceI)
+                {
+                    inletOutletPatch.refValue()[faceI][flowAxisIndex] = UxNew;
+                    inletOutletPatch.refValue()[faceI][normalAxisIndex] = UyNew;
+                }
+            }
+        }
+    }
+    // compute residuals
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
+
     wordList writeJacobians;
     daOptionPtr_->getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
     if (writeJacobians.found("dRdAOATPsi") || writeJacobians.found("all"))
@@ -4837,13 +4935,22 @@ void DASolver::calcdFdFieldAD(
         VecZeroEntries(dFdFieldPart);
         this->assignFieldGradient2Vec(fieldName, fieldType, dFdFieldPart);
 
+        // we need to add dFdFieldPart to dFdField because we want to sum
+        // all dFdFieldPart for all parts of this objFuncName.
+        VecAXPY(dFdField, 1.0, dFdFieldPart);
+
         // need to clear adjoint and tape after the computation is done!
         this->globalADTape_.clearAdjoints();
         this->globalADTape_.reset();
 
-        // we need to add dFdFieldPart to dFdField because we want to sum
-        // all dFdFieldPart for all parts of this objFuncName.
-        VecAXPY(dFdField, 1.0, dFdFieldPart);
+        // **********************************************************************************************
+        // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+        // **********************************************************************************************
+
+        this->deactivateFieldVariableInput4AD(fieldName, fieldType);
+        this->updateBoundaryConditions(fieldName, fieldType);
+        this->updateStateBoundaryConditions();
+        fRef = daObjFunc->getObjFuncValue();
 
         if (daOptionPtr_->getOption<label>("debug"))
         {
@@ -5253,13 +5360,21 @@ void DASolver::calcdFdWAD(
         VecZeroEntries(dFdWPart);
         this->assignStateGradient2Vec(dFdWPart);
 
+        // we need to add dFdWPart to dFdW because we want to sum
+        // all dFdWPart for all parts of this objFuncName.
+        VecAXPY(dFdW, 1.0, dFdWPart);
+
         // need to clear adjoint and tape after the computation is done!
         this->globalADTape_.clearAdjoints();
         this->globalADTape_.reset();
 
-        // we need to add dFdWPart to dFdW because we want to sum
-        // all dFdWPart for all parts of this objFuncName.
-        VecAXPY(dFdW, 1.0, dFdWPart);
+        // **********************************************************************************************
+        // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+        // **********************************************************************************************
+
+        this->deactivateStateVariableInput4AD();
+        this->updateStateBoundaryConditions();
+        fRef = daObjFunc->getObjFuncValue();
 
         if (daOptionPtr_->getOption<label>("debug"))
         {
@@ -5392,21 +5507,37 @@ void DASolver::calcdFdXvAD(
         VecAssemblyBegin(dFdXvPart);
         VecAssemblyEnd(dFdXvPart);
 
+        // we need to add dFd*Part to dFd* because we want to sum
+        // all dFd*Part for all parts of this objFuncName.
+        VecAXPY(dFdXv, 1.0, dFdXvPart);
+
+        VecDestroy(&dFdXvPart);
+
         // need to clear adjoint and tape after the computation is done!
         this->globalADTape_.clearAdjoints();
         this->globalADTape_.reset();
 
-        // we need to add dFd*Part to dFd* because we want to sum
-        // all dFd*Part for all parts of this objFuncName.
-        VecAXPY(dFdXv, 1.0, dFdXvPart);
+        // **********************************************************************************************
+        // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+        // **********************************************************************************************
+
+        forAll(meshPoints, i)
+        {
+            for (label j = 0; j < 3; j++)
+            {
+                this->globalADTape_.deactivateValue(meshPoints[i][j]);
+            }
+        }
+        meshPtr_->movePoints(meshPoints);
+        meshPtr_->moving(false);
+        this->updateStateBoundaryConditions();
+        fRef = daObjFunc->getObjFuncValue();
 
         if (daOptionPtr_->getOption<label>("debug"))
         {
             this->calcPrimalResidualStatistics("print");
             Info << objFuncName << ": " << fRef << endl;
         }
-
-        VecDestroy(&dFdXvPart);
     }
 
     wordList writeJacobians;
@@ -5506,33 +5637,26 @@ void DASolver::calcdRdThermalTPsiAD(
         product[i] = thermalArray[i].getGradient();
     }
 
-    // ******************************* NOTE! ******************************
-    // This step is tricky. We need to clean up the AD seeds for all variables before the next AD
-    // This is usually done by theses two calls: updateOFField(wVec) and updateOFMesh(xvVec)
-    // The above two calls will assign zeros to the gradient() part of all scalar variables because
-    // the wVec and xvVec's gradient() part is zero (they are double instead of scalar).
-    // HOWEVER, the above two calls will NOT clean the seeds for field variables's boundary values
-    // This will not cause problems for most of the functions because their boundary values always
-    // have zero seeds. However, the setThermal function actually sets non-zeros seeds for boundary
-    // values, so here we need to manually clean up the BC's gradient part. NOT doing this will
-    // mess up the AD computation for the following function.
-    // The way we clean it is that we set a zero-seed thermal var and call setThermal to propagate
-    // the zero seeds to all intermediate vars.
-    // ******************************* NOTE! ******************************
+    // clean up AD
+    this->globalADTape_.clearAdjoints();
+    this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
     for (label i = 0; i < nCouplingFaces * 2; i++)
     {
-        thermalArray[i].setGradient(0.0);
+        this->globalADTape_.deactivateValue(thermalArray[i]);
     }
     this->setThermal(thermalArray);
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
 
     delete[] thermalArray;
     delete[] volCoordsArray;
     delete[] statesArray;
 
-    // clean up AD
-    this->globalADTape_.clearAdjoints();
-    this->globalADTape_.reset();
-
 #endif
 }
 
@@ -5603,6 +5727,23 @@ void DASolver::calcdRdXvTPsiAD(
 
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    forAll(meshPoints, i)
+    {
+        for (label j = 0; j < 3; j++)
+        {
+            this->globalADTape_.deactivateValue(meshPoints[i][j]);
+        }
+    }
+    meshPtr_->movePoints(meshPoints);
+    meshPtr_->moving(false);
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
+
 #endif
 }
 
@@ -5682,6 +5823,23 @@ void DASolver::calcdForcedXvAD(
 
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    forAll(meshPoints, i)
+    {
+        for (label j = 0; j < 3; j++)
+        {
+            this->globalADTape_.deactivateValue(meshPoints[i][j]);
+        }
+    }
+    meshPtr_->movePoints(meshPoints);
+    meshPtr_->moving(false);
+    this->updateStateBoundaryConditions();
+    this->getForcesInternal(fX, fY, fZ, patchList);
+
 #endif
 }
 
@@ -5813,6 +5971,23 @@ void DASolver::calcdAcousticsdXvAD(
 
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    forAll(meshPoints, i)
+    {
+        for (label j = 0; j < 3; j++)
+        {
+            this->globalADTape_.deactivateValue(meshPoints[i][j]);
+        }
+    }
+    meshPtr_->movePoints(meshPoints);
+    meshPtr_->moving(false);
+    this->updateStateBoundaryConditions();
+    this->getAcousticDataInternal(x, y, z, nX, nY, nZ, a, fX, fY, fZ, patchList);
+
 #endif
 }
 
@@ -5876,6 +6051,15 @@ void DASolver::calcdRdFieldTPsiAD(
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
 
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    this->deactivateFieldVariableInput4AD(fieldName, fieldType);
+    this->updateBoundaryConditions(fieldName, fieldType);
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
+
     wordList writeJacobians;
     daOptionPtr_->getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
     if (writeJacobians.found("dRdFieldTPsi") || writeJacobians.found("all"))
@@ -6015,21 +6199,37 @@ void DASolver::calcdFdACTAD(
                 VecAssemblyBegin(dFdACTPart);
                 VecAssemblyEnd(dFdACTPart);
 
+                // we need to add dFd*Part to dFd* because we want to sum
+                // all dFd*Part for all parts of this objFuncName.
+                VecAXPY(dFdACT, 1.0, dFdACTPart);
+
+                VecDestroy(&dFdACTPart);
+
                 // need to clear adjoint and tape after the computation is done!
                 this->globalADTape_.clearAdjoints();
                 this->globalADTape_.reset();
 
-                // we need to add dFd*Part to dFd* because we want to sum
-                // all dFd*Part for all parts of this objFuncName.
-                VecAXPY(dFdACT, 1.0, dFdACTPart);
+                // **********************************************************************************************
+                // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+                // **********************************************************************************************
+
+                for (label i = 0; i < 10; i++)
+                {
+                    this->globalADTape_.deactivateValue(actDVList[i]);
+                }
+                for (label i = 0; i < 10; i++)
+                {
+                    fvSource.setActuatorDVs(diskName, i, actDVList[i]);
+                }
+                fvSource.updateFvSource();
+                this->updateStateBoundaryConditions();
+                fRef = daObjFunc->getObjFuncValue();
 
                 if (daOptionPtr_->getOption<label>("debug"))
                 {
                     this->calcPrimalResidualStatistics("print");
                     Info << objFuncName << ": " << fRef << endl;
                 }
-
-                VecDestroy(&dFdACTPart);
             }
         }
         else
@@ -6145,6 +6345,21 @@ void DASolver::calcdRdActTPsiAD(
 
             this->globalADTape_.clearAdjoints();
             this->globalADTape_.reset();
+
+            // **********************************************************************************************
+            // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+            // **********************************************************************************************
+
+            for (label i = 0; i < 10; i++)
+            {
+                this->globalADTape_.deactivateValue(actDVList[i]);
+            }
+            for (label i = 0; i < 10; i++)
+            {
+                fvSource.setActuatorDVs(diskName, i, actDVList[i]);
+            }
+            this->updateStateBoundaryConditions();
+            this->calcResiduals();
         }
         else
         {
@@ -6235,6 +6450,14 @@ void DASolver::calcdForcedWAD(
 
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    this->deactivateStateVariableInput4AD();
+    this->updateStateBoundaryConditions();
+    this->getForcesInternal(fX, fY, fZ, patchList);
 #endif
 }
 
@@ -6357,6 +6580,14 @@ void DASolver::calcdAcousticsdWAD(
 
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
+
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    this->deactivateStateVariableInput4AD();
+    this->updateStateBoundaryConditions();
+    this->getAcousticDataInternal(x, y, z, nX, nY, nZ, a, fX, fY, fZ, patchList);
 #endif
 }
 
@@ -6396,15 +6627,8 @@ void DASolver::calcdRdWTPsiAD(
         this->registerStateVariableInput4AD();
 
         // compute residuals
-        daResidualPtr_->correctBoundaryConditions();
-        daResidualPtr_->updateIntermediateVariables();
-        daModelPtr_->correctBoundaryConditions();
-        daModelPtr_->updateIntermediateVariables();
-        label isPC = 0;
-        dictionary options;
-        options.set("isPC", isPC);
-        daResidualPtr_->calcResiduals(options);
-        daModelPtr_->calcResiduals(options);
+        this->updateStateBoundaryConditions();
+        this->calcResiduals();
 
         this->registerResidualOutput4AD();
         this->globalADTape_.setPassive();
@@ -6484,6 +6708,13 @@ void DASolver::calcdRdWTPsiAD(
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
 
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+    this->deactivateStateVariableInput4AD();
+    this->updateStateBoundaryConditions();
+    this->calcResiduals();
+
     wordList writeJacobians;
     daOptionPtr_->getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
     if (writeJacobians.found("dRdWTPsi") || writeJacobians.found("all"))
@@ -6552,6 +6783,13 @@ void DASolver::calcdRdWOldTPsiAD(
     this->globalADTape_.clearAdjoints();
     this->globalADTape_.reset();
 
+    // **********************************************************************************************
+    // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+    // **********************************************************************************************
+
+    // ************** TODO *****************
+    // We need to deactivate the old time. NOT implemented yet.
+
     wordList writeJacobians;
     daOptionPtr_->getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
     if (writeJacobians.found("dRdWOldTPsi") || writeJacobians.found("all"))
@@ -6720,6 +6958,75 @@ void DASolver::registerStateVariableInput4AD(const label oldTimeLevel)
 #endif
 }
 
+void DASolver::deactivateStateVariableInput4AD()
+{
+#ifdef CODI_AD_REVERSE
+    /*
+    Description:
+        Deactivate all state variables as the input for reverse-mode AD
+    */
+
+    forAll(stateInfo_["volVectorStates"], idxI)
+    {
+        const word stateName = stateInfo_["volVectorStates"][idxI];
+        volVectorField& state = const_cast<volVectorField&>(
+            meshPtr_->thisDb().lookupObject<volVectorField>(stateName));
+
+        forAll(state, cellI)
+        {
+            for (label i = 0; i < 3; i++)
+            {
+                this->globalADTape_.deactivateValue(state[cellI][i]);
+            }
+        }
+    }
+
+    forAll(stateInfo_["volScalarStates"], idxI)
+    {
+        const word stateName = stateInfo_["volScalarStates"][idxI];
+        volScalarField& state = const_cast<volScalarField&>(
+            meshPtr_->thisDb().lookupObject<volScalarField>(stateName));
+
+        forAll(state, cellI)
+        {
+            this->globalADTape_.deactivateValue(state[cellI]);
+        }
+    }
+
+    forAll(stateInfo_["modelStates"], idxI)
+    {
+        const word stateName = stateInfo_["modelStates"][idxI];
+        volScalarField& state = const_cast<volScalarField&>(
+            meshPtr_->thisDb().lookupObject<volScalarField>(stateName));
+
+        forAll(state, cellI)
+        {
+
+            this->globalADTape_.deactivateValue(state[cellI]);
+        }
+    }
+
+    forAll(stateInfo_["surfaceScalarStates"], idxI)
+    {
+        const word stateName = stateInfo_["surfaceScalarStates"][idxI];
+        surfaceScalarField& state = const_cast<surfaceScalarField&>(
+            meshPtr_->thisDb().lookupObject<surfaceScalarField>(stateName));
+        forAll(state, faceI)
+        {
+            this->globalADTape_.deactivateValue(state[faceI]);
+        }
+        forAll(state.boundaryField(), patchI)
+        {
+            forAll(state.boundaryField()[patchI], faceI)
+            {
+                this->globalADTape_.deactivateValue(state.boundaryFieldRef()[patchI][faceI]);
+            }
+        }
+    }
+
+#endif
+}
+
 void DASolver::registerFieldVariableInput4AD(
     const word fieldName,
     const word fieldType)
@@ -6767,6 +7074,53 @@ void DASolver::registerFieldVariableInput4AD(
 #endif
 }
 
+void DASolver::deactivateFieldVariableInput4AD(
+    const word fieldName,
+    const word fieldType)
+{
+#ifdef CODI_AD_REVERSE
+    /*
+    Description:
+        Deacitvate field variables as the input for reverse-mode AD
+    
+    Input:
+        fieldName: the name of the flow field to register
+
+        fieldType: can be either scalar or vector
+    */
+
+    if (fieldType == "scalar")
+    {
+        volScalarField& state = const_cast<volScalarField&>(
+            meshPtr_->thisDb().lookupObject<volScalarField>(fieldName));
+
+        forAll(state, cellI)
+        {
+            this->globalADTape_.deactivateValue(state[cellI]);
+        }
+    }
+    else if (fieldType == "vector")
+    {
+        volVectorField& state = const_cast<volVectorField&>(
+            meshPtr_->thisDb().lookupObject<volVectorField>(fieldName));
+
+        forAll(state, cellI)
+        {
+            for (label i = 0; i < 3; i++)
+            {
+                this->globalADTape_.deactivateValue(state[cellI][i]);
+            }
+        }
+    }
+    else
+    {
+        FatalErrorIn("") << "fieldType not valid. Options: scalar or vector"
+                         << abort(FatalError);
+    }
+
+#endif
+}
+
 void DASolver::registerResidualOutput4AD()
 {
 #ifdef CODI_AD_REVERSE
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index 6aa09b8a..1ba84aaf 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -346,6 +346,16 @@ public:
         const word designVarName,
         Vec dFdBC);
 
+    /// assign the BC val to OF var
+    void setBCToOFVars(
+        const dictionary& dvSubDict,
+        const scalar& BC);
+    
+    /// get the BC val from OF var
+    void getBCFromOFVars(
+        const dictionary& dvSubDict,
+        scalar& BC);
+
     /// compute dRdBCAD
     void calcdRdBCTPsiAD(
         const Vec xvVec,
@@ -538,11 +548,19 @@ public:
     /// register all state variables as the input for reverse-mode AD
     void registerStateVariableInput4AD(const label oldTimeLevel = 0);
 
+    /// deactivate all state variables as the input for reverse-mode AD
+    void deactivateStateVariableInput4AD();
+
     /// register field variables as the input for reverse-mode AD
     void registerFieldVariableInput4AD(
         const word fieldName,
         const word fieldType);
 
+    /// deactivate the field variables for reverse-mode AD
+    void deactivateFieldVariableInput4AD(
+        const word fieldName,
+        const word fieldType);
+
     /// register all residuals as the output for reverse-mode AD
     void registerResidualOutput4AD();
 
@@ -904,7 +922,7 @@ public:
     void updateBoundaryConditions(
         const word fieldName,
         const word fieldType);
-    
+
     /// update the boundary conditions for all states and intermediate variables
     void updateStateBoundaryConditions();
 
@@ -997,7 +1015,7 @@ public:
 
     /// initialize tensorflow functions and interfaces for callback
     void initTensorFlowFuncs(
-        pyComputeInterface computeInterface, 
+        pyComputeInterface computeInterface,
         void* compute,
         pyJacVecProdInterface jacVecProdInterface,
         void* jacVecProd)
diff --git a/tests/refs/DAFoam_Test_DASimpleFoamkOmegaFieldInversionOmegaRef.txt b/tests/refs/DAFoam_Test_DASimpleFoamkOmegaFieldInversionOmegaRef.txt
index 5ec04f64..08ff14ca 100755
--- a/tests/refs/DAFoam_Test_DASimpleFoamkOmegaFieldInversionOmegaRef.txt
+++ b/tests/refs/DAFoam_Test_DASimpleFoamkOmegaFieldInversionOmegaRef.txt
@@ -6,7 +6,7 @@ Dictionary Key: FI
 Dictionary Key: alpha
 @value                    0 0.0001 1e-06
 Dictionary Key: alphaPorosity
-@value    7667.002672988612 0.0001 1e-06
+@value    7662.95674760256 0.0001 1e-06
 Dictionary Key: beta
 @value    252.4756985454485 0.0001 1e-06
 Dictionary Key: fail
diff --git a/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFieldInversionRef.txt b/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFieldInversionRef.txt
index d23bd943..89dcece3 100755
--- a/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFieldInversionRef.txt
+++ b/tests/refs/DAFoam_Test_DASimpleFoamkOmegaSSTFieldInversionRef.txt
@@ -10,6 +10,6 @@ Dictionary Key: alphaPorosity
 Dictionary Key: beta
 @value    63.62415240726151 0.0001 1e-06
 Dictionary Key: fvSource
-@value    1806.201838099202 0.0001 1e-06
+@value    1810.517315729984 0.0001 1e-06
 Dictionary Key: fail
 @value                    0 0.0001 1e-06
diff --git a/tests/refs/DAFoam_Test_MphysAeroFieldRef.txt b/tests/refs/DAFoam_Test_MphysAeroFieldRef.txt
index 4ffbe37c..4457daa7 100755
--- a/tests/refs/DAFoam_Test_MphysAeroFieldRef.txt
+++ b/tests/refs/DAFoam_Test_MphysAeroFieldRef.txt
@@ -1,26 +1,9 @@
 Dictionary Key: CL
+@value   0.2992762848982257 1e-06 1e-08
+Dictionary Key: CL
 Dictionary Key: alpha
-@value   0.6405997854447981 0.0001 1e-06
+@value    0.435012393483863 0.0001 1e-06
 Dictionary Key: aoa
-@value 0.000959052564689386 0.0001 1e-06
-Dictionary Key: shape
-@value   -1.359174005645564 0.0001 1e-06
-@value  -0.6692178532758711 0.0001 1e-06
-@value  -0.2149441233470428 0.0001 1e-06
-@value -0.01334827639681135 0.0001 1e-06
-@value    1.666569177875879 0.0001 1e-06
-@value     1.11477785597916 0.0001 1e-06
-@value    1.624331137829298 0.0001 1e-06
-@value   0.4017427848906134 0.0001 1e-06
-@value   0.2919486313792322 0.0001 1e-06
-@value  0.04132164500265292 0.0001 1e-06
-@value   0.9052707541190427 0.0001 1e-06
-@value   0.4881671535089137 0.0001 1e-06
-@value    1.280749801412769 0.0001 1e-06
-@value    1.203706208898332 0.0001 1e-06
-@value    1.276188815028655 0.0001 1e-06
-@value    1.160544344322793 0.0001 1e-06
-@value   -2.227350578087959 0.0001 1e-06
-@value   -2.199728454741364 0.0001 1e-06
-@value   -2.400190424831512 0.0001 1e-06
-@value   -2.371364593921126 0.0001 1e-06
\ No newline at end of file
+@value  0.09767145516644697 0.0001 1e-06
+Dictionary Key: twist
+@value  0.09790920083964022 0.0001 1e-06
\ No newline at end of file
diff --git a/tests/refs/DAFoam_Test_MphysAeroOptRef.txt b/tests/refs/DAFoam_Test_MphysAeroOptRef.txt
index ec5777f3..45adbefe 100755
--- a/tests/refs/DAFoam_Test_MphysAeroOptRef.txt
+++ b/tests/refs/DAFoam_Test_MphysAeroOptRef.txt
@@ -1,4 +1,4 @@
 Dictionary Key: cruise.aero_post.CD
-@value  0.03919004721810314 1e-06 1e-08
+@value  0.03918787297431089 1e-06 1e-08
 Dictionary Key: cruise.aero_post.CL
-@value   0.2999950298843596 1e-06 1e-08
\ No newline at end of file
+@value   0.2999891179061076 1e-06 1e-08
\ No newline at end of file
diff --git a/tests/refs/DAFoam_Test_MphysAeroPropCoupledRef.txt b/tests/refs/DAFoam_Test_MphysAeroPropCoupledRef.txt
index 272f05ff..c76c69f0 100755
--- a/tests/refs/DAFoam_Test_MphysAeroPropCoupledRef.txt
+++ b/tests/refs/DAFoam_Test_MphysAeroPropCoupledRef.txt
@@ -8,85 +8,85 @@ Dictionary Key: thrust
 @value   0.3611900799886527 1e-06 1e-10
 Dictionary Key: force_balance
 Dictionary Key: prop_shape
-@value    40.56645432033733 0.0001 1e-06
-@value   -558.1380598187881 0.0001 1e-06
-@value    467.7734184945811 0.0001 1e-06
-@value   -158.9914676337573 0.0001 1e-06
-@value   -85.02618895954103 0.0001 1e-06
-@value   -627.9557803116626 0.0001 1e-06
-@value    140.3194131387392 0.0001 1e-06
-@value  -0.1713762746667271 0.0001 1e-06
-@value    -508.127288548224 0.0001 1e-06
-@value   -540.6253032464376 0.0001 1e-06
-@value   -104.7456081289921 0.0001 1e-06
-@value   -161.1842756558842 0.0001 1e-06
-@value    60.59671558598635 0.0001 1e-06
-@value    116.2388390235131 0.0001 1e-06
-@value    354.5260914116756 0.0001 1e-06
-@value    336.3699829029644 0.0001 1e-06
-@value   -347.2898519131134 0.0001 1e-06
-@value   -130.0392459478966 0.0001 1e-06
-@value    525.1168387426317 0.0001 1e-06
-@value    298.6611927164653 0.0001 1e-06
-@value   -108.9449106018517 0.0001 1e-06
-@value   -62.29976349182109 0.0001 1e-06
-@value    341.6203318732158 0.0001 1e-06
-@value    263.7148977578634 0.0001 1e-06
+@value    40.57011447892199 0.0001 1e-06
+@value   -558.1394156909686 0.0001 1e-06
+@value    467.7771921270105 0.0001 1e-06
+@value   -158.9944668735331 0.0001 1e-06
+@value   -85.02737008055752 0.0001 1e-06
+@value   -627.9548778289284 0.0001 1e-06
+@value     140.318474446663 0.0001 1e-06
+@value  -0.1694117259289669 0.0001 1e-06
+@value   -508.1280147631104 0.0001 1e-06
+@value   -540.6248952272629 0.0001 1e-06
+@value     -104.74310426356 0.0001 1e-06
+@value    -161.184229377824 0.0001 1e-06
+@value    60.59491524505091 0.0001 1e-06
+@value    116.2393332241295 0.0001 1e-06
+@value    354.5286978848217 0.0001 1e-06
+@value    336.3699874842462 0.0001 1e-06
+@value   -347.2930367286646 0.0001 1e-06
+@value   -130.0372517970196 0.0001 1e-06
+@value    525.1221596957737 0.0001 1e-06
+@value    298.6593247134094 0.0001 1e-06
+@value   -108.9467995921807 0.0001 1e-06
+@value      -62.29917598906 0.0001 1e-06
+@value    341.6238879060646 0.0001 1e-06
+@value    263.7148030220755 0.0001 1e-06
 Dictionary Key: wing_twist
-@value  -0.5910416434157874 0.0001 1e-06
+@value  -0.6212205045129653 0.0001 1e-06
 Dictionary Key: lift
 Dictionary Key: prop_shape
-@value    7.494632028196025 0.0001 1e-06
-@value   -2.691981732712854 0.0001 1e-06
-@value    6.893016044036441 0.0001 1e-06
-@value   -5.991053822821106 0.0001 1e-06
-@value   -5.739115722020717 0.0001 1e-06
-@value     5.47681875189604 0.0001 1e-06
-@value   -4.164352933049335 0.0001 1e-06
-@value    13.88196713361107 0.0001 1e-06
-@value   -2.235091169101973 0.0001 1e-06
-@value   -3.081434681823781 0.0001 1e-06
-@value    6.324731236348367 0.0001 1e-06
-@value    9.840864367457241 0.0001 1e-06
-@value    -7.31659760544235 0.0001 1e-06
-@value   -10.80796786398239 0.0001 1e-06
-@value    5.802828854432147 0.0001 1e-06
-@value    8.951835560847094 0.0001 1e-06
-@value   -6.078515425622227 0.0001 1e-06
-@value    14.28172057520801 0.0001 1e-06
-@value    10.48958602037796 0.0001 1e-06
-@value   -9.564279220170844 0.0001 1e-06
-@value   -6.627024939462579 0.0001 1e-06
-@value   -9.842049987973438 0.0001 1e-06
-@value     9.34366016427688 0.0001 1e-06
-@value    14.36183778196542 0.0001 1e-06
+@value     7.37755063863561 0.0001 1e-06
+@value   -2.650727700451313 0.0001 1e-06
+@value    6.790617938688313 0.0001 1e-06
+@value   -5.891499618985533 0.0001 1e-06
+@value   -5.600939745650626 0.0001 1e-06
+@value    5.342123399123075 0.0001 1e-06
+@value   -4.066653677004293 0.0001 1e-06
+@value     13.5383549122466 0.0001 1e-06
+@value   -2.179837392617058 0.0001 1e-06
+@value   -2.993058240139202 0.0001 1e-06
+@value    6.178084688158777 0.0001 1e-06
+@value    9.576108764718612 0.0001 1e-06
+@value   -7.134496985715375 0.0001 1e-06
+@value   -10.51253550147928 0.0001 1e-06
+@value    5.669359902494561 0.0001 1e-06
+@value    8.708170216943536 0.0001 1e-06
+@value   -5.982348335297536 0.0001 1e-06
+@value    13.92869912246976 0.0001 1e-06
+@value     10.3269883748625 0.0001 1e-06
+@value   -9.332879662044292 0.0001 1e-06
+@value   -6.462275830980414 0.0001 1e-06
+@value   -9.569477213299155 0.0001 1e-06
+@value      9.1247009816325 0.0001 1e-06
+@value    13.97371605464054 0.0001 1e-06
 Dictionary Key: wing_twist
-@value     67.6593267198632 0.0001 1e-06
+@value      66.712745267288 0.0001 1e-06
 Dictionary Key: power
 Dictionary Key: prop_shape
-@value    21.64495476036015 0.0001 1e-06
-@value    9.726766089966382 0.0001 1e-06
-@value    -16.4582865984792 0.0001 1e-06
-@value   -66.37785975292292 0.0001 1e-06
-@value    11.59301925605034 0.0001 1e-06
-@value    29.13037119474476 0.0001 1e-06
-@value   -29.55344784858206 0.0001 1e-06
-@value   -4.154484479118967 0.0001 1e-06
-@value    22.22860938531955 0.0001 1e-06
-@value    25.97576332452039 0.0001 1e-06
-@value     19.6348794304481 0.0001 1e-06
-@value    19.12141987741848 0.0001 1e-06
-@value   -50.03884470681318 0.0001 1e-06
-@value   -32.61146584587635 0.0001 1e-06
-@value   -32.06312274498524 0.0001 1e-06
-@value    -38.1047911573768 0.0001 1e-06
-@value   0.5725780768975435 0.0001 1e-06
-@value    23.22195308385148 0.0001 1e-06
-@value    33.42823149513347 0.0001 1e-06
-@value    10.06960227413743 0.0001 1e-06
-@value    13.88047090531864 0.0001 1e-06
-@value    16.78764894295167 0.0001 1e-06
-@value    27.70306167537909 0.0001 1e-06
-@value    24.33006926407132 0.0001 1e-06
+@value    21.64495476035998 0.0001 1e-06
+@value    9.726766089966594 0.0001 1e-06
+@value   -16.45828659847884 0.0001 1e-06
+@value   -66.37785975292272 0.0001 1e-06
+@value    11.59301925605003 0.0001 1e-06
+@value    29.13037119474438 0.0001 1e-06
+@value   -29.55344784858185 0.0001 1e-06
+@value   -4.154484479118828 0.0001 1e-06
+@value    22.22860938531966 0.0001 1e-06
+@value    25.97576332452035 0.0001 1e-06
+@value    19.63487943044795 0.0001 1e-06
+@value    19.12141987741839 0.0001 1e-06
+@value   -50.03884470681285 0.0001 1e-06
+@value   -32.61146584587608 0.0001 1e-06
+@value   -32.06312274498497 0.0001 1e-06
+@value   -38.10479115737638 0.0001 1e-06
+@value   0.5725780768974928 0.0001 1e-06
+@value     23.2219530838513 0.0001 1e-06
+@value     33.4282314951332 0.0001 1e-06
+@value    10.06960227413714 0.0001 1e-06
+@value    13.88047090531867 0.0001 1e-06
+@value    16.78764894295168 0.0001 1e-06
+@value    27.70306167537875 0.0001 1e-06
+@value    24.33006926407115 0.0001 1e-06
 Dictionary Key: wing_twist
 @value                   -0 0.0001 1e-06
\ No newline at end of file
diff --git a/tests/refs/DAFoam_Test_MphysAeroPropRef.txt b/tests/refs/DAFoam_Test_MphysAeroPropRef.txt
index 2c2c16c1..ce4770fb 100755
--- a/tests/refs/DAFoam_Test_MphysAeroPropRef.txt
+++ b/tests/refs/DAFoam_Test_MphysAeroPropRef.txt
@@ -4,101 +4,101 @@ Dictionary Key: CL
 @value   -80.60966769315618 1e-06 1e-10
 Dictionary Key: CL
 Dictionary Key: prop1_axial_force
-@value 8.753887553475368e-05 0.0001 1e-06
-@value 4.296745116705247e-05 0.0001 1e-06
-@value 3.797080953198209e-05 0.0001 1e-06
-@value 2.843684634830516e-05 0.0001 1e-06
-@value 8.706391713808491e-06 0.0001 1e-06
-@value 1.012100363320951e-05 0.0001 1e-06
-@value 1.616578589963528e-05 0.0001 1e-06
-@value -1.142944108417678e-05 0.0001 1e-06
-@value -1.624692742501505e-05 0.0001 1e-06
-@value -0.0001053385574004885 0.0001 1e-06
+@value 8.706649357462623e-05 0.0001 1e-06
+@value 4.284329277164987e-05 0.0001 1e-06
+@value 3.791145792304343e-05 0.0001 1e-06
+@value 2.841656557985818e-05 0.0001 1e-06
+@value 8.750321344966871e-06 0.0001 1e-06
+@value 1.005890227350335e-05 0.0001 1e-06
+@value 1.606953345911087e-05 0.0001 1e-06
+@value -1.139506687603346e-05 0.0001 1e-06
+@value -1.615341411011944e-05 0.0001 1e-06
+@value -0.0001051186517921878 0.0001 1e-06
 Dictionary Key: prop1_integral_force
-@value 3.071965816062354e-05 0.0001 1e-06
-@value -6.636078874801797e-05 0.0001 1e-06
+@value 3.078505250423354e-05 0.0001 1e-06
+@value -6.64152399588178e-05 0.0001 1e-06
 Dictionary Key: prop1_prop_center
-@value 0.002262169443610162 0.0001 1e-06
-@value -0.0005734337364295363 0.0001 1e-06
-@value 0.000549839927615936 0.0001 1e-06
+@value 0.002262456380060131 0.0001 1e-06
+@value -0.0005725609070018943 0.0001 1e-06
+@value 0.0005508250136281442 0.0001 1e-06
 Dictionary Key: prop1_radial_center
-@value -4.673906013155489e-05 0.0001 1e-06
-@value -3.228475071137402e-05 0.0001 1e-06
-@value -1.133131854193097e-05 0.0001 1e-06
-@value 6.558823621290013e-06 0.0001 1e-06
-@value -7.862648656963834e-06 0.0001 1e-06
-@value 1.336802485581374e-05 0.0001 1e-06
-@value 2.465533238564196e-06 0.0001 1e-06
-@value 6.32653888921168e-06 0.0001 1e-06
-@value -5.200913730081363e-05 0.0001 1e-06
-@value -2.056635860143015e-05 0.0001 1e-06
+@value -4.643512838473826e-05 0.0001 1e-06
+@value -3.212302361662417e-05 0.0001 1e-06
+@value -1.11799283033605e-05 0.0001 1e-06
+@value 6.609980691156762e-06 0.0001 1e-06
+@value -7.852457740836726e-06 0.0001 1e-06
+@value 1.342589495613161e-05 0.0001 1e-06
+@value 2.449306978788844e-06 0.0001 1e-06
+@value 6.379422866753332e-06 0.0001 1e-06
+@value -5.183007541162234e-05 0.0001 1e-06
+@value -2.04580054843637e-05 0.0001 1e-06
 Dictionary Key: prop1_tangential_force
-@value -2.783537821900966e-05 0.0001 1e-06
-@value -1.727489560307937e-05 0.0001 1e-06
-@value -3.060727120155457e-05 0.0001 1e-06
-@value -3.439969536084378e-05 0.0001 1e-06
-@value 3.207596472006431e-07 0.0001 1e-06
-@value -2.804555129074096e-05 0.0001 1e-06
-@value -1.671663354771918e-05 0.0001 1e-06
-@value 8.28830364480922e-06 0.0001 1e-06
-@value -6.119725484897973e-06 0.0001 1e-06
-@value 0.0001322946710949994 0.0001 1e-06
+@value -2.774389378112688e-05 0.0001 1e-06
+@value -1.727658705445228e-05 0.0001 1e-06
+@value -3.065762494568778e-05 0.0001 1e-06
+@value -3.440949228601185e-05 0.0001 1e-06
+@value 2.554586873888304e-07 0.0001 1e-06
+@value -2.804823760239071e-05 0.0001 1e-06
+@value -1.661138587518807e-05 0.0001 1e-06
+@value 8.240985016240634e-06 0.0001 1e-06
+@value -6.058465902711757e-06 0.0001 1e-06
+@value 0.0001322211841405256 0.0001 1e-06
 Dictionary Key: prop2_axial_force
-@value 3.769467716716783e-05 0.0001 1e-06
-@value 0.0001051407693863429 0.0001 1e-06
-@value 0.0001624682213236707 0.0001 1e-06
-@value 0.0002006653437954606 0.0001 1e-06
-@value 0.0002216842576897312 0.0001 1e-06
-@value 0.0002633393102510438 0.0001 1e-06
-@value 7.290682709672559e-05 0.0001 1e-06
-@value -8.0484348477574e-05 0.0001 1e-06
-@value 0.0004665144513131048 0.0001 1e-06
-@value -0.001659864118688367 0.0001 1e-06
+@value 3.777172000862125e-05 0.0001 1e-06
+@value 0.0001052384601978654 0.0001 1e-06
+@value 0.0001626363983463342 0.0001 1e-06
+@value 0.0002007220814366768 0.0001 1e-06
+@value 0.0002214900773800584 0.0001 1e-06
+@value 0.0002630201764147135 0.0001 1e-06
+@value 7.272487136889523e-05 0.0001 1e-06
+@value -8.052064798912409e-05 0.0001 1e-06
+@value 0.0004660973215747894 0.0001 1e-06
+@value -0.001658518503947827 0.0001 1e-06
 Dictionary Key: prop2_integral_force
-@value 0.0001507599834106935 0.0001 1e-06
-@value -6.565431367562223e-07 0.0001 1e-06
+@value 0.0001507458586159512 0.0001 1e-06
+@value -6.536940084709029e-07 0.0001 1e-06
 Dictionary Key: prop2_prop_center
-@value 0.008035349668816954 0.0001 1e-06
-@value -0.003092132713889739 0.0001 1e-06
-@value -0.003272949858479228 0.0001 1e-06
+@value 0.008030203669192642 0.0001 1e-06
+@value -0.00308977207424735 0.0001 1e-06
+@value -0.003271431276926969 0.0001 1e-06
 Dictionary Key: prop2_radial_center
-@value -6.670678078067882e-05 0.0001 1e-06
-@value -0.0001453592285835322 0.0001 1e-06
-@value -0.0002331773818241109 0.0001 1e-06
-@value -0.0002601552627753417 0.0001 1e-06
-@value -0.0002161746327240625 0.0001 1e-06
-@value -0.000207624441553162 0.0001 1e-06
-@value -7.99645718944244e-05 0.0001 1e-06
-@value 3.444546730535738e-05 0.0001 1e-06
-@value  0.00196047593830362 0.0001 1e-06
-@value -0.004910846120141911 0.0001 1e-06
+@value -6.6821618253851e-05 0.0001 1e-06
+@value -0.0001454882142161054 0.0001 1e-06
+@value -0.0002333925777396506 0.0001 1e-06
+@value -0.0002602360346694545 0.0001 1e-06
+@value -0.0002160029253573562 0.0001 1e-06
+@value -0.0002073763249032774 0.0001 1e-06
+@value -7.983905119491938e-05 0.0001 1e-06
+@value 3.447463394741616e-05 0.0001 1e-06
+@value 0.001959038767328681 0.0001 1e-06
+@value -0.004906987993166209 0.0001 1e-06
 Dictionary Key: prop2_tangential_force
-@value 2.070750511408624e-07 0.0001 1e-06
-@value 1.177415520853327e-06 0.0001 1e-06
-@value 1.780246815423283e-06 0.0001 1e-06
-@value 5.691461043814424e-07 0.0001 1e-06
-@value -1.306604447675903e-06 0.0001 1e-06
-@value -6.338285115092114e-06 0.0001 1e-06
-@value -9.782960820114565e-07 0.0001 1e-06
-@value 1.307384461601081e-06 0.0001 1e-06
-@value -4.401509715648327e-06 0.0001 1e-06
-@value  1.3435558024864e-05 0.0001 1e-06
+@value 2.049849076550651e-07 0.0001 1e-06
+@value 1.173913574489503e-06 0.0001 1e-06
+@value 1.773252334941589e-06 0.0001 1e-06
+@value 5.701534304018501e-07 0.0001 1e-06
+@value -1.296653096217944e-06 0.0001 1e-06
+@value -6.322243147720968e-06 0.0001 1e-06
+@value -9.73426024981711e-07 0.0001 1e-06
+@value 1.301186301973147e-06 0.0001 1e-06
+@value -4.381349853564366e-06 0.0001 1e-06
+@value 1.338121839647789e-05 0.0001 1e-06
 Dictionary Key: shape
-@value  -0.2804055754834046 0.0001 1e-06
-@value  -0.3703737724347647 0.0001 1e-06
-@value  -0.2819511215047202 0.0001 1e-06
-@value   0.2887923696127854 0.0001 1e-06
-@value   0.3520148373231896 0.0001 1e-06
-@value   0.2883501714591917 0.0001 1e-06
-@value   0.1219989579707799 0.0001 1e-06
-@value   0.1555816296663317 0.0001 1e-06
-@value   0.1224181288529918 0.0001 1e-06
-@value   0.5503591594534432 0.0001 1e-06
-@value   0.6576414579938006 0.0001 1e-06
-@value   0.5508993251866839 0.0001 1e-06
-@value  -0.9546042735464592 0.0001 1e-06
-@value   -1.211733305356477 0.0001 1e-06
-@value  -0.9568433782447888 0.0001 1e-06
-@value   0.3358913304395638 0.0001 1e-06
-@value   0.4081272090287227 0.0001 1e-06
-@value   0.3373916756405173 0.0001 1e-06
+@value  -0.2713860581536602 0.0001 1e-06
+@value  -0.3589274885811068 0.0001 1e-06
+@value  -0.2729226892840378 0.0001 1e-06
+@value   0.2879939767344152 0.0001 1e-06
+@value    0.351155861543735 0.0001 1e-06
+@value   0.2875434218867562 0.0001 1e-06
+@value   0.1210663519143854 0.0001 1e-06
+@value   0.1544736007165708 0.0001 1e-06
+@value   0.1214818207282285 0.0001 1e-06
+@value   0.5470950724945915 0.0001 1e-06
+@value   0.6538247942485816 0.0001 1e-06
+@value   0.5476222523946299 0.0001 1e-06
+@value   -0.957447373574652 0.0001 1e-06
+@value   -1.215210144966994 0.0001 1e-06
+@value  -0.9596928903100117 0.0001 1e-06
+@value   0.3325186103593699 0.0001 1e-06
+@value   0.4041131665245409 0.0001 1e-06
+@value   0.3340096038487914 0.0001 1e-06
\ No newline at end of file
diff --git a/tests/runTests_MphysAeroField.py b/tests/runTests_MphysAeroField.py
index b452356b..58e61b1c 100755
--- a/tests/runTests_MphysAeroField.py
+++ b/tests/runTests_MphysAeroField.py
@@ -63,7 +63,7 @@
     "normalizeStates": {"U": U0, "p": U0 * U0 / 2.0, "phi": 1.0, "nuTilda": 1e-3},
     "adjPartDerivFDStep": {"State": 1e-6},
     "designVar": {
-        "shape": {"designVarType": "FFD"},
+        "twist": {"designVarType": "FFD"},
         "aoa": {"designVarType": "AOA", "patches": ["inout"], "flowAxis": "x", "normalAxis": "y"},
         "alphaPorosity": {"designVarType": "Field", "fieldName": "alphaPorosity", "fieldType": "scalar"},
     },
@@ -132,24 +132,27 @@ def alphaPorosity(val, DASolver):
         self.cruise.coupling.solver.add_dv_func("alphaPorosity", alphaPorosity)
         self.cruise.aero_post.add_dv_func("alphaPorosity", alphaPorosity)
 
-        # Select all points
-        pts = self.geometry.DVGeo.getLocalIndex(0)
-        indexList = pts[:, :, :].flatten()
-        PS = geo_utils.PointSelect("list", indexList)
-        nShapes = self.geometry.nom_addLocalDV(dvName="shape", pointSelect=PS)
+        # Set up global design variables. We dont change the root twist
+        nRefAxPts = self.geometry.nom_addRefAxis(name="wingAxis", xFraction=0.25, alignIndex="k")
+
+        def twist(val, geo):
+            for i in range(nRefAxPts):
+                geo.rot_z["wingAxis"].coef[i] = -val[0]
+
+        self.geometry.nom_addGlobalDV(dvName="twist", value=np.array([0]), func=twist)
 
         # add dvs to ivc and connect
-        self.dvs.add_output("shape", val=np.array([0] * nShapes))
+        self.dvs.add_output("twist", val=np.array([0]))
         self.dvs.add_output("aoa", val=np.array([aoa0]))
         nLocalCells = self.cruise.coupling.DASolver.solver.getNLocalCells()
         self.dvs.add_output("alphaPorosity", val=np.zeros(nLocalCells, dtype="d"))
 
-        self.connect("shape", "geometry.shape")
+        self.connect("twist", "geometry.twist")
         self.connect("aoa", "cruise.aoa")
         self.connect("alphaPorosity", "cruise.alphaPorosity")
 
         # define the design variables
-        self.add_design_var("shape", lower=-1.0, upper=1.0, scaler=1.0)
+        self.add_design_var("twist", lower=-10.0, upper=10.0, scaler=1.0)
         self.add_design_var("aoa", lower=-10.0, upper=10.0, scaler=1.0)
         self.add_design_var("alphaPorosity", lower=0, upper=1e-4, scaler=1.0)
 
@@ -191,14 +194,16 @@ def alphaPorosity(val, DASolver):
 prob.run_model()
 totals = prob.compute_totals()
 
-
 alphaNorm = np.linalg.norm(totals[("cruise.aero_post.functionals.CL", "dvs.alphaPorosity")].flatten())
 alphaNormSum = gcomm.allreduce(alphaNorm, op=MPI.SUM)
 
 if gcomm.rank == 0:
+    objFuncDict = {}
+    objFuncDict["CL"] = prob.get_val("cruise.aero_post.CL")
     derivDict = {}
     derivDict["CL"] = {}
     derivDict["CL"]["alpha"] = [alphaNormSum]
-    derivDict["CL"]["shape"] = totals[("cruise.aero_post.functionals.CL", "dvs.shape")][0]
+    derivDict["CL"]["twist"] = totals[("cruise.aero_post.functionals.CL", "dvs.twist")][0]
     derivDict["CL"]["aoa"] = totals[("cruise.aero_post.functionals.CL", "dvs.aoa")][0]
+    reg_write_dict(objFuncDict, 1e-6, 1e-8)
     reg_write_dict(derivDict, 1e-4, 1e-6)
diff --git a/tests/runTests_MphysAeroOpt.py b/tests/runTests_MphysAeroOpt.py
index b8c559d6..2ed97a7c 100755
--- a/tests/runTests_MphysAeroOpt.py
+++ b/tests/runTests_MphysAeroOpt.py
@@ -309,8 +309,8 @@ def actuator(val, DASolver):
         finalCL[key] = finalCon[key]
 
 if gcomm.rank == 0:
-    reg_write_dict(finalObj, 1e-6, 1e-8)
-    reg_write_dict(finalCL, 1e-6, 1e-8)
+    reg_write_dict(finalObj, 1e-4, 1e-6)
+    reg_write_dict(finalCL, 1e-4, 1e-6)
 
 # test the find feasible design function
 optFuncs.findFeasibleDesign(