Enabled forward mode AD. Verified machine-precision accurate adjoint …

…for DASimpleFoam

dafoam/mphys/mphys_dafoam.py

@@ -707,6 +707,43 @@ def apply_options(self, optionDict):
                 self.DASolver.setOption(key, optionDict[key])
             self.DASolver.updateDAOption()
 
+    def calcFFD2XvSeeds(self):
+        # 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
+
+        DASolver = self.DASolver
+
+        if self.DVGeo is None:
+            raise RuntimeError("calcFFD2XvSeeds is call but no DVGeo object found! Call add_dvgeo in the run script!")
+
+        dvName = DASolver.getOption("useAD")["dvName"]
+        seedIndex = DASolver.getOption("useAD")["seedIndex"]
+        # create xDVDot vec and initialize it with zeros
+        xDV = self.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(DASolver.xs0)
+        xSDot0 = DASolver.mapVector(xSDot0, DASolver.allSurfacesGroup, DASolver.designSurfacesGroup)
+
+        # get xSDot
+        xSDot = self.DVGeo.totalSensitivityProd(xDvDot, ptSetName="x_%s0" % self.discipline).reshape(xSDot0.shape)
+        # get xVDot
+        xVDot = DASolver.mesh.warpDerivFwd(xSDot)
+
+        return xVDot
+
     # calculate the residual
     def apply_nonlinear(self, inputs, outputs, residuals):
         DASolver = self.DASolver
@@ -756,14 +793,23 @@ def solve_nonlinear(self, inputs, outputs):
                     DASolver.setThermal(q_conduct)
                 else:
                     raise AnalysisError("discipline not valid!")
-            
+
             # before running the primal, we need to check if the mesh
             # quality is good
             meshOK = DASolver.solver.checkMesh()
 
             # solve the flow with the current design variable
             # if the mesh is not OK, do not run the primal
             if meshOK:
+                # if it is forward mode, we set the AD forward seeds before calling the primal
+                if DASolver.getOption("useAD")["mode"] == "forward":
+                    dvName = DASolver.getOption("useAD")["dvName"]
+                    dvType = DASolver.getOption("designVar")[dvName]["designVarType"]
+                    if dvType == "FFD":
+                        seeds = self.calcFFD2XvSeeds()
+                        DASolver.solverAD.setInputSeedForwardAD("volCoord", len(volCoords), volCoords, seeds)
+                    else:
+                        raise RuntimeError("dvType not supported for forwardAD")
                 DASolver()
             else:
                 DASolver.primalFail = 1

dafoam/pyDAFoam.py

@@ -976,45 +976,6 @@ def __call__(self):
         # update the mesh coordinates if DVGeo is set
         # add point set and update the mesh based on the DV values
 
-        if self.DVGeo is not None:
-
-            # if the point set is not in DVGeo add it first
-            if self.ptSetName not in self.DVGeo.points:
-
-                xs0 = self.mapVector(self.xs0, self.allSurfacesGroup, self.designSurfacesGroup)
-
-                self.DVGeo.addPointSet(xs0, self.ptSetName)
-                self.pointsSet = True
-
-            # set the surface coords xs
-            Info("DVGeo PointSet UpToDate: " + str(self.DVGeo.pointSetUpToDate(self.ptSetName)))
-            if not self.DVGeo.pointSetUpToDate(self.ptSetName):
-                Info("Updating DVGeo PointSet....")
-                xs = self.DVGeo.update(self.ptSetName, config=None)
-
-                # if we have surface geometry/mesh displacement computed by the structural solver,
-                # add the displace mesh here.
-                if self.surfGeoDisp is not None:
-                    xs += self.surfGeoDisp
-
-                self.setSurfaceCoordinates(xs, self.designSurfacesGroup)
-                Info("DVGeo PointSet UpToDate: " + str(self.DVGeo.pointSetUpToDate(self.ptSetName)))
-
-                # warp the mesh to get the new volume coordinates
-                Info("Warping the volume mesh....")
-                self.mesh.warpMesh()
-
-                xvNew = self.mesh.getSolverGrid()
-                self.xvFlatten2XvVec(xvNew, self.xvVec)
-
-            # if it is forward AD mode and we are computing the Xv derivatives
-            # call calcFFD2XvSeedVec
-            if self.getOption("useAD")["mode"] == "forward":
-                dvName = self.getOption("useAD")["dvName"]
-                dvType = self.getOption("designVar")[dvName]["designVarType"]
-                if dvType == "FFD":
-                    self.calcFFD2XvSeedVec()
-
         # now call the internal design var function to update DASolver parameters
         self.runInternalDVFunc()
 

