From 092d8a3726fd07adbe20a67026cf5aba858c1266 Mon Sep 17 00:00:00 2001 From: Ping He Date: Sat, 19 Jun 2021 14:46:05 -0400 Subject: [PATCH] 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 | 43 +- dafoam/pyDAFoam.py | 184 ++++- src/adjoint/DAField/DAField.C | 222 ++++-- src/adjoint/DAField/DAField.H | 3 +- .../DAFvSource/DAFvSourceActuatorDisk.C | 31 +- .../DAFvSource/DAFvSourceActuatorLine.C | 23 +- src/adjoint/DAIndex/DAIndex.C | 4 + src/adjoint/DAIntmdVar/DAIntmdVar.C | 18 +- .../DATurbulenceModel/DASpalartAllmaras.C | 8 - .../DATurbulenceModel/DASpalartAllmaras.H | 1 - .../DATurbulenceModel/DASpalartAllmarasFv3.C | 8 - .../DATurbulenceModel/DASpalartAllmarasFv3.H | 1 - .../DASpalartAllmarasFv3Beta.C | 8 - .../DASpalartAllmarasFv3Beta.H | 1 - .../DAModel/DATurbulenceModel/DAkOmega.C | 18 +- .../DAModel/DATurbulenceModel/DAkOmega.H | 2 - .../DAModel/DATurbulenceModel/DAkOmegaSST.C | 16 - .../DAModel/DATurbulenceModel/DAkOmegaSST.H | 2 - .../DAModel/DATurbulenceModel/DAkOmegaSSTLM.C | 32 - .../DAModel/DATurbulenceModel/DAkOmegaSSTLM.H | 4 - .../DAResidual/DAResidualLaplacianFoam.C | 102 +++ .../DAResidual/DAResidualLaplacianFoam.H | 79 +++ src/adjoint/DAResidual/DAResidualPimpleFoam.C | 194 ++++++ src/adjoint/DAResidual/DAResidualPimpleFoam.H | 99 +++ src/adjoint/DAResidual/DAResidualPisoFoam.C | 26 +- src/adjoint/DAResidual/DAResidualPisoFoam.H | 6 +- .../DAResidual/DAResidualRhoSimpleCFoam.C | 4 - .../DAResidual/DAResidualRhoSimpleCFoam.H | 4 - .../DAResidual/DAResidualRhoSimpleFoam.C | 4 - .../DAResidual/DAResidualRhoSimpleFoam.H | 4 - .../DAResidualScalarTransportFoam.C | 107 +++ .../DAResidualScalarTransportFoam.H | 80 +++ src/adjoint/DAResidual/DAResidualSimpleFoam.C | 3 - src/adjoint/DAResidual/DAResidualSimpleFoam.H | 3 - .../DAResidual/DAResidualSimpleTFoam.C | 4 - .../DAResidual/DAResidualSimpleTFoam.H | 4 - .../DAResidualSolidDisplacementFoam.C | 1 - .../DAResidualSolidDisplacementFoam.H | 1 - src/adjoint/DAResidual/DAResidualTurboFoam.C | 4 - src/adjoint/DAResidual/DAResidualTurboFoam.H | 4 - .../DALaplacianFoam/DALaplacianFoam.C | 229 ++++++ .../DALaplacianFoam/DALaplacianFoam.H | 76 ++ .../DALaplacianFoam/createFieldsLaplacian.H | 31 + .../DALaplacianFoam/createRefsLaplacian.H | 11 + .../DASolver/DAPimpleFoam/DAPimpleFoam.C | 264 +++++++ .../DASolver/DAPimpleFoam/DAPimpleFoam.H | 114 +++ .../DASolver/DAPimpleFoam/UEqnPimple.H | 26 + .../DAPimpleFoam/createFieldsPimple.H | 52 ++ .../DASolver/DAPimpleFoam/createRefsPimple.H | 18 + .../DASolver/DAPimpleFoam/pEqnPimple.H | 55 ++ src/adjoint/DASolver/DAPisoFoam/DAPisoFoam.C | 182 ++--- src/adjoint/DASolver/DAPisoFoam/DAPisoFoam.H | 35 +- .../DASolver/DAPisoFoam/createRefsPiso.H | 2 +- src/adjoint/DASolver/DAPisoFoam/pEqnPiso.H | 5 +- .../DARhoSimpleCFoam/createRefsRhoSimpleC.H | 2 +- .../DARhoSimpleFoam/createRefsRhoSimple.H | 2 +- .../DAScalarTransportFoam.C | 234 +++++++ .../DAScalarTransportFoam.H | 82 +++ .../createFieldsScalarTransport.H | 58 ++ .../createRefsScalarTransport.H | 12 + .../DASolver/DASimpleFoam/createRefsSimple.H | 2 +- .../DASimpleTFoam/createRefsSimpleT.H | 2 +- .../createRefsSolidDisplacement.H | 2 +- src/adjoint/DASolver/DASolver.C | 653 ++++++++++++++++-- src/adjoint/DASolver/DASolver.H | 60 +- .../DASolver/DATurboFoam/createRefsTurbo.H | 2 +- .../DAStateInfo/DAStateInfoLaplacianFoam.C | 81 +++ .../DAStateInfo/DAStateInfoLaplacianFoam.H | 55 ++ .../DAStateInfo/DAStateInfoPimpleFoam.C | 116 ++++ .../DAStateInfo/DAStateInfoPimpleFoam.H | 55 ++ .../DAStateInfoScalarTransportFoam.C | 81 +++ .../DAStateInfoScalarTransportFoam.H | 55 ++ src/adjoint/Make/files_Incompressible | 5 +- src/adjoint/Make/files_Solid | 9 + src/include/DAMacroFunctions.H | 35 +- src/include/createPimpleControlPython.H | 14 + src/pyDASolvers/DASolvers.H | 38 + src/pyDASolvers/pyDASolvers.pyx | 20 + tests/refs/DAFoam_Test_DAPisoFoamACTLRef.txt | 66 +- tests/refs/DAFoam_Test_DAPisoFoamRef.txt | 58 +- .../DAFoam_Test_DARhoSimpleCFoamADRef.txt | 18 +- .../refs/DAFoam_Test_DARhoSimpleCFoamRef.txt | 18 +- .../refs/DAFoam_Test_DARhoSimpleFoamADRef.txt | 152 ++-- .../DAFoam_Test_DARhoSimpleFoamMRFRef.txt | 18 +- tests/refs/DAFoam_Test_DARhoSimpleFoamRef.txt | 26 +- .../DAFoam_Test_DARhoSimpleFoamUBendRef.txt | 26 +- tests/refs/DAFoam_Test_DASimpleFoamADRef.txt | 18 +- tests/refs/DAFoam_Test_DASimpleFoamMRFRef.txt | 18 +- tests/refs/DAFoam_Test_DASimpleFoamRef.txt | 140 ++-- tests/refs/DAFoam_Test_DASimpleTFoamRef.txt | 26 +- ...DAFoam_Test_DASolidDisplacementFoamRef.txt | 26 +- .../DAFoam_Test_DATurboFoamSubsonicRef.txt | 34 +- .../DAFoam_Test_DATurboFoamTransonicRef.txt | 66 +- tests/refs/DAFoam_Test_IntegrationRef.txt | 42 +- tests/runTests_DAPisoFoam.py | 20 +- tests/runTests_DAPisoFoamACTL.py | 20 +- tests/runTests_DARhoSimpleCFoam.py | 2 +- tests/runTests_DARhoSimpleCFoamAD.py | 2 +- tests/runTests_DARhoSimpleFoam.py | 2 +- tests/runTests_DARhoSimpleFoamAD.py | 6 +- tests/runTests_DARhoSimpleFoamMRF.py | 2 +- tests/runTests_DARhoSimpleFoamUBend.py | 2 +- tests/runTests_DASimpleFoam.py | 4 +- tests/runTests_DASimpleFoamAD.py | 2 +- tests/runTests_DASimpleFoamMRF.py | 2 +- tests/runTests_DASimpleTFoam.py | 2 +- tests/runTests_DASolidDisplacementFoam.py | 2 +- tests/runTests_DATurboFoamSubsonic.py | 2 +- tests/runTests_DATurboFoamTransonic.py | 2 +- 109 files changed, 4052 insertions(+), 909 deletions(-) create mode 100755 src/adjoint/DAResidual/DAResidualLaplacianFoam.C create mode 100755 src/adjoint/DAResidual/DAResidualLaplacianFoam.H create mode 100755 src/adjoint/DAResidual/DAResidualPimpleFoam.C create mode 100755 src/adjoint/DAResidual/DAResidualPimpleFoam.H create mode 100755 src/adjoint/DAResidual/DAResidualScalarTransportFoam.C create mode 100755 src/adjoint/DAResidual/DAResidualScalarTransportFoam.H create mode 100755 src/adjoint/DASolver/DALaplacianFoam/DALaplacianFoam.C create mode 100755 src/adjoint/DASolver/DALaplacianFoam/DALaplacianFoam.H create mode 100755 src/adjoint/DASolver/DALaplacianFoam/createFieldsLaplacian.H create mode 100755 src/adjoint/DASolver/DALaplacianFoam/createRefsLaplacian.H create mode 100755 src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C create mode 100755 src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.H create mode 100755 src/adjoint/DASolver/DAPimpleFoam/UEqnPimple.H create mode 100755 src/adjoint/DASolver/DAPimpleFoam/createFieldsPimple.H create mode 100755 src/adjoint/DASolver/DAPimpleFoam/createRefsPimple.H create mode 100755 src/adjoint/DASolver/DAPimpleFoam/pEqnPimple.H create mode 100755 src/adjoint/DASolver/DAScalarTransportFoam/DAScalarTransportFoam.C create mode 100755 src/adjoint/DASolver/DAScalarTransportFoam/DAScalarTransportFoam.H create mode 100755 src/adjoint/DASolver/DAScalarTransportFoam/createFieldsScalarTransport.H create mode 100755 src/adjoint/DASolver/DAScalarTransportFoam/createRefsScalarTransport.H create mode 100755 src/adjoint/DAStateInfo/DAStateInfoLaplacianFoam.C create mode 100755 src/adjoint/DAStateInfo/DAStateInfoLaplacianFoam.H create mode 100755 src/adjoint/DAStateInfo/DAStateInfoPimpleFoam.C create mode 100755 src/adjoint/DAStateInfo/DAStateInfoPimpleFoam.H create mode 100755 src/adjoint/DAStateInfo/DAStateInfoScalarTransportFoam.C create mode 100755 src/adjoint/DAStateInfo/DAStateInfoScalarTransportFoam.H create mode 100755 src/include/createPimpleControlPython.H diff --git a/dafoam/optFuncs.py b/dafoam/optFuncs.py index ffd2088a..57c2da02 100755 --- a/dafoam/optFuncs.py +++ b/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() diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py index e8d40fcd..01efe1c1 100755 --- a/dafoam/pyDAFoam.py +++ b/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 diff --git a/src/adjoint/DAField/DAField.C b/src/adjoint/DAField/DAField.C index c15fbd5c..eb2d5f2d 100755 --- a/src/adjoint/DAField/DAField.C +++ b/src/adjoint/DAField/DAField.C @@ -256,6 +256,7 @@ void DAField::pointVec2OFMesh(const Vec xvVec) const // movePoints update the mesh metrics such as volume, surface area and cell centers fvMesh& mesh = const_cast(mesh_); mesh.movePoints(meshPoints); + mesh.moving(false); } void DAField::ofMesh2PointVec(Vec xvVec) const @@ -1161,7 +1162,8 @@ void DAField::ofField2List( void DAField::list2OFField( const scalarList& stateList, - const scalarList& stateBoundaryList) const + const scalarList& stateBoundaryList, + const label oldTimeLevel) const { /* Description: @@ -1170,6 +1172,8 @@ void DAField::list2OFField( Input: stateList: scalar list of states stateBoundaryList: scalar list of boundary states + oldTimeLevel: assign to oldTime field instead of the original field, this will + be used in time-accurate adjoint Output: OpenFoam field variables @@ -1188,30 +1192,66 @@ void DAField::list2OFField( label localBFaceI = 0; + if (oldTimeLevel > 2 || oldTimeLevel < 0) + { + FatalErrorIn("") << "oldTimeLevel not valid!" + << abort(FatalError); + } + forAll(stateInfo_["volVectorStates"], idxI) { // lookup state from meshDb makeState(stateInfo_["volVectorStates"][idxI], volVectorField, db); - forAll(mesh_.cells(), cellI) + label maxOldTimes = state.nOldTimes(); + + if (maxOldTimes >= oldTimeLevel) { - for (label comp = 0; comp < 3; comp++) + forAll(mesh_.cells(), cellI) { - label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp); - state[cellI][comp] = stateList[localIdx]; + for (label comp = 0; comp < 3; comp++) + { + label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI, comp); + if (oldTimeLevel == 0) + { + state[cellI][comp] = stateList[localIdx]; + } + else if (oldTimeLevel == 1) + { + state.oldTime()[cellI][comp] = stateList[localIdx]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime()[cellI][comp] = stateList[localIdx]; + } + } } - } - forAll(state.boundaryField(), patchI) - { - if (state.boundaryField()[patchI].size() > 0) + forAll(state.boundaryField(), patchI) { - forAll(state.boundaryField()[patchI], faceI) + if (state.boundaryField()[patchI].size() > 0) { - for (label comp = 0; comp < 3; comp++) + forAll(state.boundaryField()[patchI], faceI) { - state.boundaryFieldRef()[patchI][faceI][comp] = stateBoundaryList[localBFaceI]; - localBFaceI++; + for (label comp = 0; comp < 3; comp++) + { + if (oldTimeLevel == 0) + { + state.boundaryFieldRef()[patchI][faceI][comp] = + stateBoundaryList[localBFaceI]; + } + else if (oldTimeLevel == 1) + { + state.oldTime().boundaryFieldRef()[patchI][faceI][comp] = + stateBoundaryList[localBFaceI]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime().boundaryFieldRef()[patchI][faceI][comp] = + stateBoundaryList[localBFaceI]; + } + localBFaceI++; + } } } } @@ -1223,20 +1263,51 @@ void DAField::list2OFField( // lookup state from meshDb makeState(stateInfo_["volScalarStates"][idxI], volScalarField, db); - forAll(mesh_.cells(), cellI) - { - label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI); - state[cellI] = stateList[localIdx]; - } + label maxOldTimes = state.nOldTimes(); - forAll(state.boundaryField(), patchI) + if (maxOldTimes >= oldTimeLevel) { - if (state.boundaryField()[patchI].size() > 0) + + forAll(mesh_.cells(), cellI) { - forAll(state.boundaryField()[patchI], faceI) + label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI); + if (oldTimeLevel == 0) { - state.boundaryFieldRef()[patchI][faceI] = stateBoundaryList[localBFaceI]; - localBFaceI++; + state[cellI] = stateList[localIdx]; + } + else if (oldTimeLevel == 1) + { + state.oldTime()[cellI] = stateList[localIdx]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime()[cellI] = stateList[localIdx]; + } + } + + forAll(state.boundaryField(), patchI) + { + if (state.boundaryField()[patchI].size() > 0) + { + forAll(state.boundaryField()[patchI], faceI) + { + if (oldTimeLevel == 0) + { + state.boundaryFieldRef()[patchI][faceI] = + stateBoundaryList[localBFaceI]; + } + else if (oldTimeLevel == 1) + { + state.oldTime().boundaryFieldRef()[patchI][faceI] = + stateBoundaryList[localBFaceI]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime().boundaryFieldRef()[patchI][faceI] = + stateBoundaryList[localBFaceI]; + } + localBFaceI++; + } } } } @@ -1247,20 +1318,51 @@ void DAField::list2OFField( // lookup state from meshDb makeState(stateInfo_["modelStates"][idxI], volScalarField, db); - forAll(mesh_.cells(), cellI) - { - label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI); - state[cellI] = stateList[localIdx]; - } + label maxOldTimes = state.nOldTimes(); - forAll(state.boundaryField(), patchI) + if (maxOldTimes >= oldTimeLevel) { - if (state.boundaryField()[patchI].size() > 0) + + forAll(mesh_.cells(), cellI) { - forAll(state.boundaryField()[patchI], faceI) + label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, cellI); + if (oldTimeLevel == 0) { - state.boundaryFieldRef()[patchI][faceI] = stateBoundaryList[localBFaceI]; - localBFaceI++; + state[cellI] = stateList[localIdx]; + } + else if (oldTimeLevel == 1) + { + state.oldTime()[cellI] = stateList[localIdx]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime()[cellI] = stateList[localIdx]; + } + } + + forAll(state.boundaryField(), patchI) + { + if (state.boundaryField()[patchI].size() > 0) + { + forAll(state.boundaryField()[patchI], faceI) + { + if (oldTimeLevel == 0) + { + state.boundaryFieldRef()[patchI][faceI] = + stateBoundaryList[localBFaceI]; + } + else if (oldTimeLevel == 1) + { + state.oldTime().boundaryFieldRef()[patchI][faceI] = + stateBoundaryList[localBFaceI]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime().boundaryFieldRef()[patchI][faceI] = + stateBoundaryList[localBFaceI]; + } + localBFaceI++; + } } } } @@ -1271,19 +1373,51 @@ void DAField::list2OFField( // lookup state from meshDb makeState(stateInfo_["surfaceScalarStates"][idxI], surfaceScalarField, db); - forAll(mesh_.faces(), faceI) + label maxOldTimes = state.nOldTimes(); + + if (maxOldTimes >= oldTimeLevel) { - label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI); - if (faceI < daIndex_.nLocalInternalFaces) - { - state[faceI] = stateList[localIdx]; - } - else + + forAll(mesh_.faces(), faceI) { - label relIdx = faceI - daIndex_.nLocalInternalFaces; - const label& patchIdx = daIndex_.bFacePatchI[relIdx]; - const label& faceIdx = daIndex_.bFaceFaceI[relIdx]; - state.boundaryFieldRef()[patchIdx][faceIdx] = stateList[localIdx]; + label localIdx = daIndex_.getLocalAdjointStateIndex(stateName, faceI); + if (faceI < daIndex_.nLocalInternalFaces) + { + + if (oldTimeLevel == 0) + { + state[faceI] = stateList[localIdx]; + } + else if (oldTimeLevel == 1) + { + state.oldTime()[faceI] = stateList[localIdx]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime()[faceI] = stateList[localIdx]; + } + } + else + { + label relIdx = faceI - daIndex_.nLocalInternalFaces; + const label& patchIdx = daIndex_.bFacePatchI[relIdx]; + const label& faceIdx = daIndex_.bFaceFaceI[relIdx]; + if (oldTimeLevel == 0) + { + state.boundaryFieldRef()[patchIdx][faceIdx] = + stateList[localIdx]; + } + else if (oldTimeLevel == 1) + { + state.oldTime().boundaryFieldRef()[patchIdx][faceIdx] = + stateList[localIdx]; + } + else if (oldTimeLevel == 2) + { + state.oldTime().oldTime().boundaryFieldRef()[patchIdx][faceIdx] = + stateList[localIdx]; + } + } } } } diff --git a/src/adjoint/DAField/DAField.H b/src/adjoint/DAField/DAField.H index cc14ba7f..659f7966 100755 --- a/src/adjoint/DAField/DAField.H +++ b/src/adjoint/DAField/DAField.H @@ -100,7 +100,8 @@ public: /// assign the fields in OpenFOAM based on the scalar list of states void list2OFField( const scalarList& stateList, - const scalarList& stateBounaryList) const; + const scalarList& stateBounaryList, + const label oldTimeLevel) const; /// set the boundary conditions based on parameters defined in DAOption void setPrimalBoundaryConditions(const label printInfo = 1); diff --git a/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C b/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C index 13f0a94d..7d1fa422 100755 --- a/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C +++ b/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C @@ -46,6 +46,9 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource) source = rStar * sqrt(1.0 - rStar) * scale where rStar is the normalized radial location and scale is used to make sure the integral force equals the desired total thrust + + NOTE: rotDir = right means propeller rotates clockwise viewed from + the tail of the aircraft looking forward There are two options to assign the source term: 1. cylinderAnnulusToCell. Users prescribe a cylinderAnnulus and the fvSource will be @@ -161,13 +164,13 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource) vector cellC2AVecC(vector::zero); if (rotDir == "left") { - // this assumes right hand rotation of propellers - cellC2AVecC = cellC2AVecR ^ cellC2AVecA; // circ + // propeller rotates counter-clockwise viewed from the tail of the aircraft looking forward + cellC2AVecC = cellC2AVecR ^ diskDirNorm; // circ } else if (rotDir == "right") { - // this assumes left hand rotation of propellers - cellC2AVecC = cellC2AVecA ^ cellC2AVecR; // circ + // propeller rotates clockwise viewed from the tail of the aircraft looking forward + cellC2AVecC = diskDirNorm ^ cellC2AVecR; // circ } else { @@ -244,7 +247,7 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource) scalar fRMin = pow(rStarMin, expM) * pow(1.0 - rStarMin, expN); scalar fRMax = pow(rStarMax, expM) * pow(1.0 - rStarMax, expN); - label adjustThrust = diskSubDict.lookupOrDefault