From f4b0ccd7219023613b1465a042c7daeeb85595f4 Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Fri, 20 Dec 2024 09:59:30 -0600
Subject: [PATCH] Added DARhoSimpleFoam (#726)

* Fixed a bug for the distributed flag for functions

* Made printInterval and maxRes class-wide vars.

* Added the DARhoSimpleFoam tests.

* Fixed test files.

* Added the missing file.
---
 .../DAResidual/DAResidualRhoSimpleFoam.C      |   7 +-
 src/adjoint/DAResidual/Make/files             |   1 +
 .../DARhoSimpleFoam/DARhoSimpleFoam.C         |  89 +--------
 .../DARhoSimpleFoam/DARhoSimpleFoam.H         |   4 +-
 .../DASolver/DARhoSimpleFoam/EEqnRhoSimple.H  |   4 +-
 .../DASolver/DARhoSimpleFoam/UEqnRhoSimple.H  |   6 +-
 .../DASolver/DARhoSimpleFoam/pEqnRhoSimple.H  |  10 +-
 .../DASolver/DASimpleFoam/DASimpleFoam.C      |  19 +-
 .../DASolver/DASimpleFoam/UEqnSimple.H        |   4 +-
 .../DASolver/DASimpleFoam/pEqnSimple.H        |   8 +-
 src/adjoint/DASolver/DASolver.C               |  12 +-
 src/adjoint/DASolver/DASolver.H               |  13 +-
 src/adjoint/DASolver/Make/files               |   2 +-
 src/adjoint/DAStateInfo/Make/files            |   1 +
 src/include/createAdjoint.H                   |  44 +++++
 tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt |  26 +--
 tests/runRegTests_DARhoSimpleFoam.py          | 174 ++++++++++++++++++
 tests/runRegTests_DASimpleFoam.py             |   2 +-
 18 files changed, 281 insertions(+), 145 deletions(-)
 create mode 100755 src/include/createAdjoint.H
 create mode 100755 tests/runRegTests_DARhoSimpleFoam.py

diff --git a/src/adjoint/DAResidual/DAResidualRhoSimpleFoam.C b/src/adjoint/DAResidual/DAResidualRhoSimpleFoam.C
index a688018f..af6149e4 100755
--- a/src/adjoint/DAResidual/DAResidualRhoSimpleFoam.C
+++ b/src/adjoint/DAResidual/DAResidualRhoSimpleFoam.C
@@ -33,7 +33,8 @@ DAResidualRhoSimpleFoam::DAResidualRhoSimpleFoam(
       fvSourceEnergy_(const_cast<volScalarField&>(
           mesh_.thisDb().lookupObject<volScalarField>("fvSourceEnergy"))),
       fvOptions_(fv::options::New(mesh)),
-      thermo_(const_cast<fluidThermo&>(daModel.getThermo())),
+      thermo_(const_cast<fluidThermo&>(
+          mesh_.thisDb().lookupObject<fluidThermo>("thermophysicalProperties"))),
       he_(thermo_.he()),
       rho_(const_cast<volScalarField&>(
           mesh_.thisDb().lookupObject<volScalarField>("rho"))),
@@ -70,7 +71,7 @@ DAResidualRhoSimpleFoam::DAResidualRhoSimpleFoam(
         Info << "Cp " << Cp_ << endl;
     }
 
-    // this is just a dummy call because we need to run the constrain once 
+    // this is just a dummy call because we need to run the constrain once
     // to initialize fvOptions, before we can use it. Otherwise, we may get
     // a seg fault when we call fvOptions_.correct(U_) in updateIntermediateVars
     fvVectorMatrix UEqn(
@@ -193,7 +194,7 @@ void DAResidualRhoSimpleFoam::calcResiduals(const dictionary& options)
         HbyAPtr() = rAU * UEqn.H();
     }
     volVectorField& HbyA = HbyAPtr();
-    
+
     tUEqn.clear();
 
     surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho_) * fvc::flux(HbyA));
