Added options to use pre-computed PC mats and udpate part of PC mat v…

…alues when solving the adj.

dafoam/pyDAFoam.py

@@ -465,6 +465,8 @@ def __init__(self):
             "periodicity": -1.0,
             "objFuncStartTime": -1.0,
             "objFuncEndTime": -1.0,
+            "PCMatPrecomputeInterval": 100,
+            "PCMatUpdateInterval": 1,
         }
 
         ## At which iteration should we start the averaging of objective functions.
@@ -861,8 +863,8 @@ def __init__(self, comm=None, options=None):
         # initialize the adjoint vector dict
         self.adjVectors = self._initializeAdjVectors()
 
-        # initialize the dRdWOldTPsi vectors
-        self._initializeTimeAccurateAdjointVectors()
+        # initialize the dRdWTPC dict for unsteady adjoint
+        self.dRdWTPCUnsteady = None
 
         if self.getOption("tensorflow")["active"]:
             TensorFlowHelper.options = self.getOption("tensorflow")
@@ -1016,62 +1018,6 @@ 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"] == "timeAccurate":
-            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"] == "timeAccurate":
-            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
@@ -2337,6 +2283,9 @@ def solveAdjointUnsteady(self):
         self.setOption("runStatus", "solveAdjoint")
         self.updateDAOption()
 
+        PCMatPrecompute = self.getOption("unsteadyAdjoint")["PCMatPrecomputeInterval"]
+        PCMatUpdate = self.getOption("unsteadyAdjoint")["PCMatUpdateInterval"]
+
         self.adjointFail = 0
 
         # update the state to self.solver and self.solverAD
@@ -2370,13 +2319,48 @@ def solveAdjointUnsteady(self):
         # init dRdWTMF
         self.solverAD.initializedRdWTMatrixFree(self.xvVec, self.wVec)
 
-        # calc the preconditioner mat
-        self.dRdWTPC = PETSc.Mat().create(PETSc.COMM_WORLD)
-        self.solver.calcdRdWT(self.xvVec, self.wVec, 1, self.dRdWTPC)
-
-        # Initialize the KSP object
+        # precompute the KSP preconditioner Mat and save them to the self.dRdWTPCUnsteady dict
+        if self.dRdWTPCUnsteady is None:
+
+            self.dRdWTPCUnsteady = {}
+
+            # always calculate the PC mat for the endTime
+            self.solver.setTime(endTime, endTimeIndex)
+            self.solverAD.setTime(endTime, endTimeIndex)
+            # now we can read the variables
+            self.readStateVars(endTime, deltaT)
+            # calc the preconditioner mat for endTime
+            Info("Pre-Computing preconditiner mat for t = %f" % endTime)
+            dRdWTPC = PETSc.Mat().create(PETSc.COMM_WORLD)
+            self.solver.calcdRdWT(self.xvVec, self.wVec, 1, dRdWTPC)
+            # always update the PC mat values using OpenFOAM's fvMatrix
+            self.solver.calcPCMatWithFvMatrix(dRdWTPC)
+            self.dRdWTPCUnsteady[str(endTime)] = dRdWTPC
+
+            # if we define some extra PCMat in PCMatPrecomputeInterval, calculate them here
+            # and set them to the self.dRdWTPCUnsteady dict
+            if objFuncEndTimeIndex > PCMatPrecompute:
+                for timeIndex in range(objFuncEndTimeIndex - 1, 0, -1):
+                    if timeIndex % PCMatPrecompute == 0:
+                        t = timeIndex * deltaT
+                        Info("Pre-Computing preconditiner mat for t = %f" % t)
+                        # read the latest solution
+                        self.solver.setTime(t, timeIndex)
+                        self.solverAD.setTime(t, timeIndex)
+                        # now we can read the variables
+                        self.readStateVars(t, deltaT)
+
+                        # calc the preconditioner mat
+                        dRdWTPC = PETSc.Mat().create(PETSc.COMM_WORLD)
+                        self.solver.calcdRdWT(self.xvVec, self.wVec, 1, dRdWTPC)
+                        # always update the PC mat values using OpenFOAM's fvMatrix
+                        self.solver.calcPCMatWithFvMatrix(dRdWTPC)
+                        self.dRdWTPCUnsteady[str(t)] = dRdWTPC
+
+        # Initialize the KSP object using the PCMat from the endTime
+        PCMat = self.dRdWTPCUnsteady[str(endTime)]
         ksp = PETSc.KSP().create(PETSc.COMM_WORLD)
-        self.solverAD.createMLRKSPMatrixFree(self.dRdWTPC, ksp)
+        self.solverAD.createMLRKSPMatrixFree(PCMat, ksp)
 
         objFuncDict = self.getOption("objFunc")
         designVarDict = self.getOption("designVar")
