Added pimple adjoint and verified totals (#737)

* Added DAPimpleFoam and primal working.

* Updated the pimple ref.

* Added pimple and verified derivs

* Fixed a critical bug about normRes.

* Added the updatePC call back.

* Updated ref.

* Updated simpleFoam refs.

* Added shape var for pimple. Totals not accurate yet.

* Tested pimple shape with forward and the totals are accurate.

dafoam/pyDAFoam.py

@@ -589,7 +589,7 @@ def __init__(self):
         self.printInterval = 100
 
         ## The print interval of unsteady primal solvers, e.g., for DAPisoFoam
-        self.printIntervalUnsteady = 500
+        self.printIntervalUnsteady = 1
 
         ## Users can adjust primalMinResTolDiff to tweak how much difference between primalMinResTol
         ## and the actual primal convergence is consider to be fail=True for the primal solution.
@@ -1526,7 +1526,68 @@ def readStateVars(self, timeVal, deltaT):
 
         # assign the state from OF field to wVec so that the wVec
         # is update to date for unsteady adjoint
-        self.solver.ofField2StateVec(self.wVec)
+        # self.solver.ofField2StateVec(self.wVec)
+
+    def set_solver_input(self, inputs, DVGeo=None):
+        """
+        Set solver input. If it is forward mode, we also set the seeds
+        """
+        inputDict = self.getOption("solverInput")
+
+        for inputName in list(inputDict.keys()):
+            inputType = inputDict[inputName]["type"]
+            input = inputs[inputName]
+            inputSize = len(input)
+            seeds = np.zeros(inputSize)
+            if self.getOption("useAD")["mode"] == "forward":
+                if inputType == "volCoord":
+                    seeds = self.calcFFD2XvSeeds(DVGeo)
+                elif inputName == self.getOption("useAD")["dvName"]:
+                    seedIndex = self.getOption("useAD")["seedIndex"]
+                    seeds[seedIndex] = 1.0
+            # here we need to update the solver input for both solver and solverAD
+            self.solver.setSolverInput(inputName, inputType, inputSize, input, seeds)
+            self.solverAD.setSolverInput(inputName, inputType, inputSize, input, seeds)
+
+    def calcFFD2XvSeeds(self, DVGeo=None):
+        # Calculate the FFD2XvSeed array:
+        # Given a FFD seed xDvDot, run pyGeo and IDWarp and propagate the seed to Xv seed xVDot:
+        #     xSDot = \\frac{dX_{S}}{dX_{DV}}\\xDvDot
+        #     xVDot = \\frac{dX_{V}}{dX_{S}}\\xSDot
+
+        # Then, we assign this vector to FFD2XvSeed in the mphys_dafoam
+        # This will be used in forward mode AD runs
+
+        discipline = self.getOption("discipline")
+
+        if DVGeo is None:
+            raise Error("calcFFD2XvSeeds is call but no DVGeo object found! Call add_dvgeo in the run script!")
+
+        if self.mesh is None:
+            raise Error("calcFFD2XvSeeds is call but no mesh object found!")
+
+        dvName = self.getOption("useAD")["dvName"]
+        seedIndex = self.getOption("useAD")["seedIndex"]
+        # create xDVDot vec and initialize it with zeros
+        xDV = DVGeo.getValues()
+
+        # create a copy of xDV and set the seed to 1.0
+        # the dv and index depends on dvName and seedIndex
+        xDvDot = {}
+        for key in list(xDV.keys()):
+            xDvDot[key] = np.zeros_like(xDV[key])
+        xDvDot[dvName][seedIndex] = 1.0
+
+        # get the original surf coords
+        xSDot0 = np.zeros_like(self.xs0)
+        xSDot0 = self.mapVector(xSDot0, self.allSurfacesGroup, self.designSurfacesGroup)
+
+        # get xSDot
+        xSDot = DVGeo.totalSensitivityProd(xDvDot, ptSetName="x_%s0" % discipline).reshape(xSDot0.shape)
+        # get xVDot
+        xVDot = self.mesh.warpDerivFwd(xSDot)
+
+        return xVDot
 
     def solveAdjointUnsteady(self):
         """
@@ -3016,6 +3077,41 @@ def convertMPIVec2SeqArray(self, mpiVec):
 
         return array1
 
+    def arrayVal2Vec(self, array1, vec):
+        """
+        Assign the values from array1 to vec
+        """
+
+        size = len(array1)
+
+        Istart, Iend = vec.getOwnershipRange()
+
+        if (Iend - Istart) != size:
+            raise Error("array and vec's sizes are not consistent")
+
+        for i in range(Istart, Iend):
+            iRel = i - Istart
+            vec[i] = array1[iRel]
+
+        vec.assemblyBegin()
+        vec.assemblyEnd()
+
+    def vecVal2Array(self, vec, array1):
+        """
+        Assign the values from vec to array1
+        """
+
+        size = len(array1)
+
+        Istart, Iend = vec.getOwnershipRange()
+
+        if (Iend - Istart) != size:
+            raise Error("array and vec's sizes are not consistent")
+
+        for i in range(Istart, Iend):
+            iRel = i - Istart
+            array1[iRel] = vec[i]
+
     def vec2Array(self, vec):
         """
         Convert a Petsc vector to numpy array

src/adjoint/DAInput/DAInputVolCoord.C

@@ -57,6 +57,9 @@ void DAInputVolCoord::run(const scalarList& input)
         }
     }
     mesh_.movePoints(meshPoints);
