From 09d98c2d7ec5be4449a17d9f091bc2330e288315 Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Mon, 15 Jan 2024 21:11:17 -0600
Subject: [PATCH] Updated the solid solver.

---
 .../DAObjFunc/DAObjFuncVonMisesStressKS.C     |  9 +--
 .../DAResidualSolidDisplacementFoam.C         | 75 ++-----------------
 .../DAResidualSolidDisplacementFoam.H         |  8 +-
 .../DASolidDisplacementFoam.C                 | 56 ++++----------
 .../DASolidDisplacementFoam.H                 | 20 -----
 .../calculateStressSolidDisplacement.H        | 40 +++++-----
 .../createControlsSolidDisplacement.H         |  7 --
 .../createFieldsSolidDisplacement.H           | 38 ----------
 .../createRefsSolidDisplacement.H             |  2 -
 ...eadMechanicalPropertiesSolidDisplacement.H | 16 ++--
 .../tractionDisplacementFvPatchVectorField.C  | 37 ++-------
 ...DAFoam_Test_DASolidDisplacementFoamRef.txt | 26 +++----
 tests/runTests_DASolidDisplacementFoam.py     |  9 ++-
 13 files changed, 81 insertions(+), 262 deletions(-)
 delete mode 100644 src/adjoint/DASolver/DASolidDisplacementFoam/createControlsSolidDisplacement.H

diff --git a/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C b/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
index de67494e..111c05bb 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
@@ -46,7 +46,6 @@ DAObjFuncVonMisesStressKS::DAObjFuncVonMisesStressKS(
     objFuncDict_.readEntry<scalar>("scale", scale_);
 
     objFuncDict_.readEntry<scalar>("coeffKS", coeffKS_);
-
 }
 
 /// calculate the value of objective function