src/adjoint/DASolver/DASolver.C

@@ -4799,12 +4799,47 @@ void DASolver::normalizeJacTVecProduct(
 #endif
 }
 
+void DASolver::setInputSeedForwardAD(
+    const word inputType,
+    const int inputSize,
+    const double* input,
+    const double* seed)
+{
+    /*
+    Description:
+        Set seeds for forward mode AD using the DAInput class
+    */
+#ifdef CODI_ADF
+    // initialize the input and output objects
+    autoPtr<DAInput> daInput(
+        DAInput::New(
+            inputType,
+            meshPtr_(),
+            daOptionPtr_(),
+            daModelPtr_(),
+            daIndexPtr_()));
+
+    scalarList inputList(inputSize, 0.0);
+
+    // assign the input array to the input list and the seed to its gradient()
+    // Note: we need to use scalarList for AD
+    forAll(inputList, idxI)
+    {
+        inputList[idxI] = input[idxI];
+        inputList[idxI].gradient() = seed[idxI];
+    }
+
+    // call daInput->run to assign inputList to OF variables
+    daInput->run(inputList);
+#endif
+}
+
 void DASolver::calcJacTVecProduct(
-    const word inputName,
+    const word inputType,
     const int inputSize,
     const int distributedInput,
     const double* input,
-    const word outputName,
+    const word outputType,
     const int outputSize,
     const int distributedOutput,
     const double* seed,
@@ -4816,30 +4851,38 @@ void DASolver::calcJacTVecProduct(
         Calculate the Jacobian-matrix-transposed and vector product for [dOutput/dInput]^T * psi
     
     Input:
-        inputName: name of the input
+        inputType: type of the input. This should be consistent with the child class names in DAInput
 
-        outputName: name of the output
+        inputSize: size of the input array
+
+        distributedInput: whether the input is distributed across processors in parallel. For example, the state variable is a distributed input, while the angle of attack is not a distributed input
 
         input: the actual value of the input array
 
-        psi: the vector for the mat-vec product
+        outputType: type of the output. This should be consistent with the child class names in DAOutput
+
+        outputSize: size of the output array
+
+        distributedInput: whether the output is distributed across processors in parallel. For example, residual is a distributed output, while the function (drag and lift) is not a distributed output
+
+        seed: the seed array
     
     Output:
-        product: the mat-vec product
+        product: the mat-vec product array
     */
 
     // initialize the input and output objects
     autoPtr<DAInput> daInput(
         DAInput::New(
-            inputName,
+            inputType,
             meshPtr_(),
             daOptionPtr_(),
             daModelPtr_(),
             daIndexPtr_()));
 
     autoPtr<DAOutput> daOutput(
         DAOutput::New(
-            outputName,
+            outputType,
             meshPtr_(),
             daOptionPtr_(),
             daModelPtr_(),
@@ -4915,8 +4958,8 @@ void DASolver::calcJacTVecProduct(
         }
     }
 
-    // we need to normalize the jacobian vector product if inputName == stateVar
-    this->normalizeJacTVecProduct(inputName, product);
+    // we need to normalize the jacobian vector product if inputType == stateVar
+    this->normalizeJacTVecProduct(inputType, product);
 
     // need to clear adjoint and tape after the computation is done!
     this->globalADTape_.clearAdjoints();

src/adjoint/DASolver/DASolver.H

@@ -431,15 +431,22 @@ public:
 
     /// calculate the Jacobian-matrix-transposed and vector product for product = [dOutput/dInput]^T * seed
     void calcJacTVecProduct(
-        const word inputName,
+        const word inputType,
         const int inputSize,
         const int distributedInput,
         const double* input,
-        const word outputName,
+        const word outputType,
         const int outputSize,
         const int distributedOutput,
         const double* seed,
         double* product);
+
+    void setInputSeedForwardAD(
+        const word inputType,
+        const int inputSize,
+        const double* input,
+        const double* seed
+    );
 
     /// compute dRdBCAD
     void calcdRdBCTPsiAD(

src/pyDASolvers/DASolvers.H

@@ -68,28 +68,42 @@ public:
 
     /// calculate the Jacobian-matrix-transposed and vector product for product = [dOutput/dInput]^T * seed
     void calcJacTVecProduct(
-        const word inputName,
+        const word inputType,
         const int inputSize,
         const int distributedInput,
         const double* input,
-        const word outputName,
+        const word outputType,
         const int outputSize,
         const int distributedOutput,
         const double* seed,
         double* product)
     {
         DASolverPtr_->calcJacTVecProduct(
-            inputName,
+            inputType,
             inputSize,
             distributedInput,
             input,
-            outputName,
+            outputType,
             outputSize,
             distributedOutput,
             seed,
             product);
     }
 
+    void setInputSeedForwardAD(
+        const word inputType,
+        const int inputSize,
+        const double* input,
+        const double* seed
+    )
+    {
+        DASolverPtr_->setInputSeedForwardAD(
+            inputType,
+            inputSize,
+            input,
+            seed);
+    }
+
     /// compute dRdWT
     void calcdRdWT(
         const label isPC,

src/pyDASolvers/pyDASolvers.pyx

@@ -48,6 +48,7 @@ cdef extern from "DASolvers.H" namespace "Foam":
         void initSolver()
         int solvePrimal()
         void calcJacTVecProduct(char *, int, int, double *, char *, int, int, double *, double *)
+        void setInputSeedForwardAD(char *, int, double *, double *)
         void calcdRdWT(int, PetscMat)
         void calcdRdWTPsiAD(PetscVec, PetscVec, PetscVec, PetscVec)
         void initializedRdWTMatrixFree()
@@ -200,12 +201,30 @@ cdef class pyDASolvers:
     def solvePrimal(self):
         return self._thisptr.solvePrimal()
 
+    def setInputSeedForwardAD(self,
+            inputType,
+            inputSize,
+            np.ndarray[double, ndim=1, mode="c"] inputs,
+            np.ndarray[double, ndim=1, mode="c"] seeds):
+
+        assert len(inputs) == inputSize, "invalid input array size!"
+        assert len(seeds) == inputSize, "invalid seed array size!"
+
+        cdef double *inputs_data = <double*>inputs.data
+        cdef double *seeds_data = <double*>seeds.data
+
+        self._thisptr.setInputSeedForwardAD(
+            inputType.encode(),
+            inputSize,
+            inputs_data,
+            seeds_data)
+
     def calcJacTVecProduct(self,
-            inputName,
+            inputType,
             inputSize,
             distributedInput,
             np.ndarray[double, ndim=1, mode="c"] inputs,
-            outputName,
+            outputType,
             outputSize,
             distributedOutput,
             np.ndarray[double, ndim=1, mode="c"] seeds,
@@ -220,11 +239,11 @@ cdef class pyDASolvers:
         cdef double *product_data = <double*>product.data
 
         self._thisptr.calcJacTVecProduct(
-            inputName.encode(),
+            inputType.encode(),
             inputSize,
             distributedInput,
             inputs_data,
-            outputName.encode(),
+            outputType.encode(),
             outputSize,
             distributedOutput,
             seeds_data,