+    // set this to false to avoid V0 not found error
+    // TODO. We need to change this for dynamic mesh cases
+    mesh_.moving(false);
 
     return;
 }

src/adjoint/DAResidual/Make/files

@@ -1,6 +1,7 @@
 DAResidual.C
 DAResidualSimpleFoam.C
 DAResidualSimpleTFoam.C
+DAResidualPimpleFoam.C
 DAResidualRhoSimpleFoam.C
 DAResidualRhoSimpleCFoam.C
 DAResidualHeatTransferFoam.C

src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C

@@ -67,13 +67,22 @@ void DAPimpleFoam::initSolver()
     fvMesh& mesh = meshPtr_();
 #include "createPimpleControlPython.H"
 #include "createFieldsPimple.H"
-#include "createAdjointIncompressible.H"
-    // initialize checkMesh
-    daCheckMeshPtr_.reset(new DACheckMesh(daOptionPtr_(), runTime, mesh));
 
-    daLinearEqnPtr_.reset(new DALinearEqn(mesh, daOptionPtr_()));
-
-    this->setDAObjFuncList();
+    // 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_()));
+
+#include "createAdjoint.H"
 
     // initialize fvSource and the source term
     const dictionary& allOptions = daOptionPtr_->getAllOptions();
@@ -101,69 +110,25 @@ void DAPimpleFoam::initSolver()
     }
 }
 
-label DAPimpleFoam::solvePrimal(
-    const Vec xvVec,
-    Vec wVec)
+label DAPimpleFoam::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 "createRefsPimple.H"
 
-    // change the run status
-    daOptionPtr_->setOption<word>("runStatus", "solvePrimal");
-
-    runTime.setTime(0.0, 0);
-    // if readZeroFields, we need to read in the states from the 0 folder every time 
-    // we start the primal here we read in all time levels
-    label readZeroFields = daOptionPtr_->getAllOptions().subDict("unsteadyAdjoint").getLabel("readZeroFields");
-    if (readZeroFields)
-    {
-        this->readStateVars(0.0, 0);
-        this->readStateVars(0.0, 1);
-    }
-
     // call correctNut, this is equivalent to turbulence->validate();
     daTurbulenceModelPtr_->updateIntermediateVariables();
 
     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"
-
-    // We need to set the mesh moving to false, otherwise we will get V0 not found error.
-    // Need to dig into this issue later
-    // NOTE: we have commented this out. Setting mesh.moving(false) has been done
-    // right after mesh.movePoints() calls.
-    //mesh.moving(false);
-
-    label printInterval = daOptionPtr_->getOption<label>("printIntervalUnsteady");
-    label printToScreen = 0;
     label pimplePrintToScreen = 0;
 
     // reset the unsteady obj func to zeros
-    this->initUnsteadyObjFuncs();
+    this->initUnsteadyFunctions();
 
     // we need to reduce the number of files written to the disk to minimize the file IO load
     label reduceIO = daOptionPtr_->getAllOptions().subDict("unsteadyAdjoint").getLabel("reduceIO");
