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.

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));

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)

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

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();
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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();
 }

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_);

src/adjoint/DASolver/DARhoSimpleFoam/pEqnRhoSimple.H

@@ -47,15 +47,15 @@ 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())
     {
         phi = phiHbyA + pEqn.flux();
     }
 }
 
-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_);

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

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_);

src/adjoint/DASolver/DASimpleFoam/pEqnSimple.H

@@ -49,15 +49,15 @@
         // and final residuals
         SolverPerformance<scalar> solverP = pEqn.solve();
 
-        DAUtility::primalResidualControl(solverP, printToScreen, "p", primalMaxRes);
+        DAUtility::primalResidualControl(solverP, printToScreen_, "p", primalMaxRes_);
 
         if (simple.finalNonOrthogonalIter())
         {
             phi = phiHbyA - pEqn.flux();
         }
     }
 
-    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);
 }