diff --git a/src/adjoint/DAResidual/Make/files b/src/adjoint/DAResidual/Make/files
index b1750277..47b28e89 100755
--- a/src/adjoint/DAResidual/Make/files
+++ b/src/adjoint/DAResidual/Make/files
@@ -1,4 +1,5 @@
 DAResidual.C
 DAResidualSimpleFoam.C
+DAResidualRhoSimpleFoam.C
 
 LIB = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedLibs/libDAResidual$(WM_CODI_AD_LIB_POSTFIX)
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
index b4e7e68e..bcbebcea 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.C
@@ -73,13 +73,13 @@ void DARhoSimpleFoam::initSolver()
     argList& args = argsPtr_();
 #include "createSimpleControlPython.H"
 #include "createFieldsRhoSimple.H"
-#include "createAdjointCompressible.H"
+#include "createAdjoint.H"
     // initialize checkMesh
     daCheckMeshPtr_.reset(new DACheckMesh(daOptionPtr_(), runTime, mesh));
 
     daLinearEqnPtr_.reset(new DALinearEqn(mesh, daOptionPtr_()));
 
-    this->setDAObjFuncList();
+    this->setDAFunctionList();
 
     // initialize fvSource and compute the source term
     const dictionary& allOptions = daOptionPtr_->getAllOptions();
@@ -94,19 +94,11 @@ void DARhoSimpleFoam::initSolver()
     }
 }
 
