Time accurate adjoint (#152)

* Changed the objFuncHist filename to objFuncTimeSeries to avoid confusion.

* Changed the hybridAdjoint flag such that it will not save states when inactive.

* Added an option to save intermediate variables for time-accurate adjoint.

* Added the dt residual terms and also fixed the missing ddtCorr terms for pisoFoam.

* Checked if both hybrid and time-accurate adjoints are active.

* Updated pisoFoam tests

* Added mesh.moving(false) right after mesh.movePoints() calls.

* Updated tests.

* Assigned fields and their oldTime(). Adjoint residual matched flow.

* Removed the time and timeIndex lists.

* Reverted back the time list.

* The time accurate residual worked for JacobianFree.

* Fixed the issue that the initial residual (ddt part) in the AD code is not correctly computed.

* Assigned all oldTimes even a variable does not need them.

* Added the calculation of dRdWOldTPsi using AD.

* Implemented the off-diagonal block terms. Also fixed a problem for first two step residuals. The last step residual is still not quite right.

* Updated the residual for pisoFoam with ts.

* Moved the unsteady adjoint functions from PisoFoam to DASolver.

* Updated tests.

* Added DAPimpleFoam primal.

* Added residual computation for pimpleFoam.

* Added the initOldTimes function to fix the issue that the first setTimeInstance call could not set the correct oldTime field.

* Updated the unsteadyAdjoint option.

* Reverted the adjoint calling sequence for hybridAdjoint.

* Updated test.

* Fixed an issue in assigning dRdWOldTPsi. Added zero timeAccurateVecs function.

* Fixed an issue for computing mean field.

* Fixed an issue for oldTime derivative by not registering oldTime fields if we do not need them in residual computation. Not doing this will get completely off adjoint derivative, especially for the calcdRdWTOldPsi function.

* Added the DALaplacianFoam solver.

* Removed all ResPartDeriv variables because they are not used.

* Fixed the missing part in solid MakeFiles.

* Added DAScalarTransportFoam.

* Fixed bugs for Laplacian and ScalarTransportFoam. Derivs not accurate.

* Fixed a bug in the actuatorDisk and actuatorLine models that impacted their circ force computation.

* Fixed a bug in the fCirc calculation for actuators to make the profile more smooth.

* Fixed a minor issue.

* Updated the tests.

* Relaxed the test tolerance.

* Changed the version to 226

dafoam/optFuncs.py

@@ -136,10 +136,10 @@ def calcObjFuncValuesMP(xDV):
     return funcsMP, fail
 
 
-def calcObjFuncValuesHybridAdjoint(xDV):
+def calcObjFuncValuesUnsteady(xDV):
     """
     Update the design surface and run the primal solver to get objective function values.
-    This is the hybrid adjoint version of calcObjFuncValues
+    This is the unsteady adjoint version of calcObjFuncValues
     """
 
     Info("\n")
@@ -164,9 +164,12 @@ def calcObjFuncValuesHybridAdjoint(xDV):
     # Solve the CFD problem
     DASolver()
 
-    # Set values for the hybrid adjoint objectives. This function needs to be
+    # Set values for the unsteady adjoint objectives. This function needs to be
     # implemented in run scripts
-    setHybridAdjointObjFuncs(DASolver, funcs, evalFuncs)
+    setObjFuncsUnsteady(DASolver, funcs, evalFuncs)
+
+    # assign state lists to mats
+    DASolver.setTimeInstanceVar(mode="list2Mat")
 
     b = time.time()
 
@@ -217,7 +220,7 @@ def calcObjFuncSens(xDV, funcs):
         Info("Objective Function Sensitivity: ")
         Info(funcsSens)
         Info("Adjoint Runtime: %g s" % (b - a))
-    
+
     # write the sensitivity values to file
     DASolver.writeTotalDeriv("totalDerivHist.txt", funcsSens, evalFuncs)
 
@@ -299,10 +302,10 @@ def calcObjFuncSensMP(xDV, funcs):
     return funcsSensMP, fail
 
 
-def calcObjFuncSensHybridAdjoint(xDV, funcs):
+def calcObjFuncSensUnsteady(xDV, funcs):
     """
     Run the adjoint solver and get objective function sensitivities.
-    This is the hybrid adjoint version of calcObjFuncSens
+    This is the unsteady adjoint version of calcObjFuncSens
     """
 
     Info("\n")
@@ -320,14 +323,32 @@ def calcObjFuncSensHybridAdjoint(xDV, funcs):
     # write the deform FFDs
     DASolver.writeDeformedFFDs()
 
+    # assign the state mats to lists
+    DASolver.setTimeInstanceVar(mode="mat2List")
+
     # Setup an empty dictionary for the evaluated derivative values
     funcsSensCombined = {}
 
     funcsSensAllInstances = []
 
-    nTimeInstances = DASolver.getOption("hybridAdjoint")["nTimeInstances"]
+    mode = DASolver.getOption("unsteadyAdjoint")["mode"]
+    nTimeInstances = DASolver.getOption("unsteadyAdjoint")["nTimeInstances"]
+    if mode == "hybridAdjoint":
+        iEnd = -1
+    elif mode == "timeAccurateAdjoint":
+        iEnd = 0
+
+    # NOTE: calling calcRes here is critical because it will setup the correct
+    # old time levels for setTimeInstanceField. Otherwise, the residual for the
+    # first adjoint time instance will be incorrect because the residuals have
+    # not been computed and the old time levels will be zeros for all variables,
+    # this will create issues for the setTimeInstanceField call (nOldTimes)
+    DASolver.calcPrimalResidualStatistics("calc")
+
+    # set these vectors zeros
+    DASolver.zeroTimeAccurateAdjointVectors()
 
-    for i in range(nTimeInstances):
+    for i in range(nTimeInstances - 1, iEnd, -1):
 
         Info("--Solving Adjoint for Time Instance %d--" % i)
 
@@ -355,13 +376,13 @@ def calcObjFuncSensHybridAdjoint(xDV, funcs):
 
         funcsSensAllInstances.append(funcsSens)
 
-    setHybridAdjointObjFuncsSens(DASolver, funcs, funcsSensAllInstances, funcsSensCombined)
+    setObjFuncsSensUnsteady(DASolver, funcs, funcsSensAllInstances, funcsSensCombined)
 
     funcsSensCombined["fail"] = fail
 
     # Print the current solution to the screen
     with np.printoptions(precision=16, threshold=5, suppress=True):
-        Info("Objective Function Sensitivity Hybrid Adjoint: ")
+        Info("Objective Function Sensitivity Unsteady Adjoint: ")
         Info(funcsSensCombined)
 
     b = time.time()

dafoam/pyDAFoam.py

@@ -11,7 +11,7 @@
 
 """
 
-__version__ = "2.2.5"
+__version__ = "2.2.6"
 
 import subprocess
 import os
@@ -57,6 +57,7 @@ class DAOPTION(object):
     ## - DASimpleFoam:            Incompressible steady-state flow solver for Navier-Stokes equations
     ## - DASimpleTFoam:           Incompressible steady-state flow solver for Navier-Stokes equations with temperature
     ## - DAPisoFoam:              Incompressible transient flow solver for Navier-Stokes equations
+    ## - DAPimpleFoam:            Incompressible transient flow solver for Navier-Stokes equations
     ## - DARhoSimpleFoam:         Compressible steady-state flow solver for Navier-Stokes equations (subsonic)
     ## - DARhoSimpleCFoam:        Compressible steady-state flow solver for Navier-Stokes equations (transonic)
     ## - DATurboFoam:             Compressible steady-state flow solver for Navier-Stokes equations (turbomachinery)
@@ -228,6 +229,20 @@ class DAOPTION(object):
     ##                "addToAdjoint": True,
     ##            }
     ##        },
+    ##        "THRUST": {
+    ##            "part1": {
+    ##                "type": "variableVolSum",
+    ##                "source": "boxToCell",
+    ##                "min": [-50.0, -50.0, -50.0],
+    ##                "max": [50.0, 50.0, 50.0],
+    ##                "varName": "fvSource",
+    ##                "varType": "vector",
+    ##                "component": 0,
+    ##                "isSquare": 0,
+    ##                "scale": 1.0,
+    ##                "addToAdjoint": True,
+    ##            },
+    ##        },
     ##        "FI": {
     ##            "part1": {
     ##                "type": "stateErrorNorm",
@@ -397,9 +412,10 @@ class DAOPTION(object):
     ## This is used only for transonic solvers such as DARhoSimpleCFoam
     transonicPCOption = -1
 
-    ## Options for hybrid adjoint. Here nTimeInstances is the number of time instances
-    ## periodicity is the periodicity of flow oscillation
-    hybridAdjoint = {"active": False, "nTimeInstances": -1, "periodicity": -1.0}
+    ## Options for unsteady adjoint. mode can be hybridAdjoint or timeAccurateAdjoint
+    ## Here nTimeInstances is the number of time instances and periodicity is the 
+    ## periodicity of flow oscillation (hybrid adjoint only)
+    unsteadyAdjoint = {"mode": "None", "nTimeInstances": -1, "periodicity": -1.0}
 
     ## At which iteration should we start the averaging of objective functions.
     ## This is only used for unsteady solvers
@@ -703,6 +719,12 @@ def __init__(self, comm=None, options=None):
         # initialize the adjoint vector dict
         self.adjVectors = self._initializeAdjVectors()
 
+        # initialize the dRdWOldTPsi vectors
+        self._initializeTimeAccurateAdjointVectors()
+
+        # check if the combination of options is valid.
+        self._checkOptions()
+
         Info("pyDAFoam initialization done!")
 
         return
@@ -714,9 +736,9 @@ def _solverRegistry(self):
         """
 
         self.solverRegistry = {
-            "Incompressible": ["DASimpleFoam", "DASimpleTFoam", "DAPisoFoam"],
+            "Incompressible": ["DASimpleFoam", "DASimpleTFoam", "DAPisoFoam", "DAPimpleFoam"],
             "Compressible": ["DARhoSimpleFoam", "DARhoSimpleCFoam", "DATurboFoam"],
-            "Solid": ["DASolidDisplacementFoam"],
+            "Solid": ["DASolidDisplacementFoam", "DALaplacianFoam", "DAScalarTransportFoam"],
         }
 
     def __call__(self):
@@ -831,6 +853,62 @@ def _initializeAdjTotalDeriv(self):
 
         return adjTotalDeriv
 
+    def _initializeTimeAccurateAdjointVectors(self):
+        """
+        Initialize the dRdWTPsi vectors for time accurate adjoint.
+        Here we need to initialize current time step and two previous
+        time steps 0 and 00 for both state and residuals. This is
+        because the backward ddt scheme depends on U, U0, and U00
+        """
+        if self.getOption("unsteadyAdjoint")["mode"] == "timeAccurateAdjoint":
+            objFuncDict = self.getOption("objFunc")
+            wSize = self.solver.getNLocalAdjointStates()
+            self.dRdW0TPsi = {}
+            self.dRdW00TPsi = {}
+            self.dR0dW0TPsi = {}
+            self.dR0dW00TPsi = {}
+            self.dR00dW0TPsi = {}
+            self.dR00dW00TPsi = {}
+            for objFuncName in objFuncDict:
+                if objFuncName in self.objFuncNames4Adj:
+                    vecA = PETSc.Vec().create(PETSc.COMM_WORLD)
+                    vecA.setSizes((wSize, PETSc.DECIDE), bsize=1)
+                    vecA.setFromOptions()
+                    vecA.zeroEntries()
+                    self.dRdW0TPsi[objFuncName] = vecA
+
+                    vecB = vecA.duplicate()
+                    vecB.zeroEntries()
+                    self.dRdW00TPsi[objFuncName] = vecB
+
+                    vecC = vecA.duplicate()
+                    vecC.zeroEntries()
+                    self.dR0dW0TPsi[objFuncName] = vecC
+
+                    vecD = vecA.duplicate()
+                    vecD.zeroEntries()
+                    self.dR0dW00TPsi[objFuncName] = vecD
+
+                    vecE = vecA.duplicate()
+                    vecE.zeroEntries()
+                    self.dR00dW0TPsi[objFuncName] = vecE
+
+                    vecF = vecA.duplicate()
+                    vecF.zeroEntries()
+                    self.dR00dW00TPsi[objFuncName] = vecF
+
+    def zeroTimeAccurateAdjointVectors(self):
+        if self.getOption("unsteadyAdjoint")["mode"] == "timeAccurateAdjoint":
+            objFuncDict = self.getOption("objFunc")
+            for objFuncName in objFuncDict:
+                if objFuncName in self.objFuncNames4Adj:
+                    self.dRdW0TPsi[objFuncName].zeroEntries()
+                    self.dRdW00TPsi[objFuncName].zeroEntries()
+                    self.dR0dW0TPsi[objFuncName].zeroEntries()
+                    self.dR0dW00TPsi[objFuncName].zeroEntries()
+                    self.dR00dW0TPsi[objFuncName].zeroEntries()
+                    self.dR00dW00TPsi[objFuncName].zeroEntries()
+
     def _calcObjFuncNames4Adj(self):
         """
         Compute the objective function names for which we solve the adjoint equation
@@ -857,6 +935,20 @@ def _calcObjFuncNames4Adj(self):
                     raise Error("addToAdjoint can be either True or False")
         return objFuncList
 
+    def _checkOptions(self):
+        """
+        Check if the combination of options are valid.
+        For example, if timeAccurateAdjoint is active, we have to set
+        adjJacobianOption = JacobianFree
+        NOTE: we should add all possible checks here!
+        """
+        # check time accurate adjoint
+        if self.getOption("unsteadyAdjoint")["mode"] == "timeAccurateAdjoint":
+            if not self.getOption("adjJacobianOption") == "JacobianFree":
+                raise Error("timeAccurateAdjoint only supports adjJacobianOption=JacobianFree!")
+
+        # check other combinations...
+
     def saveMultiPointField(self, indexMP):
         """
         Save the state variable vector to self.wVecMPList
@@ -885,18 +977,76 @@ def setMultiPointField(self, indexMP):
 
         return
 
+    def calcPrimalResidualStatistics(self, mode):
+        if self.getOption("adjJacobianOption") == "JacobianFD":
+            self.solver.calcPrimalResidualStatistics(mode.encode())
+        elif self.getOption("adjJacobianOption") == "JacobianFree":
+            self.solverAD.calcPrimalResidualStatistics(mode.encode())
+
     def setTimeInstanceField(self, instanceI):
         """
         Set the OpenFOAM state variables based on instance index
         """
 
-        self.solver.setTimeInstanceField(instanceI)
+        if self.getOption("adjJacobianOption") == "JacobianFD":
+            solver = self.solver
+        elif self.getOption("adjJacobianOption") == "JacobianFree":
+            solver = self.solverAD
+
+        solver.setTimeInstanceField(instanceI)
         # NOTE: we need to set the OF field to wVec vector here!
         # this is because we will assign self.wVec to the solveAdjoint function later
-        self.solver.ofField2StateVec(self.wVec)
+        solver.ofField2StateVec(self.wVec)
 
         return
 
+    def initTimeInstanceMats(self):
+
+        nLocalAdjointStates = self.solver.getNLocalAdjointStates()
+        nLocalAdjointBoundaryStates = self.solver.getNLocalAdjointBoundaryStates()
+        nTimeInstances = -99999
+        adjMode = self.getOption("unsteadyAdjoint")["mode"]
+        if adjMode == "hybridAdjoint" or adjMode == "timeAccurateAdjoint":
+            nTimeInstances = self.getOption("unsteadyAdjoint")["nTimeInstances"]
+
+        self.stateMat = PETSc.Mat().create(PETSc.COMM_WORLD)
+        self.stateMat.setSizes(((nLocalAdjointStates, None), (None, nTimeInstances)))
+        self.stateMat.setFromOptions()
+        self.stateMat.setPreallocationNNZ((nTimeInstances, nTimeInstances))
+        self.stateMat.setUp()
+
+        self.stateBCMat = PETSc.Mat().create(PETSc.COMM_WORLD)
+        self.stateBCMat.setSizes(((nLocalAdjointBoundaryStates, None), (None, nTimeInstances)))
+        self.stateBCMat.setFromOptions()
+        self.stateBCMat.setPreallocationNNZ((nTimeInstances, nTimeInstances))
+        self.stateBCMat.setUp()
+
+        self.timeVec = PETSc.Vec().createSeq(nTimeInstances, bsize=1, comm=PETSc.COMM_SELF)
+        self.timeIdxVec = PETSc.Vec().createSeq(nTimeInstances, bsize=1, comm=PETSc.COMM_SELF)
+
+    def initOldTimes(self):
+        # No need to initialize oldTimes for FD
+        if self.getOption("adjJacobianOption") == "JacobianFree":
+            self.solverAD.initOldTimes()
+
+    def setTimeInstanceVar(self, mode):
+
+        if mode == "list2Mat":
+            self.solver.setTimeInstanceVar(mode.encode(), self.stateMat, self.stateBCMat, self.timeVec, self.timeIdxVec)
+        elif mode == "mat2List":
+            if self.getOption("adjJacobianOption") == "JacobianFD":
+                self.solver.setTimeInstanceVar(
+                    mode.encode(), self.stateMat, self.stateBCMat, self.timeVec, self.timeIdxVec
+                )
+            elif self.getOption("adjJacobianOption") == "JacobianFree":
+                self.solverAD.setTimeInstanceVar(
+                    mode.encode(), self.stateMat, self.stateBCMat, self.timeVec, self.timeIdxVec
+                )
+            else:
+                raise Error("adjJacobianOption can only be either JacobianFD or JacobianFree!")
+        else:
+            raise Error("mode can only be either mat2List or list2Mat!")
+
     def writeDesignVariable(self, fileName, xDV):
         """
         Write the design variable history to files in the json format
@@ -1703,12 +1853,26 @@ def solveAdjoint(self):
                 elif self.getOption("adjJacobianOption") == "JacobianFree":
                     self.solverAD.calcdFdWAD(self.xvVec, self.wVec, objFuncName.encode(), dFdW)
 
+                # if it is time accurate adjoint, add extra terms for dFdW
+                if self.getOption("unsteadyAdjoint")["mode"] == "timeAccurateAdjoint":
+                    # first copy the vectors from previous residual time step level
+                    self.dR0dW0TPsi[objFuncName].copy(self.dR00dW0TPsi[objFuncName])
+                    self.dR0dW00TPsi[objFuncName].copy(self.dR00dW00TPsi[objFuncName])
+                    self.dRdW0TPsi[objFuncName].copy(self.dR0dW0TPsi[objFuncName])
+                    self.dRdW00TPsi[objFuncName].copy(self.dR0dW00TPsi[objFuncName])
+                    dFdW.axpy(-1.0, self.dR0dW0TPsi[objFuncName])
+                    dFdW.axpy(-1.0, self.dR00dW00TPsi[objFuncName])
+
                 # Initialize the adjoint vector psi and solve for it
                 if self.getOption("adjJacobianOption") == "JacobianFD":
                     self.adjointFail = self.solver.solveLinearEqn(ksp, dFdW, self.adjVectors[objFuncName])
                 elif self.getOption("adjJacobianOption") == "JacobianFree":
                     self.adjointFail = self.solverAD.solveLinearEqn(ksp, dFdW, self.adjVectors[objFuncName])
 
+                if self.getOption("unsteadyAdjoint")["mode"] == "timeAccurateAdjoint":
+                    self.solverAD.calcdRdWOldTPsiAD(1, self.adjVectors[objFuncName], self.dRdW0TPsi[objFuncName])
+                    self.solverAD.calcdRdWOldTPsiAD(2, self.adjVectors[objFuncName], self.dRdW00TPsi[objFuncName])
+
                 dFdW.destroy()
 
         ksp.destroy()
@@ -2206,6 +2370,10 @@ def _initSolver(self):
         if self.getOption("printDAOptions"):
             self.solver.printAllOptions()
 
+        adjMode = self.getOption("unsteadyAdjoint")["mode"]
+        if adjMode == "hybridAdjoint" or adjMode == "timeAccurateAdjoint":
+            self.initTimeInstanceMats()
+
         self.solverInitialized = 1
 
         return