From e8ac6479a2c8032bd52f6feab7c8ae1ca12d9bac Mon Sep 17 00:00:00 2001 From: Ping He Date: Fri, 20 Dec 2024 22:00:02 -0600 Subject: [PATCH] Enabled forward mode AD. Verified machine-precision accurate adjoint for DASimpleFoam --- dafoam/mphys/mphys_dafoam.py | 48 ++++++++++++++++++++++++- dafoam/pyDAFoam.py | 39 -------------------- src/adjoint/DASolver/DASolver.C | 63 +++++++++++++++++++++++++++------ src/adjoint/DASolver/DASolver.H | 11 ++++-- src/pyDASolvers/DASolvers.H | 22 +++++++++--- src/pyDASolvers/pyDASolvers.pyx | 27 +++++++++++--- 6 files changed, 150 insertions(+), 60 deletions(-) diff --git a/dafoam/mphys/mphys_dafoam.py b/dafoam/mphys/mphys_dafoam.py index f6306543..555a43b0 100644 --- a/dafoam/mphys/mphys_dafoam.py +++ b/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,7 +793,7 @@ 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() @@ -764,6 +801,15 @@ def solve_nonlinear(self, inputs, outputs): # 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 diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py index f3d3ff99..431b6d10 100755 --- a/dafoam/pyDAFoam.py +++ b/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() diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C index 4bbfbee1..1f134b9f 100644 --- a/src/adjoint/DASolver/DASolver.C +++ b/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::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,22 +4851,30 @@ 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::New( - inputName, + inputType, meshPtr_(), daOptionPtr_(), daModelPtr_(), @@ -4839,7 +4882,7 @@ void DASolver::calcJacTVecProduct( autoPtr 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(); diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H index 3b7fa0da..f306947a 100644 --- a/src/adjoint/DASolver/DASolver.H +++ b/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( diff --git a/src/pyDASolvers/DASolvers.H b/src/pyDASolvers/DASolvers.H index be9c340c..7050f57c 100755 --- a/src/pyDASolvers/DASolvers.H +++ b/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, diff --git a/src/pyDASolvers/pyDASolvers.pyx b/src/pyDASolvers/pyDASolvers.pyx index 2d5e1194..bc15c0cd 100755 --- a/src/pyDASolvers/pyDASolvers.pyx +++ b/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 = inputs.data + cdef double *seeds_data = 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 = product.data self._thisptr.calcJacTVecProduct( - inputName.encode(), + inputType.encode(), inputSize, distributedInput, inputs_data, - outputName.encode(), + outputType.encode(), outputSize, distributedOutput, seeds_data,