-label DARhoSimpleFoam::solvePrimal(
-    const Vec xvVec,
-    Vec wVec)
+label DARhoSimpleFoam::solvePrimal()
 {
     /*
     Description:
         Call the primal solver to get converged state variables
-
-    Input:
-        xvVec: a vector that contains all volume mesh coordinates
-
-    Output:
-        wVec: state variable vector
     */
 
 #include "createRefsRhoSimple.H"
@@ -121,44 +113,10 @@ label DARhoSimpleFoam::solvePrimal(
     Info << "\nStarting time loop\n"
          << endl;
 
-    // deform the mesh based on the xvVec
-    this->pointVec2OFMesh(xvVec);
-
-    // check mesh quality
-    label meshOK = this->checkMesh();
-
-    if (!meshOK)
-    {
-        this->writeFailedMesh();
-        return 1;
-    }
-
-    // if the forwardModeAD is active, we need to set the seed here
-#include "setForwardADSeeds.H"
-
-    word divUScheme = "div(phi,U)";
-    if (daOptionPtr_->getSubDictOption<label>("runLowOrderPrimal4PC", "active"))
-    {
-        if (daOptionPtr_->getSubDictOption<label>("runLowOrderPrimal4PC", "isPC"))
-        {
-            Info << "Using low order div scheme for primal solution .... " << endl;
-            divUScheme = "div(pc)";
-        }
-    }
-
-    // if useMeanStates is used, we need to zero meanStates before the primal run
-    this->zeroMeanStates();
-
-    label printInterval = daOptionPtr_->getOption<label>("printInterval");
-    label printToScreen = 0;
-    label regModelFail = 0;
     while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
     {
-        DAUtility::primalMaxInitRes_ = -1e16;
 
-        printToScreen = this->isPrintTime(runTime, printInterval);
-
-        if (printToScreen)
+        if (printToScreen_)
         {
             Info << "Time = " << runTime.timeName() << nl << endl;
         }
@@ -171,60 +129,29 @@ label DARhoSimpleFoam::solvePrimal(
 #include "EEqnRhoSimple.H"
 #include "pEqnRhoSimple.H"
 
-        daTurbulenceModelPtr_->correct(printToScreen);
-
-        // update the output field value at each iteration, if the regression model is active
-        regModelFail = daRegressionPtr_->compute();
+        daTurbulenceModelPtr_->correct(printToScreen_, primalMaxRes_);
 
-        if (this->validateStates())
-        {
-            // write data to files and quit
-            runTime.writeNow();
-            mesh.write();
-            return 1;
-        }
-
-        if (printToScreen)
+        if (printToScreen_)
         {
             daTurbulenceModelPtr_->printYPlus();
 
-            this->printAllObjFuncs();
-
-            daRegressionPtr_->printInputInfo();
+            this->printAllFunctions();
 
             Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                  << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                  << nl << endl;
         }
 
-        // if useMeanStates is used, we need to calculate the meanStates
-        this->calcMeanStates();
-
         runTime.write();
     }
 
-    if (regModelFail != 0)
-    {
-        return 1;
-    }
-
-    // if useMeanStates is used, we need to assign meanStates to states right after the case converges
-    this->assignMeanStatesToStates();
-
-    this->writeAssociatedFields();
-
-    this->calcPrimalResidualStatistics("print");
-
-    // primal converged, assign the OpenFoam fields to the state vec wVec
-    this->ofField2StateVec(wVec);
-
     // write the mesh to files
     mesh.write();
 
     Info << "End\n"
          << endl;
 
-    return this->checkResidualTol();
+    return this->checkResidualTol(primalMaxRes_);
 }
 
 } // End namespace Foam
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.H b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.H
index 46761168..884b8bb9 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.H
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/DARhoSimpleFoam.H
@@ -123,9 +123,7 @@ public:
     virtual void initSolver();
 
     /// solve the primal equations
-    virtual label solvePrimal(
-        const Vec xvVec,
-        Vec wVec);
+    virtual label solvePrimal();
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/EEqnRhoSimple.H b/src/adjoint/DASolver/DARhoSimpleFoam/EEqnRhoSimple.H
index f17ac6de..0123a74d 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/EEqnRhoSimple.H
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/EEqnRhoSimple.H
@@ -21,10 +21,10 @@
     // and final residuals
     SolverPerformance<scalar> solverE = EEqn.solve();
 
-    DAUtility::primalResidualControl(solverE, printToScreen, "he");
+    DAUtility::primalResidualControl(solverE, printToScreen_, "he", primalMaxRes_);
 
     // bound he
-    DAUtility::boundVar(allOptions, he, printToScreen);
+    DAUtility::boundVar(allOptions, he, printToScreen_);
 
     thermo.correct();
 }
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H b/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H
index 10c04ba7..494767a2 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/UEqnRhoSimple.H
@@ -8,7 +8,7 @@ if (hasFvSource_)
 }
 
 tmp<fvVectorMatrix> tUEqn(
-    fvm::div(phi, U, divUScheme)
+    fvm::div(phi, U)
     + MRF.DDt(rho, U)
     + turbulencePtr_->divDevRhoReff(U)
     - fvSource
@@ -23,9 +23,9 @@ fvOptions.constrain(UEqn);
 // and final residuals
 SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-DAUtility::primalResidualControl(solverU, printToScreen, "U");
+DAUtility::primalResidualControl(solverU, printToScreen_, "U", primalMaxRes_);
 
 fvOptions.correct(U);
 
 // bound U
-DAUtility::boundVar(allOptions, U, printToScreen);
+DAUtility::boundVar(allOptions, U, printToScreen_);
diff --git a/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H b/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H
index e100626e..91e3035b 100755
--- a/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H
+++ b/src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H
@@ -47,7 +47,7 @@ while (simple.correctNonOrthogonal())
     // and final residuals
     SolverPerformance<scalar> solverP = pEqn.solve();
 
-    DAUtility::primalResidualControl(solverP, printToScreen, "p");
+    DAUtility::primalResidualControl(solverP, printToScreen_, "p", primalMaxRes_);
 
     if (simple.finalNonOrthogonalIter())
     {
@@ -55,7 +55,7 @@ while (simple.correctNonOrthogonal())
     }
 }
 
-if (printToScreen)
+if (printToScreen_)
 {
 #include "continuityErrsPython.H"
 }
@@ -64,11 +64,11 @@ if (printToScreen)
 p.relax();
 
 // bound p
-DAUtility::boundVar(allOptions, p, printToScreen);
+DAUtility::boundVar(allOptions, p, printToScreen_);
 
 U = HbyA - rAU * fvc::grad(p);
 // bound U