@@ -87,16 +86,16 @@ void DAObjFuncVonMisesStressKS::calcObjFunc(
     objFuncValue = 0.0;
 
     const objectRegistry& db = mesh_.thisDb();
-    //const volVectorField& D = db.lookupObject<volVectorField>("D");
+    const volTensorField& gradD = db.lookupObject<volTensorField>("gradD");
     const volScalarField& lambda = db.lookupObject<volScalarField>("solid:lambda");
     const volScalarField& mu = db.lookupObject<volScalarField>("solid:mu");
     const volScalarField& rho = db.lookupObject<volScalarField>("solid:rho");
-    const volTensorField& gradD = db.lookupObject<volTensorField>("gradD");
 
-    volSymmTensorField sigma = rho * (mu * twoSymm(gradD) + lambda * (I * tr(gradD)));
+    // volTensorField gradD(fvc::grad(D));
+    volSymmTensorField sigma = rho * (mu * twoSymm(gradD) + (lambda * I) * tr(gradD));
 
     // NOTE: vonMises stress is scaled by scale_ provided in the objFunc dict
-    volScalarField vonMises = scale_* sqrt((3.0 / 2.0) * magSqr(dev(sigma)));
+    volScalarField vonMises = scale_ * sqrt((3.0 / 2.0) * magSqr(dev(sigma)));
 
     scalar objValTmp = 0.0;
     forAll(objFuncCellSources, idxI)
diff --git a/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.C b/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.C
index 9488ee43..4e061d70 100755
--- a/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.C
+++ b/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.C
@@ -28,47 +28,11 @@ DAResidualSolidDisplacementFoam::DAResidualSolidDisplacementFoam(
       // these are intermediate variables or objects
       gradD_(const_cast<volTensorField&>(
           mesh.thisDb().lookupObject<volTensorField>("gradD"))),
-      sigmaD_(const_cast<volSymmTensorField&>(
-          mesh.thisDb().lookupObject<volSymmTensorField>("sigmaD"))),
-      divSigmaExp_(const_cast<volVectorField&>(
-          mesh.thisDb().lookupObject<volVectorField>("divSigmaExp"))),
       lambda_(const_cast<volScalarField&>(
           mesh.thisDb().lookupObject<volScalarField>("solid:lambda"))),
       mu_(const_cast<volScalarField&>(
           mesh.thisDb().lookupObject<volScalarField>("solid:mu")))
 {
-
-    IOdictionary thermalProperties(
-        IOobject(
-            "thermalProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE));
-
-    Switch thermalStress(thermalProperties.lookup("thermalStress"));
-    if (thermalStress)
-    {
-        FatalErrorIn("") << "thermalStress=true not supported" << abort(FatalError);
-    }
-
-    const dictionary& stressControl = mesh.solutionDict().subDict("stressAnalysis");
-
-    Switch compactNormalStress(stressControl.lookup("compactNormalStress"));
-
-    if (!compactNormalStress)
-    {
-        FatalErrorIn("") << "compactNormalStress=false not supported" << abort(FatalError);
-    }
-
-    isTractionDisplacementBC_ = 0;
-    forAll(D_.boundaryField(), patchI)
-    {
-        if (D_.boundaryField()[patchI].type() == "tractionDisplacement")
-        {
-            isTractionDisplacementBC_ = 1;
-        }
-    }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -99,42 +63,19 @@ void DAResidualSolidDisplacementFoam::calcResiduals(const dictionary& options)
         URes_, pRes_, TRes_, phiRes_: residual field variables
     */
 
+    volSymmTensorField sigmaD = mu_ * twoSymm(gradD_) + (lambda_ * I) * tr(gradD_);
+
+    volVectorField divSigmaExp = fvc::div(sigmaD - (2 * mu_ + lambda_) * gradD_, "div(sigmaD)");
+
     fvVectorMatrix DEqn(
         fvm::d2dt2(D_)
         == fvm::laplacian(2 * mu_ + lambda_, D_, "laplacian(DD,D)")
-            + divSigmaExp_);
+            + divSigmaExp);
 
     DRes_ = DEqn & D_;
     normalizeResiduals(DRes);
 }
 
-void DAResidualSolidDisplacementFoam::updateDAndGradD()
-{
-    /*
-    Description:
-        Update D and gradD.
-
-        NOTE: we need to update D boundary conditions iteratively if tractionDisplacement BC 
-        is used this is because tractionDisplacement BC is dependent on gradD, while gradD
-        is dependent on the D bc values.
-    */
-
-    // this will be called after doing perturbStates
-    if (isTractionDisplacementBC_)
-    {
-        for (label i = 0; i < daOption_.getOption<label>("maxTractionBCIters"); i++)
-        {
-            D_.correctBoundaryConditions();
-            gradD_ = fvc::grad(D_);
-        }
-    }
-    else
-    {
-        D_.correctBoundaryConditions();
-        gradD_ = fvc::grad(D_);
-    }
-}
-
 void DAResidualSolidDisplacementFoam::updateIntermediateVariables()
 {
     /* 
@@ -142,11 +83,7 @@ void DAResidualSolidDisplacementFoam::updateIntermediateVariables()
         Update the intermediate variables that depend on the state variables
     */
 
-    this->updateDAndGradD();
-
-    sigmaD_ = mu_ * twoSymm(gradD_) + (lambda_ * I) * tr(gradD_);
-
-    divSigmaExp_ = fvc::div(sigmaD_ - (2 * mu_ + lambda_) * gradD_, "div(sigmaD)");
+    gradD_ = fvc::grad(D_);
 }
 
 void DAResidualSolidDisplacementFoam::correctBoundaryConditions()
diff --git a/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.H b/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.H
index d3d3361c..a47ca9ef 100755
--- a/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.H
+++ b/src/adjoint/DAResidual/DAResidualSolidDisplacementFoam.H
@@ -28,22 +28,16 @@ class DAResidualSolidDisplacementFoam
 {
 
 protected:
+
     /// \name These are state variables, state residuals, and partial derivatives
     //@{
     volVectorField& D_;
     volVectorField DRes_;
     //@}
-
     volTensorField& gradD_;
-    volSymmTensorField& sigmaD_;
-    volVectorField& divSigmaExp_;
     volScalarField& lambda_;
     volScalarField& mu_;
 
-    label isTractionDisplacementBC_;
-
-    void updateDAndGradD();
-
 public:
     TypeName("DASolidDisplacementFoam");
     // Constructors
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C b/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C
index c1405ea3..21f18f3d 100755
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C
+++ b/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.C
@@ -42,16 +42,13 @@ DASolidDisplacementFoam::DASolidDisplacementFoam(
     char* argsAll,
     PyObject* pyOptions)
     : DASolver(argsAll, pyOptions),
-      mechanicalPropertiesPtr_(nullptr),
       rhoPtr_(nullptr),
       muPtr_(nullptr),
       lambdaPtr_(nullptr),
       EPtr_(nullptr),
       nuPtr_(nullptr),
       DPtr_(nullptr),
-      sigmaDPtr_(nullptr),
-      gradDPtr_(nullptr),
-      divSigmaExpPtr_(nullptr)
+      gradDPtr_(nullptr)
 {
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -66,7 +63,6 @@ void DASolidDisplacementFoam::initSolver()
     Info << "Initializing fields for DASolidDisplacementFoam" << endl;
     Time& runTime = runTimePtr_();
     fvMesh& mesh = meshPtr_();
-#include "createControlsSolidDisplacement.H"
 #include "createFieldsSolidDisplacement.H"
 #include "createAdjointSolid.H"
     // initialize checkMesh
@@ -125,43 +121,21 @@ label DASolidDisplacementFoam::solvePrimal(
             Info << "Iteration = " << runTime.value() << nl << endl;
         }
 
-        label iCorr = 0;
-        scalar initialResidual = 0;
+        gradD = fvc::grad(D);
+        volSymmTensorField sigmaD = mu * twoSymm(gradD) + (lambda * I) * tr(gradD);
 
-        do
-        {
-            fvVectorMatrix DEqn(
-                fvm::d2dt2(D)
-                == fvm::laplacian(2 * mu + lambda, D, "laplacian(DD,D)")
-                    + divSigmaExp);
-
-            // get the solver performance info such as initial
-            // and final residuals
-            SolverPerformance<vector> solverU = DEqn.solve();
-            initialResidual = solverU.max().initialResidual();
-
-            this->primalResidualControl<vector>(solverU, printToScreen, printInterval, "U");
-
-            if (!compactNormalStress_)
-            {
-                divSigmaExp = fvc::div(DEqn.flux());
-            }
-
-            gradD = fvc::grad(D);
-            sigmaD = mu * twoSymm(gradD) + (lambda * I) * tr(gradD);
-
-            if (compactNormalStress_)
-            {
-                divSigmaExp = fvc::div(
-                    sigmaD - (2 * mu + lambda) * gradD,
-                    "div(sigmaD)");
-            }
-            else
-            {
-                divSigmaExp += fvc::div(sigmaD);
-            }
-
-        } while (initialResidual > convergenceTolerance_ && ++iCorr < nCorr_);
+        volVectorField divSigmaExp = fvc::div(sigmaD - (2 * mu + lambda) * gradD, "div(sigmaD)");
+
+        fvVectorMatrix DEqn(
+            fvm::d2dt2(D)
+            == fvm::laplacian(2 * mu + lambda, D, "laplacian(DD,D)")
+                + divSigmaExp);
+
+        // get the solver performance info such as initial
+        // and final residuals
+        SolverPerformance<vector> solverD = DEqn.solve();
+
+        this->primalResidualControl<vector>(solverD, printToScreen, printInterval, "D");
 
         if (this->validateStates())
         {
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.H b/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.H
index c54f8ab5..f428e369 100755
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.H
+++ b/src/adjoint/DASolver/DASolidDisplacementFoam/DASolidDisplacementFoam.H
@@ -52,18 +52,6 @@ class DASolidDisplacementFoam
 
 protected:
 
-    /// whether to use compact normal stress
-    Switch compactNormalStress_;
-
-    /// number of correctors
-    label nCorr_;
-
-    /// convergence tolerance for D
-    scalar convergenceTolerance_;
-
-    /// mechanicalProperties pointer
-    autoPtr<IOdictionary> mechanicalPropertiesPtr_;
-
     /// density field pointer
     autoPtr<volScalarField> rhoPtr_;
 
@@ -81,17 +69,10 @@ protected:
 
     /// displacement field pointer
     autoPtr<volVectorField> DPtr_;
-    
-    /// stress field pointer
-    autoPtr<volSymmTensorField> sigmaDPtr_;
 
     /// displacement gradient pointer
     autoPtr<volTensorField> gradDPtr_;
 
-    /// explicit part of div(sigma) pointer
-    autoPtr<volVectorField> divSigmaExpPtr_;
-
-    
 public:
     TypeName("DASolidDisplacementFoam");
     // Constructors
@@ -113,7 +94,6 @@ public:
     virtual label solvePrimal(
         const Vec xvVec,
         Vec wVec);
-
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/calculateStressSolidDisplacement.H b/src/adjoint/DASolver/DASolidDisplacementFoam/calculateStressSolidDisplacement.H
index caed0ec8..7fb9cf08 100644
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/calculateStressSolidDisplacement.H
+++ b/src/adjoint/DASolver/DASolidDisplacementFoam/calculateStressSolidDisplacement.H
@@ -1,23 +1,25 @@
+{
+    volSymmTensorField sigmaD = mu * twoSymm(gradD) + (lambda * I) * tr(gradD);
 
-volSymmTensorField sigma(
-    IOobject(
-        "sigma",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE),
-    rho* sigmaD);
+    volSymmTensorField sigma(
+        IOobject(
+            "sigma",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE),
+        rho * sigmaD);
 
-volScalarField sigmaEq(
-    IOobject(
-        "sigmaEq",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE),
-    sqrt((3.0 / 2.0) * magSqr(dev(sigma))));
+    volScalarField sigmaEq(
+        IOobject(
+            "sigmaEq",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE),
+        sqrt((3.0 / 2.0) * magSqr(dev(sigma))));
 
-Info << "Max sigmaEq = " << max(sigmaEq).value()
-     << endl;
+    Info << "Max sigmaEq = " << gMax(sigmaEq) << endl;
 
-runTime.write();
+    runTime.write();
+}
\ No newline at end of file
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/createControlsSolidDisplacement.H b/src/adjoint/DASolver/DASolidDisplacementFoam/createControlsSolidDisplacement.H
deleted file mode 100644
index a8a428bf..00000000
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/createControlsSolidDisplacement.H
+++ /dev/null
@@ -1,7 +0,0 @@
-const dictionary& stressControl = mesh.solutionDict().subDict("stressAnalysis");
-
-stressControl.readEntry<Switch>("compactNormalStress", compactNormalStress_);
-
-nCorr_ = stressControl.lookupOrDefault<label>("nCorrectors", 1);
-
-convergenceTolerance_ = readScalar(stressControl.lookup("D"));
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/createFieldsSolidDisplacement.H b/src/adjoint/DASolver/DASolidDisplacementFoam/createFieldsSolidDisplacement.H
index 056c5cbf..9f9a9380 100644
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/createFieldsSolidDisplacement.H
+++ b/src/adjoint/DASolver/DASolidDisplacementFoam/createFieldsSolidDisplacement.H
@@ -13,20 +13,6 @@ DPtr_.reset(
         mesh));
 volVectorField& D = DPtr_();
 
-Info << "Calculating stress field sigmaD\n"
-     << endl;
-
-sigmaDPtr_.reset(
-    new volSymmTensorField(
-        IOobject(
-            "sigmaD",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE),
-        muPtr_() * twoSymm(fvc::grad(D)) + lambdaPtr_() * (I * tr(fvc::grad(D)))));
-volSymmTensorField& sigmaD = sigmaDPtr_();
-
 // gradD is used in the tractionDisplacement BC
 gradDPtr_.reset(
     new volTensorField(
@@ -37,27 +23,3 @@ gradDPtr_.reset(
             IOobject::NO_READ,
             IOobject::NO_WRITE),
         fvc::grad(D)));
-
-Info << "Calculating explicit part of div(sigma) divSigmaExp\n"
-     << endl;
-divSigmaExpPtr_.reset(
-    new volVectorField(
-    IOobject(
-        "divSigmaExp",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE),
-    fvc::div(sigmaD)));
-volVectorField& divSigmaExp = divSigmaExpPtr_();
-
-if (compactNormalStress_)
-{
-    divSigmaExp -= fvc::laplacian(2 * muPtr_() + lambdaPtr_(), D, "laplacian(DD,D)");
-}
-else
-{
-    divSigmaExp -= fvc::div((2 * muPtr_() + lambdaPtr_()) * fvc::grad(D), "div(sigmaD)");
-}
-
-mesh.setFluxRequired(D.name());
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/createRefsSolidDisplacement.H b/src/adjoint/DASolver/DASolidDisplacementFoam/createRefsSolidDisplacement.H
index 952ace22..cd458fb1 100755
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/createRefsSolidDisplacement.H
+++ b/src/adjoint/DASolver/DASolidDisplacementFoam/createRefsSolidDisplacement.H
@@ -6,6 +6,4 @@ volScalarField& rho = rhoPtr_();
 volScalarField& mu = muPtr_();
 volScalarField& lambda = lambdaPtr_();
 volVectorField& D = DPtr_();
-volSymmTensorField& sigmaD = sigmaDPtr_();
 volTensorField& gradD = gradDPtr_();
-volVectorField& divSigmaExp = divSigmaExpPtr_();
diff --git a/src/adjoint/DASolver/DASolidDisplacementFoam/readMechanicalPropertiesSolidDisplacement.H b/src/adjoint/DASolver/DASolidDisplacementFoam/readMechanicalPropertiesSolidDisplacement.H
index 15e4f3d9..9a921b0e 100644
--- a/src/adjoint/DASolver/DASolidDisplacementFoam/readMechanicalPropertiesSolidDisplacement.H
+++ b/src/adjoint/DASolver/DASolidDisplacementFoam/readMechanicalPropertiesSolidDisplacement.H
@@ -1,15 +1,13 @@
 Info << "Reading mechanical properties\n"
      << endl;
 
-mechanicalPropertiesPtr_.reset(
-    new IOdictionary(
-        IOobject(
-            "mechanicalProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE)));
-IOdictionary& mechanicalProperties = mechanicalPropertiesPtr_();
+IOdictionary mechanicalProperties(
+    IOobject(
+        "mechanicalProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE));
 
 const dictionary& rhoDict(mechanicalProperties.subDict("rho"));
 word rhoType(rhoDict.lookup("type"));
diff --git a/src/adjoint/boundaryConditions/tractionDisplacement/tractionDisplacementFvPatchVectorField.C b/src/adjoint/boundaryConditions/tractionDisplacement/tractionDisplacementFvPatchVectorField.C
index 4b214265..21288f26 100644
--- a/src/adjoint/boundaryConditions/tractionDisplacement/tractionDisplacementFvPatchVectorField.C
+++ b/src/adjoint/boundaryConditions/tractionDisplacement/tractionDisplacementFvPatchVectorField.C
@@ -63,7 +63,6 @@ tractionDisplacementFvPatchVectorField::
 {
 }
 
-
 tractionDisplacementFvPatchVectorField::
     tractionDisplacementFvPatchVectorField(
         const fvPatch& p,
@@ -105,7 +104,6 @@ void tractionDisplacementFvPatchVectorField::autoMap(
     fixedGradientFvPatchVectorField::autoMap(m);
     traction_.autoMap(m);
     pressure_.autoMap(m);
-
 }
 
 void tractionDisplacementFvPatchVectorField::rmap(
@@ -120,7 +118,6 @@ void tractionDisplacementFvPatchVectorField::rmap(
 
     traction_.rmap(dmptf.traction_, addr);
     pressure_.rmap(dmptf.pressure_, addr);
-
 }
 
 void tractionDisplacementFvPatchVectorField::updateCoeffs()
@@ -130,39 +127,21 @@ void tractionDisplacementFvPatchVectorField::updateCoeffs()
         return;
     }
 
-    const dictionary& mechanicalProperties =
-        db().lookupObject<IOdictionary>("mechanicalProperties");
-
     const fvPatchField<scalar>& rho =
         patch().lookupPatchField<volScalarField, scalar>("solid:rho");
 
-    const fvPatchField<scalar>& rhoE =
-        patch().lookupPatchField<volScalarField, scalar>("E");
-
-    const fvPatchField<scalar>& nu =
-        patch().lookupPatchField<volScalarField, scalar>("solid:nu");
-
-    scalarField E(rhoE / rho);
-    scalarField mu(E / (2.0 * (1.0 + nu)));
-    scalarField lambda(nu * E / ((1.0 + nu) * (1.0 - 2.0 * nu)));
-
-    Switch planeStress(mechanicalProperties.lookup("planeStress"));
+    const fvPatchField<scalar>& lambda =
+        patch().lookupPatchField<volScalarField, scalar>("solid:lambda");
 
-    if (planeStress)
-    {
-        lambda = nu * E / ((1.0 + nu) * (1.0 - nu));
-    }
+    const fvPatchField<scalar>& mu =
+        patch().lookupPatchField<volScalarField, scalar>("solid:mu");
 
-    scalarField twoMuLambda(2 * mu + lambda);
-
-    vectorField n(patch().nf());
-
-    // NOTE: we can't use the built-in gradient() computation because
-    // it is designed for transient problem, i.e., snGrad() is actually
-    // related to gradient()
-    // Here we implement our owned BC for steady-state
     const fvPatchField<tensor>& gradD =
         patch().lookupPatchField<volTensorField, tensor>("gradD");
+
+    // here we use the BC implementation from:
+    // Tang, Tian: Implementation of solid body stress analysis in OpenFOAM, 2013
+    vectorField n(patch().nf());
     gradient() = ((traction_ - pressure_ * n) / rho
                   - (n & (mu * gradD.T() - (mu + lambda) * gradD))
                   - n * tr(gradD) * lambda)
diff --git a/tests/refs/DAFoam_Test_DASolidDisplacementFoamRef.txt b/tests/refs/DAFoam_Test_DASolidDisplacementFoamRef.txt
index ce299941..e5a7f7fa 100755
--- a/tests/refs/DAFoam_Test_DASolidDisplacementFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DASolidDisplacementFoamRef.txt
@@ -1,26 +1,26 @@
 Dictionary Key: M
 @value    11950.18306038372 1e-08 1e-10
 Dictionary Key: VMS
-@value    29161.94589567093 1e-08 1e-10
+@value    29161.94589567075 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
 Dictionary Key: M
 Dictionary Key: shapex
-@value   -134.4035951515252 0.0001 1e-06
-@value   -355.6352710947976 0.0001 1e-06
-@value   -134.3171787448227 0.0001 1e-06
+@value    -39.4467911094272 0.0001 1e-06
+@value   -104.4851626918655 0.0001 1e-06
+@value   -39.44679110942718 0.0001 1e-06
 Dictionary Key: shapey
-@value    133.5780999106646 0.0001 1e-06
-@value    351.9087686872808 0.0001 1e-06
-@value    133.8342145372735 0.0001 1e-06
+@value    39.44679110942718 0.0001 1e-06
+@value    104.4851626918655 0.0001 1e-06
+@value    39.44679110942717 0.0001 1e-06
 Dictionary Key: VMS
 Dictionary Key: shapex
-@value   -992.7386180744397 0.01 1e-06
-@value   -2587.884010719518 0.01 1e-06
-@value   -1264.630576690301 0.01 1e-06
+@value   -981.9681824399274 0.0001 1e-06
+@value   -1120.804160309113 0.0001 1e-06
+@value   -981.8260104013304 0.0001 1e-06
 Dictionary Key: shapey
-@value   -740.4371760828722 0.01 1e-06
-@value   -6358.686218955751 0.01 1e-06
-@value   -1162.331263670126 0.01 1e-06
+@value   -368.8831651062262 0.0001 1e-06
+@value   -5129.234283204046 0.0001 1e-06
+@value   -368.8587906435345 0.0001 1e-06
 Dictionary Key: fail
 @value                    0 0.0001 1e-06
\ No newline at end of file
diff --git a/tests/runTests_DASolidDisplacementFoam.py b/tests/runTests_DASolidDisplacementFoam.py
index 15f74939..a524394d 100755
--- a/tests/runTests_DASolidDisplacementFoam.py
+++ b/tests/runTests_DASolidDisplacementFoam.py
@@ -22,16 +22,19 @@
 gcomm = MPI.COMM_WORLD
 
 os.chdir("./input/PlateHole")
+if gcomm.rank == 0:
+    os.system("rm -rf processor*")
 
 # test incompressible solvers
 aeroOptions = {
     "debug": True,
-    "maxTractionBCIters": 20,
     "solverName": "DASolidDisplacementFoam",
-    "useAD": {"mode": "fd"},
+    "useAD": {"mode": "reverse"},
     "designSurfaces": ["hole"],
     "primalMinResTol": 1e-10,
     "primalMinResTolDiff": 1e10,
+    "hasIterativeBC": True,
+    "maxCorrectBCCalls": 20,
     "objFunc": {
         "VMS": {
             "part1": {
@@ -57,7 +60,7 @@
     },
     "normalizeStates": {"D": 1.0e-7},
     "adjPartDerivFDStep": {"State": 1e-5, "FFD": 1e-3},
-    "adjEqnOption": {"gmresRelTol": 1.0e-15, "gmresAbsTol": 1.0e-15, "pcFillLevel": 1, "jacMatReOrdering": "rcm"},
+    "adjEqnOption": {"gmresRelTol": 1.0e-12, "gmresAbsTol": 1.0e-12, "pcFillLevel": 1, "jacMatReOrdering": "rcm"},
     # Design variable setup
     "designVar": {
         "shapey": {"designVarType": "FFD"},