@@ -178,24 +143,21 @@ label DAPimpleFoam::solvePrimal(
     label nInstances = round(endTime / deltaT);
 
     // main loop
-    label regModelFail = 0;
-    label fail = 0;
     for (label iter = 1; iter <= nInstances; iter++)
     {
-
         ++runTime;
 
-        printToScreen = this->isPrintTime(runTime, printInterval);
+        printToScreen_ = this->isPrintTime(runTime, printIntervalUnsteady_);
 
-        if (printToScreen)
+        if (printToScreen_)
         {
             Info << "Time = " << runTime.timeName() << nl << endl;
         }
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            if (pimple.finalIter() && printToScreen)
+            if (pimple.finalIter() && printToScreen_)
             {
                 pimplePrintToScreen = 1;
             }
@@ -213,39 +175,25 @@ label DAPimpleFoam::solvePrimal(
             }
 
             laminarTransport.correct();
-            daTurbulenceModelPtr_->correct(pimplePrintToScreen);
-
-            // update the output field value at each iteration, if the regression model is active
-            fail = daRegressionPtr_->compute();
-        }
-
-        regModelFail += fail;
-
-        if (this->validateStates())
-        {
-            // write data to files and quit
-            runTime.writeNow();
-            mesh.write();
-            return 1;
+            // primalMaxRes_ is just a dummy input, we will not use it
+            daTurbulenceModelPtr_->correct(pimplePrintToScreen, primalMaxRes_);
         }
 
-        this->calcUnsteadyObjFuncs();
+        this->calcUnsteadyFunctions();
 
-        if (printToScreen)
+        if (printToScreen_)
         {
 #include "CourantNo.H"
 
             daTurbulenceModelPtr_->printYPlus();
 
-            this->printAllObjFuncs();
+            this->printAllFunctions();
 
             if (daOptionPtr_->getOption<label>("debug"))
             {
                 this->calcPrimalResidualStatistics("print");
             }
 
-            daRegressionPtr_->printInputInfo();
-
             Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                  << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                  << nl << endl;
@@ -254,25 +202,13 @@ label DAPimpleFoam::solvePrimal(
         if (reduceIO && iter < nInstances)
         {
             this->writeAdjStates(reduceIOWriteMesh_, additionalOutput);
-            daRegressionPtr_->writeFeatures();
         }
         else
         {
             runTime.write();
-            daRegressionPtr_->writeFeatures();
         }
     }
 
-    if (regModelFail > 0)
-    {
-        return 1;
-    }
-
-    this->calcPrimalResidualStatistics("print");
-
-    // primal converged, assign the OpenFoam fields to the state vec wVec
-    this->ofField2StateVec(wVec);
-
     // write the mesh to files
     mesh.write();
 

src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.H

@@ -119,9 +119,7 @@ public:
     virtual void initSolver();
 
     /// solve the primal equations
-    virtual label solvePrimal(
-        const Vec xvVec,
-        Vec wVec);
+    virtual label solvePrimal();
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DASolver/DAPimpleFoam/UEqnPimple.H

@@ -22,8 +22,8 @@ if (pimple.momentumPredictor())
     // and final residuals
     SolverPerformance<vector> solverU = solve(UEqn == -fvc::grad(p));
 
-    DAUtility::primalResidualControl(solverU, pimplePrintToScreen, "U");
+    DAUtility::primalResidualControl(solverU, pimplePrintToScreen, "U", primalMaxRes_);
 }
 
 // bound U
-DAUtility::boundVar(allOptions, U, printToScreen);
+DAUtility::boundVar(allOptions, U, pimplePrintToScreen);

src/adjoint/DASolver/DAPimpleFoam/pEqnPimple.H

@@ -52,7 +52,7 @@ while (pimple.correctNonOrthogonal())
     // and final residuals
     SolverPerformance<scalar> solverP = pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
-    DAUtility::primalResidualControl(solverP, pimplePrintToScreen, "p");
+    DAUtility::primalResidualControl(solverP, pimplePrintToScreen, "p", primalMaxRes_);
 
     if (pimple.finalNonOrthogonalIter())
     {
@@ -69,9 +69,9 @@ if (pimplePrintToScreen)
 p.relax();
 
 // bound p
-DAUtility::boundVar(allOptions, p, printToScreen);
+DAUtility::boundVar(allOptions, p, pimplePrintToScreen);
 
 U = HbyA - rAtU() * fvc::grad(p);
 // bound U
-DAUtility::boundVar(allOptions, U, printToScreen);
+DAUtility::boundVar(allOptions, U, pimplePrintToScreen);
 U.correctBoundaryConditions();