-DAUtility::boundVar(allOptions, U, printToScreen);
+DAUtility::boundVar(allOptions, U, printToScreen_);
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 
@@ -92,4 +92,4 @@ rho = thermo.rho();
 rho.relax();
 
 // bound rho
-DAUtility::boundVar(allOptions, rho, printToScreen);
+DAUtility::boundVar(allOptions, rho, printToScreen_);
diff --git a/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
index 2604acec..e28d65f7 100755
--- a/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
+++ b/src/adjoint/DASolver/DASimpleFoam/DASimpleFoam.C
@@ -69,7 +69,7 @@ void DASimpleFoam::initSolver()
     fvMesh& mesh = meshPtr_();
 #include "createSimpleControlPython.H"
 #include "createFieldsSimple.H"
-#include "createAdjointIncompressible.H"
+#include "createAdjoint.H"
     // initialize checkMesh
     daCheckMeshPtr_.reset(new DACheckMesh(daOptionPtr_(), runTime, mesh));
 
@@ -109,16 +109,11 @@ label DASimpleFoam::solvePrimal()
     Info << "\nStarting time loop\n"
          << endl;
 
-    label printInterval = daOptionPtr_->getOption<label>("printInterval");
-    label printToScreen = 0;
-    scalar primalMaxRes = 0.0;
-    while (this->loop(runTime, primalMaxRes)) // using simple.loop() will have seg fault in parallel
+    
+    while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
     {
-        primalMaxRes = -1e10;
-
-        printToScreen = this->isPrintTime(runTime, printInterval);
 
-        if (printToScreen)
+        if (printToScreen_)
         {
             Info << "Time = " << runTime.timeName() << nl << endl;
         }
@@ -132,9 +127,9 @@ label DASimpleFoam::solvePrimal()
         }
 
         laminarTransport.correct();
-        daTurbulenceModelPtr_->correct(printToScreen, primalMaxRes);
+        daTurbulenceModelPtr_->correct(printToScreen_, primalMaxRes_);
 
-        if (printToScreen)
+        if (printToScreen_)
         {
             daTurbulenceModelPtr_->printYPlus();
 
@@ -154,7 +149,7 @@ label DASimpleFoam::solvePrimal()
     Info << "End\n"
          << endl;
 
-    return this->checkResidualTol(primalMaxRes);
+    return this->checkResidualTol(primalMaxRes_);
 }
 
 // ************ the following are functions for consistent fixed-point adjoint
diff --git a/src/adjoint/DASolver/DASimpleFoam/UEqnSimple.H b/src/adjoint/DASolver/DASimpleFoam/UEqnSimple.H
index 01a81dd9..ba769839 100755
--- a/src/adjoint/DASolver/DASimpleFoam/UEqnSimple.H
+++ b/src/adjoint/DASolver/DASimpleFoam/UEqnSimple.H
@@ -26,10 +26,10 @@ if (simple.momentumPredictor())
     // and final residuals
     SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-    DAUtility::primalResidualControl(solverU, printToScreen, "U", primalMaxRes);
+    DAUtility::primalResidualControl(solverU, printToScreen_, "U", primalMaxRes_);
 
     fvOptions.correct(U);
 }
 
 // bound U
-DAUtility::boundVar(allOptions, U, printToScreen);
+DAUtility::boundVar(allOptions, U, printToScreen_);
diff --git a/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H b/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H
index 5aab2da5..ce2c273e 100755
--- a/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H
+++ b/src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H
@@ -49,7 +49,7 @@
         // and final residuals
         SolverPerformance<scalar> solverP = pEqn.solve();
 
-        DAUtility::primalResidualControl(solverP, printToScreen, "p", primalMaxRes);
+        DAUtility::primalResidualControl(solverP, printToScreen_, "p", primalMaxRes_);
 
         if (simple.finalNonOrthogonalIter())
         {
@@ -57,7 +57,7 @@
         }
     }
 