@@ -2443,6 +2427,16 @@ def solveAdjointUnsteady(self):
                     else:
                         raise Error("ddtSchemeOrder not valid!" % ddtSchemeOrder)
 
+                    # check if we need to update the PC Mat vals or use the pre-computed PC matrix
+                    if str(timeVal) in list(self.dRdWTPCUnsteady.keys()):
+                        Info("Using pre-computed KSP PC mat for %f" % timeVal)
+                        PCMat = self.dRdWTPCUnsteady[str(timeVal)]
+                        self.solverAD.updateKSPPCMat(PCMat, ksp)
+                    if n % PCMatUpdate == 0 and n < endTimeIndex:
+                        # udpate part of the PC mat
+                        Info("Updating dRdWTPC mat value using OF fvMatrix")
+                        self.solver.calcPCMatWithFvMatrix(PCMat)
+
                     # now solve the adjoint eqn
                     self.adjointFail = self.solverAD.solveLinearEqn(ksp, dFdW, self.adjVectors[objFuncName])
 

src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C

@@ -695,6 +695,47 @@ void DASpalartAllmarasFv3::calcLduResidualTurb(volScalarField& nuTildaRes)
     // Below is not necessary, but it doesn't hurt
     nuTildaRes.correctBoundaryConditions();
 }
+
+void DASpalartAllmarasFv3::getFvMatrixFields(
+    const word varName,
+    scalarField& diag,
+    scalarField& upper,
+    scalarField& lower)
+{
+    /* 
+    Description:
+        return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
+        this will be use to compute the preconditioner matrix
+    */
+
+    if (varName != "nuTilda")
+    {
+        FatalErrorIn(
+            "varName not valid. It has to be nuTilda")
+            << exit(FatalError);
+    }
+
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+
+    const volScalarField Stilda(
+        this->fv3(chi, fv1) * ::sqrt(2.0) * mag(skew(fvc::grad(U_)))
+        + this->fv2(chi, fv1) * nuTilda_ / sqr(kappa_ * y_));
+
+    fvScalarMatrix nuTildaEqn(
+        fvm::ddt(phase_, rho_, nuTilda_)
+            + fvm::div(phaseRhoPhi_, nuTilda_, "div(pc)")
+            - fvm::laplacian(phase_ * rho_ * DnuTildaEff(), nuTilda_)
+            - Cb2_ / sigmaNut_ * phase_ * rho_ * magSqr(fvc::grad(nuTilda_))
+        == Cb1_ * phase_ * rho_ * Stilda * nuTilda_ * betaFI_
+            - fvm::Sp(Cw1_ * phase_ * rho_ * fw(Stilda) * nuTilda_ / sqr(y_), nuTilda_));
+
+    nuTildaEqn.relax();
+
+    diag = nuTildaEqn.D();
+    upper = nuTildaEqn.upper();
+    lower = nuTildaEqn.lower();
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H

@@ -161,6 +161,13 @@ public:
 
     /// calculate the turbulence residual using LDU matrix
     virtual void calcLduResidualTurb(volScalarField& nuTildaRes);
+
+    /// return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
+    virtual void getFvMatrixFields(
+        const word varName,
+        scalarField& diag,
+        scalarField& upper,
+        scalarField& lower);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C

@@ -542,6 +542,18 @@ void DATurbulenceModel::rhsSolvePseudoNuTildaEqn(const volScalarField& nuTildaSo
         << abort(FatalError);
 }
 
+/// return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
+void DATurbulenceModel::getFvMatrixFields(
+    const word varName,
+    scalarField& diag,
+    scalarField& upper,
+    scalarField& lower)
+{
+    FatalErrorIn("DATurbulenceModel::getFvMatrixFields")
+        << "Child class not implemented!"
+        << abort(FatalError);
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H

@@ -296,6 +296,13 @@ public:
 
     /// calculate the turbulence residual using LDU matrix
     virtual void calcLduResidualTurb(volScalarField& nuTildaRes);
+
+    /// return the diag(), upper(), and lower() scalarFields from the turbulence model's fvMatrix
+    virtual void getFvMatrixFields(
+        const word varName,
+        scalarField& diag,
+        scalarField& upper,
+        scalarField& lower);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAResidual/DAResidual.C

@@ -155,6 +155,13 @@ void DAResidual::masterFunction(
     }
 }
 
+void DAResidual::calcPCMatWithFvMatrix(Mat PCMat)
+{
+    FatalErrorIn("DAResidual::calcPCMatWithFvMatrix")
+        << "Child class not implemented!"
+        << abort(FatalError);
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/adjoint/DAResidual/DAResidual.H

@@ -120,6 +120,9 @@ public:
         const Vec xvVec,
         const Vec wVec,
         Vec resVec);
+
+    /// calculating the adjoint preconditioner matrix using fvMatrix
+    virtual void calcPCMatWithFvMatrix(Mat PCMat);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //