DartAdjointInterface.py

import numpy as np

class DartAdjointInterface():
    def __init__(self, _coupledAdjointSolver, _iSolverAPI, _vSolver):
        self.coupledAdjointSolver = _coupledAdjointSolver
        self.iSolverAPI = _iSolverAPI
        self.vSolver = _vSolver

    def computeDerivatives(self):
        nDim = self.iSolverAPI.solver.pbl.nDim
        nNd = self.iSolverAPI.solver.pbl.msh.nodes.size()
        nEs = len(self.iSolverAPI.iBnd[0].elemsCoord[:,0])
        nNs = len(self.iSolverAPI.iBnd[0].nodesCoord[:,0])
        print(f"{nDim=}")
        print(f"{nNd=}")
        print(f"{nEs=}")
        print(f"{nNs=}")
        dRblw_M = np.zeros((nEs, nNs))
        dRblw_rho = np.zeros((nEs, nNs))
        dRblw_v = np.zeros((nEs, nDim*nNs))

        delta = 1e-4

        saveBlw = np.zeros(nEs)
        for i in range(nEs):
            saveBlw[i] = self.iSolverAPI.blw[0].getU(i)
        
        saveMach = np.zeros(nNs)
        saverho = np.zeros(nNs)
        savev = np.zeros(nDim*nNs)
        for i in range(nNs):
            saveMach[i] = self.iSolverAPI.solver.M[i]
            saverho[i] = self.iSolverAPI.solver.rho[i]
            for j in range(nDim):
                savev[i*nDim + j] = self.iSolverAPI.solver.U[i][j]
        
        # Mach
        for iNo in range(nNs):
            self.iSolverAPI.solver.M[i] += delta
            blowing = self.__runViscous()
            dRblw_M[:,iNo] = (blowing-saveBlw)/delta
            self.iSolverAPI.solver.M[i] = saveMach[i]
        
        # Rho
        for iNo in range(nNs):
            self.iSolverAPI.solver.rho[i] += delta
            blowing = self.__runViscous()
            dRblw_rho[:,iNo] = (blowing-saveBlw)/delta
            self.iSolverAPI.solver.rho[i] = saverho[i]
        
        # V
        for iNo in range(nNs):
            for iDim in range(nDim):
                self.iSolverAPI.solver.U[iNo][iDim] += delta
                blowing = self.__runViscous()
                dRblw_v[:,iNo*nDim + iDim] = (blowing-saveBlw)/delta
                self.iSolverAPI.solver.U[iNo][iDim] = savev[iNo*nDim + iDim]

        # dRblw_M[abs(dRblw_M) < 1e-10] = 0.
        # dRblw_rho[abs(dRblw_rho) < 1e-10] = 0.
        # dRblw_v[abs(dRblw_v) < 1e-10] = 0.
        
        self.coupledAdjointSolver.setdRblw_M(dRblw_M)
        self.coupledAdjointSolver.setdRblw_rho(dRblw_rho)
        self.coupledAdjointSolver.setdRblw_v(dRblw_v)

    def __runViscous(self):
        self.iSolverAPI.imposeInvBC(1)
        # Run viscous solver.
        exitCode = self.vSolver.run(1)
        # Impose blowing boundary condition in the inviscid solver.
        self.iSolverAPI.imposeBlwVel()

        blowing = np.zeros(len(self.iSolverAPI.iBnd[0].elemsCoord[:,0]))
        for i in range(len(self.iSolverAPI.iBnd[0].elemsCoord[:,0])):
            blowing[i] = self.iSolverAPI.blw[0].getU(i)
        return blowing




































    #     nRun = 0

    #     # Store blowing velocity
    #     saveBlw = self.__getBlowingVelocity()

    #     saveMach, saveRho, saveV = self.__getInviscidState()

    #     dRblw_dM = np.zeros((len(saveBlw), 0))
    #     deltaMach = 1e-4
        
    #     # Finite difference
    #     for i in range(len(saveMach)):
    #         self.iSolverAPI.iBnd[0].M[i] = saveMach[i] + deltaMach
    #         self.iSolverAPI.setViscousSolver(1)
    #         self.vSolver.run(1)
    #         nRun+=1
    #         self.iSolverAPI.imposeBlwVel()
    #         newBlowing = self.__getBlowingVelocity()
    #         dRblw_dM = np.column_stack((dRblw_dM, (newBlowing - saveBlw) / deltaMach))
    #         self.iSolverAPI.iBnd[0].M[i] = saveMach[i]

    #     dRblw_dRho = np.zeros((len(saveBlw), 0))
    #     deltaRho = 1e-4
        
    #     # Finite difference
    #     for i in range(len(saveRho)):
    #         self.iSolverAPI.iBnd[0].Rho[i] = saveRho[i] + deltaRho
    #         self.iSolverAPI.setViscousSolver(1)
    #         self.vSolver.run(1)
    #         nRun+=1
    #         self.iSolverAPI.imposeBlwVel()
    #         newBlowing = self.__getBlowingVelocity()
    #         dRblw_dRho = np.column_stack((dRblw_dRho, (newBlowing - saveBlw) / deltaRho))
    #         self.iSolverAPI.iBnd[0].Rho[i] = saveRho[i]

    #     dRblw_dv = np.zeros((len(saveBlw), 0))
    #     deltaV = 1e-4
        
    #     # Finite difference
    #     for j in range(self.iSolverAPI.solver.pbl.nDim):
    #         for i in range(len(saveV[:,0])):
    #             self.iSolverAPI.iBnd[0].V[i,j] = saveV[i,j] + deltaV
    #             self.iSolverAPI.setViscousSolver(1)
    #             self.vSolver.run(1)
    #             nRun+=1
    #             self.iSolverAPI.imposeBlwVel()
    #             newBlowing = self.__getBlowingVelocity()
    #             dRblw_dv = np.column_stack((dRblw_dv, (newBlowing - saveBlw) / deltaV))
    #             self.iSolverAPI.iBnd[0].V[i,j] = saveV[i,j]
        
    #     dRblw_dM[abs(dRblw_dM) < 1e-10] = 0.
    #     dRblw_dRho[abs(dRblw_dRho) < 1e-10] = 0.
    #     dRblw_dv[abs(dRblw_dv) < 1e-10] = 0.

    #     print('Total runs for finite difference', nRun)
    #     self.coupledAdjointSolver.setdRblw_M(dRblw_dM)
    #     self.coupledAdjointSolver.setdRblw_v(dRblw_dv)
    #     self.coupledAdjointSolver.setdRblw_rho(dRblw_dRho)

    # def __getBlowingVelocity(self):
    #     blowingVelocity = np.zeros(0)
    #     for k in range(len(self.iSolverAPI.blw[0].tag.elems)):
    #         blowingVelocity = np.append(blowingVelocity, self.iSolverAPI.blw[0].getU(k))
    #     return blowingVelocity
    
    # def __getInviscidState(self):
    #     V = np.zeros((len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]), 3))
    #     M = np.zeros(len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]))
    #     Rho = np.zeros(len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]))
    #     for n in range(2):
    #         for iRow, row in enumerate(self.iSolverAPI.iBnd[n].nodesCoord[:,3]):
    #             row=int(row)
    #             for iDim in range(3):
    #                 V[iRow,iDim] = self.iSolverAPI.solver.U[row][iDim]
    #             M[iRow] = self.iSolverAPI.solver.M[row]
    #             Rho[iRow] = self.iSolverAPI.solver.rho[row]
    #     return M, Rho, V