-    if (printToScreen)
+    if (printToScreen_)
     {
 #include "continuityErrsPython.H"
     }
@@ -66,12 +66,12 @@
     p.relax();
 
     // bound p
-    DAUtility::boundVar(allOptions, p, printToScreen);
+    DAUtility::boundVar(allOptions, p, printToScreen_);
 
     // Momentum corrector
     U = HbyA - rAtU() * fvc::grad(p);
     // bound U
-    DAUtility::boundVar(allOptions, U, printToScreen);
+    DAUtility::boundVar(allOptions, U, printToScreen_);
     U.correctBoundaryConditions();
     fvOptions.correct(U);
 }
diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index d5e65a57..4bbfbee1 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -64,6 +64,7 @@ DASolver::DASolver(
 
     primalMinResTol_ = daOptionPtr_->getOption<scalar>("primalMinResTol");
     primalMinIters_ = daOptionPtr_->getOption<label>("primalMinIters");
+    printInterval_ = daOptionPtr_->getOption<label>("printInterval");
 
     // initialize the objStd variables.
     this->initObjStd();
@@ -116,9 +117,7 @@ autoPtr<DASolver> DASolver::New(
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-label DASolver::loop(
-    Time& runTime,
-    const scalar& primalMaxRes)
+label DASolver::loop(Time& runTime)
 {
     /*
     Description:
@@ -145,11 +144,11 @@ label DASolver::loop(
     this->calcObjStd(runTime);
 
     // check exit condition, we need to satisfy both the residual and function std condition
-    if (primalMaxRes < primalMinResTol_ && runTime.timeIndex() > primalMinIters_ && primalObjStd_ < primalObjStdTol_)
+    if (primalMaxRes_ < primalMinResTol_ && runTime.timeIndex() > primalMinIters_ && primalObjStd_ < primalObjStdTol_)
     {
         Info << "Time = " << t << endl;
 
-        Info << "Minimal residual " << primalMaxRes << " satisfied the prescribed tolerance " << primalMinResTol_ << endl
+        Info << "Minimal residual " << primalMaxRes_ << " satisfied the prescribed tolerance " << primalMinResTol_ << endl
              << endl;
 
         if (primalObjStdActive_)
@@ -175,6 +174,9 @@ label DASolver::loop(
     else
     {
         ++runTime;
+        // initialize primalMaxRes_ with a small value for this iteration
+        primalMaxRes_ = -1e10;
+        printToScreen_ = this->isPrintTime(runTime, printInterval_);
         return 1;
     }
 }
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index 88f20a7e..3b7fa0da 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -189,6 +189,15 @@ protected:
     /// primal residual tolerance
     scalar primalMinResTol_ = 0.0;
 
+    /// primal maximal residual computed from all the governing equations
+    scalar primalMaxRes_ = 0.0;
+
+    /// whether to print primal information to the screen
+    label printToScreen_ = 0;
+
+    /// how frequent do you want to print the primal info default is every 100 steps  
+    label printInterval_ = 100;
+
     /// primal min number of iterations
     label primalMinIters_ = -1;
 
@@ -657,9 +666,7 @@ public:
     void initializeGlobalADTape4dRdWT();
 
     /// return whether to loop the primal solution, similar to runTime::loop() except we don't do file IO
-    label loop(
-        Time& runTime,
-        const scalar& primalMaxRes);
+    label loop(Time& runTime);
 
     /// basically, we call DAIndex::getGlobalXvIndex
     label getGlobalXvIndex(
diff --git a/src/adjoint/DASolver/Make/files b/src/adjoint/DASolver/Make/files
index 5b6f8030..7e08f0eb 100755
--- a/src/adjoint/DASolver/Make/files
+++ b/src/adjoint/DASolver/Make/files
@@ -1,5 +1,5 @@
 DASolver.C
 DASimpleFoam/DASimpleFoam.C
-
+DARhoSimpleFoam/DARhoSimpleFoam.C
 
 LIB = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedLibs/libDASolver$(WM_CODI_AD_LIB_POSTFIX)
diff --git a/src/adjoint/DAStateInfo/Make/files b/src/adjoint/DAStateInfo/Make/files
index cfd16011..f8f17958 100755
--- a/src/adjoint/DAStateInfo/Make/files
+++ b/src/adjoint/DAStateInfo/Make/files
@@ -1,4 +1,5 @@
 DAStateInfo.C
 DAStateInfoSimpleFoam.C
+DAStateInfoRhoSimpleFoam.C
 
 LIB = $(DAFOAM_ROOT_PATH)/OpenFOAM/sharedLibs/libDAStateInfo$(WM_CODI_AD_LIB_POSTFIX)
diff --git a/src/include/createAdjoint.H b/src/include/createAdjoint.H
new file mode 100755
index 00000000..81399478
--- /dev/null
+++ b/src/include/createAdjoint.H
@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+
+    DAFoam  : Discrete Adjoint with OpenFOAM
+    Version : v4
+
+\*---------------------------------------------------------------------------*/
+
+// adjoint
+
+// read the RAS model from constant/turbulenceProperties
+const word turbModelName(
+    IOdictionary(
+        IOobject(
+            "turbulenceProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false))
+        .subDict("RAS")
+        .lookup("RASModel"));
+daTurbulenceModelPtr_.reset(DATurbulenceModel::New(turbModelName, mesh, daOptionPtr_()));
+
+daModelPtr_.reset(new DAModel(mesh, daOptionPtr_()));
+
+daRegressionPtr_.reset(new DARegression(mesh, daOptionPtr_(), daModelPtr_()));
+
+word solverName = daOptionPtr_->getOption<word>("solverName");
+daStateInfoPtr_.reset(DAStateInfo::New(solverName, mesh, daOptionPtr_(), daModelPtr_()));
+stateInfo_ = daStateInfoPtr_->getStateInfo();
+
+daIndexPtr_.reset(new DAIndex(mesh, daOptionPtr_(), daModelPtr_()));
+
+daIndexPtr_->printIndices();
+
+daFieldPtr_.reset(new DAField(mesh, daOptionPtr_(), daModelPtr_(), daIndexPtr_()));
+
+daResidualPtr_.reset(DAResidual::New(solverName, mesh, daOptionPtr_(), daModelPtr_(), daIndexPtr_()));
+
+// calculate the initial (default) value for validate state func
+this->getInitStateVals(initStateVals_);
+
+// initialize the mean states if needed
+this->initMeanStates();
diff --git a/tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt b/tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt
index 30e75f1a..c01a3189 100755
--- a/tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt
@@ -1,24 +1,10 @@
 Dictionary Key: CD
-@value  0.09084626022005476 1e-08 1e-10
+@value   0.0139069043904323 1e-10 1e-12
 Dictionary Key: CL
-@value   0.3304878903763182 1e-08 1e-10
-Dictionary Key: fail
-@value                    0 1e-08 1e-10
+@value   0.2551422843394988 1e-10 1e-12
 Dictionary Key: CD
-Dictionary Key: shapey
-@value   0.2125351330069156 1e-04 1e-06
-@value   0.1780152462220735 1e-04 1e-06
-@value    2.376356190880319 1e-04 1e-06
-@value     1.88949857276801 1e-04 1e-06
-Dictionary Key: uin
-@value 0.003747752666385211 1e-04 1e-06
+Dictionary Key: twist
+@value 0.001053063900009023 1e-08 1e-12
 Dictionary Key: CL
-Dictionary Key: shapey
-@value   0.5355908807229974 1e-04 1e-06
-@value   0.4425505807722847 1e-04 1e-06
-@value   0.4449741622842366 1e-04 1e-06
-@value   0.3619240288380361 1e-04 1e-06
-Dictionary Key: uin
-@value  0.01343529283859845 1e-04 1e-06
-Dictionary Key: fail
-@value                    0 1e-04 1e-06
\ No newline at end of file
+Dictionary Key: twist
+@value   0.1261578915912536 1e-08 1e-12
\ No newline at end of file
diff --git a/tests/runRegTests_DARhoSimpleFoam.py b/tests/runRegTests_DARhoSimpleFoam.py
new file mode 100755
index 00000000..2ef3de6b
--- /dev/null
+++ b/tests/runRegTests_DARhoSimpleFoam.py
@@ -0,0 +1,174 @@
+#!/usr/bin/env python
+"""
+Run Python tests for optimization integration
+"""
+
+from mpi4py import MPI
+import os
+import numpy as np
+from testFuncs import *
+
+import openmdao.api as om
+from mphys.multipoint import Multipoint
+from dafoam.mphys import DAFoamBuilder, OptFuncs
+from mphys.scenario_aerodynamic import ScenarioAerodynamic
+from pygeo.mphys import OM_DVGEOCOMP
+
+gcomm = MPI.COMM_WORLD
+
+os.chdir("./reg_test_files-main/NACA0012")
+if gcomm.rank == 0:
+    os.system("rm -rf 0 processor* *.bin")
+    os.system("cp -r 0.compressible 0")
+    os.system("cp -r system.subsonic system")
+    os.system("cp -r constant/turbulenceProperties.sa constant/turbulenceProperties")
+    replace_text_in_file("system/fvSchemes", "meshWave;", "meshWaveFrozen;")
+
+# aero setup
+U0 = 100.0
+p0 = 101325.0
+T0 = 300.0
+A0 = 0.1
+twist0 = 2.0
+LRef = 1.0
+nuTilda0 = 4.5e-5
+
+daOptions = {
+    "designSurfaces": ["wing"],
+    "solverName": "DARhoSimpleFoam",
+    "primalMinResTol": 1.0e-11,
+    "primalMinResTolDiff": 1e4,
+    "primalBC": {
+        "U0": {"variable": "U", "patches": ["inout"], "value": [U0, 0.0, 0.0]},
+        "T0": {"variable": "T", "patches": ["inout"], "value": [T0]},
+        "p0": {"variable": "p", "patches": ["inout"], "value": [p0]},
+        "nuTilda0": {"variable": "nuTilda", "patches": ["inout"], "value": [nuTilda0]},
+        "useWallFunction": True,
+    },
+    "function": {
+        "CD": {
+            "part1": {
+                "type": "force",
+                "source": "patchToFace",
+                "patches": ["wing"],
+                "directionMode": "fixedDirection",
+                "direction": [1.0, 0.0, 0.0],
+                "scale": 1.0 / (0.5 * U0 * U0 * A0),
+                "addToAdjoint": True,
+            }
+        },
+        "CL": {
+            "part1": {
+                "type": "force",
+                "source": "patchToFace",
+                "patches": ["wing"],
+                "directionMode": "fixedDirection",
+                "direction": [0.0, 1.0, 0.0],
+                "scale": 1.0 / (0.5 * U0 * U0 * A0),
+                "addToAdjoint": True,
+            }
+        },
+    },
+    "adjEqnOption": {
+        "gmresRelTol": 1.0e-11,
+        "pcFillLevel": 1,
+        "jacMatReOrdering": "rcm",
+        "dynAdjustTol": False
+    },
+    "normalizeStates": {"U": U0, "p": p0, "phi": 1.0, "T": T0, "nuTilda": 1e-3},
+    "designVar": {
+        "twist": {"designVarType": "FFD"},
+    },
+}
+
+meshOptions = {
+    "gridFile": os.getcwd(),
+    "fileType": "OpenFOAM",
+    # point and normal for the symmetry plane
+    "symmetryPlanes": [[[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]], [[0.0, 0.0, 0.1], [0.0, 0.0, 1.0]]],
+}
+
+
+class Top(Multipoint):
+    def setup(self):
+        dafoam_builder = DAFoamBuilder(daOptions, meshOptions, scenario="aerodynamic")
+        dafoam_builder.initialize(self.comm)
+
+        ################################################################################
+        # MPHY setup
+        ################################################################################
+
+        # ivc to keep the top level DVs
+        self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"])
+
+        # create the mesh and cruise scenario because we only have one analysis point
+        self.add_subsystem("mesh", dafoam_builder.get_mesh_coordinate_subsystem())
+
+        # add the geometry component, we dont need a builder because we do it here.
+        self.add_subsystem("geometry", OM_DVGEOCOMP(file="FFD/wingFFD.xyz", type="ffd"))
+
+        self.mphys_add_scenario("cruise", ScenarioAerodynamic(aero_builder=dafoam_builder))
+
+        self.connect("mesh.x_aero0", "geometry.x_aero_in")
+        self.connect("geometry.x_aero0", "cruise.x_aero")
+
+        self.add_subsystem("LoD", om.ExecComp("val=CL/CD"))
+
+    def configure(self):
+
+        self.cruise.aero_post.mphys_add_funcs()
+
+        # create geometric DV setup
+        points = self.mesh.mphys_get_surface_mesh()
+
+        # add pointset
+        self.geometry.nom_add_discipline_coords("aero", points)
+
+        # geometry setup
+        self.geometry.nom_addRefAxis(name="wingAxis", xFraction=0.25, alignIndex="k")
+
+        # Set up global design variables. We dont change the root twist
+        def twist(val, geo):
+            for i in range(2):
+                geo.rot_z["wingAxis"].coef[i] = -val[0]
+
+        self.geometry.nom_addGlobalDV(dvName="twist", value=np.ones(1) * twist0, func=twist)
+
+        # add the design variables to the dvs component's output
+        self.dvs.add_output("twist", val=np.ones(1) * twist0)
+        # manually connect the dvs output to the geometry and cruise
+        self.connect("twist", "geometry.twist")
+
+        # define the design variables to the top level
+        self.add_design_var("twist", lower=-10.0, upper=10.0, scaler=1.0)
+
+        # add constraints and the objective
+        self.add_objective("cruise.aero_post.CD", scaler=1.0)
+        self.add_constraint("cruise.aero_post.CL", equals=0.3)
+
+
+prob = om.Problem()
+prob.model = Top()
+
+prob.setup(mode="rev")
+om.n2(prob, show_browser=False, outfile="mphys_aero.html")
+
+optFuncs = OptFuncs(daOptions, prob)
+
+# verify the total derivatives against the finite-difference
+prob.run_model()
+results = prob.check_totals(
+    of=["cruise.aero_post.CD", "cruise.aero_post.CL"], wrt=["twist"], compact_print=True, step=1e-3, form="central", step_calc="abs"
+)
+
+if gcomm.rank == 0:
+    funcDict = {}
+    funcDict["CD"] = prob.get_val("cruise.aero_post.CD")
+    funcDict["CL"] = prob.get_val("cruise.aero_post.CL")
+    derivDict = {}
+    derivDict["CD"] = {}
+    derivDict["CD"]["twist"] = results[("cruise.aero_post.CD", "twist")]["J_fwd"][0]
+    derivDict["CL"] = {}
+    derivDict["CL"]["twist"] = results[("cruise.aero_post.CL", "twist")]["J_fwd"][0]
+    reg_write_dict(funcDict, 1e-10, 1e-12)
+    reg_write_dict(derivDict, 1e-8, 1e-12)
diff --git a/tests/runRegTests_DASimpleFoam.py b/tests/runRegTests_DASimpleFoam.py
index e6b7dfd4..5e453e34 100755
--- a/tests/runRegTests_DASimpleFoam.py
+++ b/tests/runRegTests_DASimpleFoam.py
@@ -18,7 +18,7 @@
 
 os.chdir("./reg_test_files-main/NACA0012")
 if gcomm.rank == 0:
-    os.system("rm -rf 0 processor*")
+    os.system("rm -rf 0 processor* *.bin")
     os.system("cp -r 0.incompressible 0")
     os.system("cp -r system.incompressible system")
     os.system("cp -r constant/turbulenceProperties.sa constant/turbulenceProperties")