diff --git a/blast/interfaces/DAdjointInterface.py b/blast/interfaces/DAdjointInterface.py
deleted file mode 100644
index 7eee2fe9576c1630f05d6fa85b04df1833a7edd6..0000000000000000000000000000000000000000
--- a/blast/interfaces/DAdjointInterface.py
+++ /dev/null
@@ -1,566 +0,0 @@
-import numpy as np
-from matplotlib import pyplot as plt
-
-class AdjointInterface():
- def __init__(self, _coupledAdjointSolver, _iSolverAPI, _vSolver):
- self.coupledAdjointSolver = _coupledAdjointSolver
- self.isol = _iSolverAPI
- self.vsol = _vSolver
-
- def build(self):
- nDim = self.isol.solver.pbl.nDim
- nNd = self.isol.solver.pbl.msh.nodes.size()
-
- nNs_v = 0
- for isec in self.isol.sec[0]:
- nNs_v += isec.nNodes
- nEs_v = 0
- for isec in self.isol.sec[0]:
- nEs_v += isec.nElms
-
- dRdStar_dStar = np.zeros((nNs_v, nNs_v))
- dRdStar_M = np.zeros((nNs_v, nNs_v))
- dRdStar_rho = np.zeros((nNs_v, nNs_v))
- dRdStar_v = np.zeros((nNs_v, nNs_v*nDim))
- dRdStar_x = np.zeros((nNs_v, nNd*nDim))
-
- dRblw_rho = np.zeros((nEs_v, nNs_v))
- dRblw_v = np.zeros((nEs_v, nNs_v*nDim))
- dRblw_dStar = np.zeros((nEs_v, nNs_v))
- dRblw_x = np.zeros((nEs_v, nNd*nDim))
-
- delta = 1e-6
- delta_x = 1e-8
-
- # Compute dVt/dV
- dVt_v = np.zeros((nNs_v, nNs_v*nDim))
- v = np.zeros((nNs_v, nDim))
- for d in range(nDim):
- v[:self.isol.sec[0][0].nNodes, d] = self.isol.vBnd[0][0].getVelocity("upper", d)
- v[self.isol.sec[0][0].nNodes:self.isol.sec[0][0].nNodes+self.isol.sec[0][1].nNodes, d] = self.isol.vBnd[0][0].getVelocity("lower", d)
- v[self.isol.sec[0][0].nNodes+self.isol.sec[0][1].nNodes:, d] = self.isol.vBnd[0][1].getVelocity("wake", d)
-
- i = 0
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- for d in range(nDim):
- dVt_v[i, i*nDim+d] = 1/np.sqrt(v[i, 0]**2 + v[i, 1]**2)*v[i, d]
- i += 1
-
- """# Finite difference
- i = 0
- offset = 0
- dVt_v_fd = np.zeros((nNs_v, nNs_v*nDim))
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- for d in range(nDim):
- save = v[i, d]
- v[i, d] = save + delta
- vt1 = np.sqrt(v[:,0]**2 + v[:,1]**2)
- v[i, d] = save - delta
- vt2 = np.sqrt(v[:,0]**2 + v[:,1]**2)
- dVt_v_fd[:, i*nDim+d] = (vt1 - vt2)/(2*delta)
- v[i, d] = save
- i += 1"""
-
- # dRdStar_M
- i = 0
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- mach_lo = isec.Me[j]
- isec.Me[j] = mach_lo + delta
- deltaStar1 = self.__runViscous()
- isec.Me[j] = mach_lo - delta
- deltaStar2 = self.__runViscous()
- isec.Me[j] = mach_lo
- dRdStar_M[:, i] = -(deltaStar1 - deltaStar2)/(2.*delta)
- i += 1
-
- # dRdStar_dStar
- dRdStar_dStar_test = np.zeros((nNs_v, nNs_v))
- i = 0
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- dStar_lo = isec.deltaStarExt[j]
- isec.deltaStarExt[j] = dStar_lo + delta
- deltaStar1 = self.__runViscous()
- isec.deltaStarExt[j] = dStar_lo - delta
- deltaStar2 = self.__runViscous()
- isec.deltaStarExt[j] = dStar_lo
- dRdStar_dStar_test[:, i] = (deltaStar1 - deltaStar2)/(2.*delta)
- i += 1
-
- # Set dRdStar_dStar to identity
- dRdStar_dStar = np.eye(nNs_v) - dRdStar_dStar_test
-
- # dRdStar_vt
- dRdStar_vt = np.zeros((nNs_v, nNs_v))
- i = 0
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- save = isec.vt[j]
- isec.vt[j] = save + delta
- #isec.vtExt[j] = save + delta
- deltaStar1 = self.__runViscous()
- isec.vt[j] = save - delta
- #isec.vtExt[j] = save - delta
- deltaStar2 = self.__runViscous()
- isec.vt[j] = save
- #isec.vtExt[j] = save
- dRdStar_vt[:, i] = -(deltaStar1 - deltaStar2)/(2.*delta)
- i += 1
-
- dRdStar_v = dRdStar_vt @ dVt_v
-
- """# dRdStar_v fd
- i = 0
- offset = 0
- dRdStar_v_fd = np.zeros((nNs_v, nNs_v*nDim))
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- for d in range(nDim):
- save = v[i, d]
- v[i, d] = save + delta
- # Recompute vt
- vt = np.sqrt(v[:, 0]**2 + v[:, 1]**2)
- # Change vt in c++ object
- # Resplit vt in the sections
- for nodesSec in range(isec.nNodes):
- #print('nodesSec', nodesSec, 'isec.vt', isec.vt[nodesSec], 'vt', vt[offset + nodesSec])
- isec.vt[nodesSec] = vt[offset + nodesSec]
- deltaStar1 = self.__runViscous()
- v[i, d] = save - delta
- # Recompute vt
- vt = np.sqrt(v[:, 0]**2 + v[:, 1]**2)
- # Change vt in c++ object
- for nodesSec in range(isec.nNodes):
- isec.vt[nodesSec] = vt[offset + nodesSec]
- deltaStar2 = self.__runViscous()
- v[i, d] = save
- dRdStar_v_fd[:, i*nDim+d] = -(deltaStar1 - deltaStar2)/(2*delta)
- i += 1
- offset += isec.nNodes
- print('diff', dRdStar_v - dRdStar_v_fd)
- print('max diff dRdstar_v', np.max(np.abs(dRdStar_v - dRdStar_v_fd)))
- quit()"""
-
- # dRdStar_x
- print('dRdStar_x')
- for isec in self.isol.sec[0]:
- if isec.name == 'upper' or isec.name == 'lower':
- iReg = 0
- elif isec.name == 'wake':
- iReg = 1
- for j in range(isec.nNodes):
- rowj = self.isol.blw[iReg].nodes[isec.rows[j]].row
- save = isec.x[j]
-
- isec.x[j] = save + delta_x
- isec.computeLocalCoordinates()
- isec.xxExt = isec.loc
-
- deltaStar1 = self.__runViscous()
-
- isec.x[j] = save - delta_x
- isec.computeLocalCoordinates()
- isec.xxExt = isec.loc
-
- deltaStar2 = self.__runViscous()
-
- diff = (deltaStar1 - deltaStar2)/(2.*delta_x)
-
- isec.x[j] = save
- isec.computeLocalCoordinates()
- isec.xxExt = isec.loc
-
- dRdStar_x[:, rowj*nDim+0] -= diff
-
- save = isec.y[j]
- isec.y[j] = save + delta_x
- isec.computeLocalCoordinates()
- isec.xxExt = isec.loc
-
- deltaStar1 = self.__runViscous()
-
- isec.y[j] = save - delta_x
- isec.computeLocalCoordinates()
- isec.xxExt = isec.loc
-
- deltaStar2 = self.__runViscous()
-
- diff = (deltaStar1 - deltaStar2)/(2.*delta_x)
-
- isec.y[j] = save
- isec.computeLocalCoordinates()
- isec.xxExt = isec.loc
- dRdStar_x[:, rowj*nDim+1] -= diff
-
- print('done')
- """# Finite difference dRdStar_x_fd_2
- dRdStar_x_fd = np.zeros((nNs_v, nNd*nDim))
- saveNodes = np.zeros((len(self.isol.iobj['pbl'].msh.nodes), self.isol.getnDim()))
- for n in self.isol.iobj['pbl'].msh.nodes:
- for i in range(self.isol.getnDim()):
- saveNodes[n.row, i] = n.pos[i]
- for ib, b in enumerate(self.isol.blw):
- for inod, n in enumerate(b.nodes):
- # Find all the nodes in the viscous mesh corresponding to the current node. We also save the boundary index (0: airfoil, 1: wake)
- vidx = []
- try:
- vidx.append((np.where(self.isol.vBnd[0][0].nodesCoord[:,3] == n.row)[0][0], 0))
- except IndexError:
- pass
- try:
- vidx.append((np.where(self.isol.vBnd[0][1].nodesCoord[:,3] == n.row)[0][0], 1))
- except IndexError:
- pass
- for idim in range(self.isol.getnDim()):
- ddStar = [np.zeros(nNs_v), np.zeros(nNs_v)]
- for isgn, sign in enumerate([-1, 1]):
- # Modify all the nodes of the viscous mesh (contained in vidx)
- for (ii, boundary) in vidx:
- self.isol.vBnd[0][boundary].nodesCoord[ii, idim] = saveNodes[n.row, idim] + sign * delta
- self.setvMesh(self.isol.vBnd, self.isol.sec[0])
-
- ddStar[isgn] = self.__runViscous()
-
- # Reset mesh
- for vb in self.isol.vBnd[0]:
- for jj, jrow in enumerate(vb.nodesCoord[:,3]):
- for d in range(self.isol.getnDim()):
- vb.nodesCoord[jj, d] = saveNodes[int(jrow), d]
- dRdStar_x_fd[:, n.row*self.isol.getnDim()+idim] = -(ddStar[1] - ddStar[0])/(2*delta)
- np.savetxt('dRdStar_x_fd.txt', dRdStar_x_fd)
- np.savetxt('dRdStar_x.txt', dRdStar_x)
- print('max diff', np.max(np.abs(dRdStar_x - dRdStar_x_fd)))
- quit()"""
-
- # dRblw_rho
- offsetNodes = 0
- offsetElems = 0
- for isec in self.isol.sec[0]:
- test = np.array(isec.evalGradwrtRho())
- dRblw_rho[offsetElems:offsetElems+isec.nElms, offsetNodes:offsetNodes+isec.nNodes] = -test
- offsetNodes += isec.nNodes
- offsetElems += isec.nElms
-
- """# Finite difference
- start = time.time()
- i = 0
- offset = 0
- dRblw_rho_fd = np.zeros((nEs_v, nNs_v))
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- save = isec.rhoe[j]
- isec.rhoe[j] = save + delta
- isec.computeBlowingVelocity()
- blowing_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing1 = np.concatenate([blowing_part1, blowing_part2, blowing_part3])
- isec.rhoe[j] = save - delta
- isec.computeBlowingVelocity()
- blowing_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing2 = np.concatenate([blowing_part1, blowing_part2, blowing_part3])
- isec.rhoe[j] = save
- dRblw_rho_fd[:, i] = (blowing1 - blowing2)/(2*delta)
- i += 1
- offset += isec.nNodes
- end = time.time()
- print('done dRblw_rho_fd in', end-start)"""
-
- # dRblw_vt
- offsetNodes = 0
- offsetElems = 0
- dRblw_vt = np.zeros((nEs_v, nNs_v))
- for isec in self.isol.sec[0]:
- test = np.array(isec.evalGradwrtVt())
- dRblw_vt[offsetElems:offsetElems+isec.nElms, offsetNodes:offsetNodes+isec.nNodes] = -test
- offsetNodes += isec.nNodes
- offsetElems += isec.nElms
-
- #print('dRblw_vt', dRblw_vt)
-
- # dRblw_v = dRblw_vt*dVt_v
- dRblw_v = dRblw_vt @ dVt_v
-
- # Validate by finite difference
- """i = 0
- offset = 0
- dRblw_v_fd = np.zeros((nEs_v, nNs_v*nDim))
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- for d in range(nDim):
- save = v[i, d]
- v[i, d] = save + delta
- # Recompute vt
- vt = np.sqrt(v[:, 0]**2 + v[:, 1]**2)
- # Change vt in c++ object
- # Resplit vt in the sections
- for nodesSec in range(isec.nNodes):
- #print('nodesSec', nodesSec, 'isec.vt', isec.vt[nodesSec], 'vt', vt[offset + nodesSec])
- isec.vt[nodesSec] = vt[offset + nodesSec]
-
- isec.computeBlowingVelocity()
- blowing_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing1 = np.concatenate([blowing_part1, blowing_part2, blowing_part3])
- v[i, d] = save - delta
- # Recompute vt
- vt = np.sqrt(v[:, 0]**2 + v[:, 1]**2)
- # Change vt in c++ object
- for nodesSec in range(isec.nNodes):
- isec.vt[nodesSec] = vt[offset + nodesSec]
- isec.computeBlowingVelocity()
- blowing_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing2 = np.concatenate([blowing_part1, blowing_part2, blowing_part3])
- v[i, d] = save
- dRblw_v_fd[:, i*nDim+d] = (blowing1 - blowing2)/(2*delta)
- i += 1
- offset += isec.nNodes"""
-
- # dRblw_dStar
- offsetNodes = 0
- offsetElems = 0
- for isec in self.isol.sec[0]:
- test = np.array(isec.evalGradwrtDeltaStar())
- dRblw_dStar[offsetElems:offsetElems+isec.nElms, offsetNodes:offsetNodes+isec.nNodes] = -test
- offsetNodes += isec.nNodes
- offsetElems += isec.nElms
-
- """# Finite difference
- dRblw_dStar_fd = np.zeros((nEs_v, nNs_v))
- i = 0
- for isec in self.isol.sec[0]:
- for j in range(isec.nNodes):
- save = isec.rhoe[j]
- isec.deltaStar[j] = save + delta
- isec.computeBlowingVelocity()
- blowing_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing1 = np.concatenate([blowing_part1, blowing_part2, blowing_part3])
- isec.deltaStar[j] = save - delta
- isec.computeBlowingVelocity()
- blowing_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing2 = np.concatenate([blowing_part1, blowing_part2, blowing_part3])
- isec.deltaStar[j] = save
- dRblw_dStar_fd[:, i] = -(blowing1 - blowing2)/(2*delta)
- i += 1
- print('diff', dRblw_dStar - dRblw_dStar_fd)
- print('max diff', np.max(np.abs(dRblw_dStar - dRblw_dStar_fd)))
- quit()"""
-
-
- """print('Max error dRblw_v', np.max(np.abs(dRblw_v - dRblw_v_fd)))
- print('Max error dRblw_rho', np.max(np.abs(dRblw_rho - dRblw_rho_fd)))
- print('Max error dRblw_dStar', np.max(np.abs(dRblw_dStar - dRblw_dStar_fd)))"""
-
- # dRblw_x
- offsetElems = 0
- dRblw_x = np.zeros((nEs_v, nNd*nDim))
- for isec in self.isol.sec[0]:
- if isec.name == 'upper' or isec.name == 'lower':
- iReg = 0
- elif isec.name == 'wake':
- iReg = 1
- test_x = -np.array(isec.evalGradwrtX())
- test_y = -np.array(isec.evalGradwrtY())
- for j in range(isec.nNodes):
- rowj = self.isol.blw[iReg].nodes[isec.rows[j]].row
- dRblw_x[offsetElems:offsetElems+isec.nElms, rowj*nDim+0] = test_x[:, j]
- dRblw_x[offsetElems:offsetElems+isec.nElms, rowj*nDim+1] = test_y[:, j]
- offsetElems += isec.nElms
-
-
- """# Finite difference dRblw_x_fd_2
- dRblw_x_fd = np.zeros((nEs_v, nNd*nDim))
- saveNodes = np.zeros((len(self.isol.iobj['pbl'].msh.nodes), self.isol.getnDim()))
- for n in self.isol.iobj['pbl'].msh.nodes:
- for i in range(self.isol.getnDim()):
- saveNodes[n.row, i] = n.pos[i]
- for ib, b in enumerate(self.isol.blw):
- for inod, n in enumerate(b.nodes):
- # Find all the nodes in the viscous mesh corresponding to the current node. We also save the boundary index (0: airfoil, 1: wake)
- vidx = []
- try:
- vidx.append((np.where(self.isol.vBnd[0][0].nodesCoord[:,3] == n.row)[0][0], 0))
- except IndexError:
- pass
- try:
- vidx.append((np.where(self.isol.vBnd[0][1].nodesCoord[:,3] == n.row)[0][0], 1))
- except IndexError:
- pass
- for idim in range(self.isol.getnDim()):
- dblw = [np.zeros(nEs_v), np.zeros(nEs_v)]
- for isgn, sign in enumerate([-1, 1]):
- # Modify all the nodes of the viscous mesh (contained in vidx)
- for (ii, boundary) in vidx:
- self.isol.vBnd[0][boundary].nodesCoord[ii, idim] = saveNodes[n.row, idim] + sign * delta
- self.setvMesh(self.isol.vBnd, self.isol.sec[0])
- for reg in self.isol.sec[0]:
- reg.computeBlowingVelocity()
-
- blowing1_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing1_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing1_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- dblw[isgn] = np.concatenate([blowing1_part1, blowing1_part2, blowing1_part3])
-
- # Reset mesh
- for vb in self.isol.vBnd[0]:
- for jj, jrow in enumerate(vb.nodesCoord[:,3]):
- for d in range(self.isol.getnDim()):
- vb.nodesCoord[jj, d] = saveNodes[int(jrow), d]
- dRblw_x_fd[:, n.row*self.isol.getnDim()+idim] = -(dblw[1] - dblw[0])/(2*delta)
- # Write to file
- np.savetxt('dRblw_x_fd.txt', dRblw_x_fd)
- np.savetxt('dRblw_x.txt', dRblw_x)
- print('max diff', np.max(np.abs(dRblw_x - dRblw_x_fd)))
- quit()"""
-
- """# Finite difference dRblw_x_fd
- offsetElems = 0
- dRblw_x_fd = np.zeros((nEs_v, nNd*nDim))
- for isec in self.isol.sec[0]:
- if isec.name == 'upper' or isec.name == 'lower':
- iReg = 0
- elif isec.name == 'wake':
- iReg = 1
- for j in range(isec.nNodes):
- rowj = self.isol.blw[iReg].nodes[isec.rows[j]].row
- save = isec.x[j]
- isec.x[j] = save + delta
- isec.computeLocalCoordinates()
- isec.computeBlowingVelocity()
-
- blowing1_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing1_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing1_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing1 = np.concatenate([blowing1_part1, blowing1_part2, blowing1_part3])
-
- isec.x[j] = save - delta
- isec.computeLocalCoordinates()
- isec.computeBlowingVelocity()
-
- blowing2_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing2_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing2_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing2 = np.concatenate([blowing2_part1, blowing2_part2, blowing2_part3])
- bb = (blowing1 - blowing2)/(2*delta)
-
- isec.x[j] = save
- dRblw_x_fd[offsetElems:offsetElems+isec.nElms, rowj*nDim+0] = -bb[offsetElems:offsetElems+isec.nElms]
-
- if rowj == 0:
- print(isec.name, j, -(blowing1 - blowing2)/(2*delta))
-
- save = isec.y[j]
- isec.y[j] = save + delta
- isec.computeLocalCoordinates()
- isec.computeBlowingVelocity()
-
- blowing1_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing1_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing1_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing1 = np.concatenate([blowing1_part1, blowing1_part2, blowing1_part3])
-
- isec.y[j] = save - delta
- isec.computeLocalCoordinates()
- isec.computeBlowingVelocity()
-
- blowing2_part1 = np.asarray(self.isol.sec[0][0].blowingVelocity)
- blowing2_part2 = np.asarray(self.isol.sec[0][1].blowingVelocity)
- blowing2_part3 = np.asarray(self.isol.sec[0][2].blowingVelocity)
- blowing2 = np.concatenate([blowing2_part1, blowing2_part2, blowing2_part3])
- bb = (blowing1 - blowing2)/(2*delta)
-
- isec.y[j] = save
- dRblw_x_fd[offsetElems:offsetElems+isec.nElms, rowj*nDim+1] = -bb[offsetElems:offsetElems+isec.nElms]
- offsetElems += isec.nElms
- for row in range(len(dRblw_x[:,0])):
- print(dRblw_x[row,:] - dRblw_x_fd[row,:])
- print(
- 'max diff', np.max(np.abs(dRblw_x - dRblw_x_fd)))
- quit()"""
-
- print('Imposing matrix')
- self.coupledAdjointSolver.setdRdStar_M(dRdStar_M)
- print('dRdStar_M done')
- self.coupledAdjointSolver.setdRdStar_v(dRdStar_v)
- print('dRdStar_v done')
- self.coupledAdjointSolver.setdRdStar_dStar(dRdStar_dStar)
- print('dRdStar_dStar done')
- self.coupledAdjointSolver.setdRdStar_x(dRdStar_x)
- print('dRdStar_x done')
-
-
- self.coupledAdjointSolver.setdRblw_rho(dRblw_rho)
- print('dRblw_rho done')
- self.coupledAdjointSolver.setdRblw_v(dRblw_v)
- print('dRblw_v done')
- self.coupledAdjointSolver.setdRblw_dStar(dRblw_dStar)
- print('dRblw_dStar done')
- self.coupledAdjointSolver.setdRblw_x(dRblw_x)
- print('dRblw_x done')
-
- def __getDeltaStar(self):
- deltaStar_part1 = np.asarray(self.isol.sec[0][0].deltaStar)
- deltaStar_part2 = np.asarray(self.isol.sec[0][1].deltaStar)
- deltaStar_part3 = np.asarray(self.isol.sec[0][2].deltaStar)
- deltaStar = np.concatenate([deltaStar_part1, deltaStar_part2, deltaStar_part3])
- #deltaStar = np.concatenate([deltaStar_part1, deltaStar_part2])
- return deltaStar
-
- def setvMesh(self, vBnd, secs):
- import tbox
- # Compute section chord
- xMin = np.min(vBnd[0][0].nodesCoord[:,0])
- xMax = np.max(vBnd[0][0].nodesCoord[:,0])
- chord = xMax - xMin
- for reg in secs:
- if reg.name == 'wake':
- iReg = 1
- elif reg.name == 'lower' or reg.name == 'upper':
- iReg = 0
- else:
- raise RuntimeError('Invalid region', reg.name)
-
- nodeRows = tbox.std_vector_int()
- for val in vBnd[0][iReg].getNodesRow(reg.name):
- nodeRows.push_back(int(val))
- connectElems = tbox.std_vector_int()
- for val in vBnd[0][iReg].getConnectList(reg.name):
- connectElems.push_back(int(val))
- # x, y, z
- xv = tbox.std_vector_double()
- for val in vBnd[0][iReg].getNodesCoord(reg.name, 0):
- xv.push_back(val)
- yv = tbox.std_vector_double()
- for val in vBnd[0][iReg].getNodesCoord(reg.name, 1):
- yv.push_back(val)
- zv = tbox.std_vector_double()
- for val in vBnd[0][iReg].getNodesCoord(reg.name, 2):
- zv.push_back(val)
-
- # Set mesh
- reg.setMesh(xv,
- yv,
- zv,
- chord,
- xMin,
- nodeRows,
- connectElems)
-
- def __runViscous(self):
- ext = self.vsol.run()
- if ext != 0:
- raise ValueError('Viscous solver did not converge')
- deltaStar = self.__getDeltaStar()
- return deltaStar
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index e71203795e1e1969ead3c77129565d7cbe52116f..8ea9318add45b1a8ec63c746970f578903e819fb 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -72,11 +72,6 @@ void CoupledAdjoint::reset()
for (auto bBC : isol->pbl->bBCs)
nEs += bBC->uE.size(); // Number of blowing elements
- std::cout << "nDim " << nDim << std::endl;
- std::cout << "nNd " << nNd << std::endl;
- std::cout << "nNs " << nNs << std::endl;
- std::cout << "nEs " << nEs << std::endl;
-
dRphi_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNd);
dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
@@ -130,8 +125,12 @@ void CoupledAdjoint::reset()
dRphi_AoA = Eigen::VectorXd::Zero(nNd);
// Objective functions
+ // Cl = f(phi) & Cd = Cdp + Cdf = f(phi, M, v)
dCl_phi = Eigen::VectorXd::Zero(nNd);
dCd_phi = Eigen::VectorXd::Zero(nNd);
+ dCd_M = Eigen::VectorXd::Zero(nNs);
+ dCd_v = Eigen::VectorXd::Zero(nDim*nNs);
+ dCd_dStar = Eigen::VectorXd::Zero(nNs);
// Angle of attack
dCl_AoA = 0.0;
@@ -142,7 +141,8 @@ void CoupledAdjoint::reset()
// Mesh
dRx_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNd, nDim*nNd);
dCl_x = Eigen::VectorXd::Zero(nDim*nNd);
- dCd_x = Eigen::VectorXd::Zero(nDim*nNd);
+ dCdp_x = Eigen::VectorXd::Zero(nDim*nNd);
+ dCdf_x = Eigen::VectorXd::Zero(nDim*nNd);
tdCl_x.resize(nNd, Eigen::Vector3d::Zero());
tdCd_x.resize(nNd, Eigen::Vector3d::Zero());
@@ -192,17 +192,35 @@ void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMat
void CoupledAdjoint::run()
{
+ tms["0-adj_total"].start();
+ tms["1-adj_inviscid"].start();
gradientswrtInviscidFlow();
+ tms["1-adj_inviscid"].stop();
+ tms["2-adj_mach"].start();
gradientswrtMachNumber();
+ tms["2-adj_mach"].stop();
+ tms["3-adj_density"].start();
gradientswrtDensity();
+ tms["3-adj_density"].stop();
+ tms["4-adj_velocity"].start();
gradientswrtVelocity();
+ tms["4-adj_velocity"].stop();
+ tms["5-adj_dStar"].start();
gradientswrtDeltaStar();
+ tms["5-adj_dStar"].stop();
+ tms["6-adj_blw"].start();
gradientswrtBlowingVelocity();
+ tms["6-adj_blw"].stop();
+ tms["7-adj_aoa"].start();
gradientwrtAoA();
+ tms["7-adj_aoa"].stop();
+ tms["8-adj_mesh"].start();
gradientwrtMesh();
+ tms["8-adj_mesh"].stop();
transposeMatrices();
+ tms["9-build"].start();
// Store pointers to the matrices in a 2D vector
std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> matrices = {
{&dRphi_phi, &dRM_phi, &dRrho_phi, &dRv_phi, &dRdStar_phi, &dRblw_phi},
@@ -218,6 +236,7 @@ void CoupledAdjoint::run()
for (const auto& mat1 : matrices[0]) cols += mat1->cols(); // All matrices in a column have the same number of columns
Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols);
buildAdjointMatrix(matrices, A_adjoint);
+ tms["9-build"].stop();
size_t phiIdx = 0;
size_t machIdx = dRphi_phi.cols();
@@ -229,13 +248,15 @@ void CoupledAdjoint::run()
//***** Gradient wrt angle of attack *****//
//****************************************//
// CL
+ tms["10-solveCxAoA"].start();
std::vector<double> rhsCl(A_adjoint.cols(), 0.0);
- for (auto i = 0; i<dCl_phi.size(); ++i)
+ for (auto i = phiIdx; i<phiIdx+dCl_phi.size(); ++i)
rhsCl[i] = dCl_phi[i]; // Only dCl_dphi is non-zero
Eigen::VectorXd rhsCl_ = Eigen::Map<const Eigen::VectorXd>(rhsCl.data(), rhsCl.size());
Eigen::VectorXd lambdaCl(A_adjoint.cols()); // Solution vector for lambdasCl
Eigen::Map<Eigen::VectorXd> lambdaCl_(lambdaCl.data(), lambdaCl.size());
+
alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()), lambdaCl_);
// Total gradient
@@ -243,17 +264,40 @@ void CoupledAdjoint::run()
// CD
std::vector<double> rhsCd(A_adjoint.cols(), 0.0);
- for (auto i = 0; i<dCd_phi.size(); ++i)
- rhsCd[i] = dCd_phi[i];
+ int jj = 0;
+ for (auto i = phiIdx; i<phiIdx+dCd_phi.size(); ++i)
+ {
+ rhsCd[i] = dCd_phi[jj];
+ jj++;
+ }
+ jj = 0;
+ for (auto i = machIdx; i<machIdx+dCd_M.size(); ++i)
+ {
+ rhsCd[i] = dCd_M[jj];
+ jj++;
+ }
+ jj = 0;
+ for (auto i = vIdx; i<vIdx+dCd_v.size(); ++i)
+ {
+ rhsCd[i] = dCd_v[jj];
+ jj++;
+ }
+ jj = 0;
+ for (auto i = dStarIdx; i<dStarIdx+dCd_dStar.size(); ++i)
+ {
+ rhsCd[i] = dCd_dStar[jj];
+ jj++;
+ }
Eigen::VectorXd rhsCd_ = Eigen::Map<const Eigen::VectorXd>(rhsCd.data(), rhsCd.size());
Eigen::VectorXd lambdaCd(A_adjoint.cols()); // Solution vector for lambdasCd
Eigen::Map<Eigen::VectorXd> lambdaCd_(lambdaCd.data(), lambdaCd.size());
- // Solve linear system
+
alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCd_.data(), rhsCd_.size()), lambdaCd_);
// Total gradient
tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Total gradient
+ tms["10-solveCxAoA"].stop();
//***** Gradient wrt mesh coordinates *****//
//*****************************************//
@@ -261,6 +305,7 @@ void CoupledAdjoint::run()
// Solution lambdaCl_x of (from augmented Lagrangian, eqn d/dx )
// "dCk_x + dRphi_x * lambdaCk_phi + dRdStar_x * lambdaCk_x + dRblw_x * lambdaCk_blw + dRx_x * lambdaCk_x = 0"
+ tms["12-solveCxMesh"].start();
Eigen::VectorXd rhsCl_x = dCl_x;
rhsCl_x -= dRphi_x.transpose() * lambdaCl.segment(phiIdx, dRphi_phi.cols()); // Potential contribution
rhsCl_x -= dRM_x.transpose() * lambdaCl.segment(machIdx, dRM_M.cols()); // Mach number contribution
@@ -273,7 +318,8 @@ void CoupledAdjoint::run()
Eigen::Map<Eigen::VectorXd> lambdaCl_x_(lambdaCl_x.data(), lambdaCl_x.size());
alinsol->compute(dRx_x.transpose(), Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size()), lambdaCl_x_);
- Eigen::VectorXd rhsCd_x = dCd_x;
+ Eigen::VectorXd rhsCd_x = dCdp_x; // Pressure drag coefficient contribution
+ rhsCd_x += dCdf_x; // Friction drag coefficient contribution
rhsCd_x -= dRphi_x.transpose() * lambdaCd.segment(phiIdx, dRphi_phi.cols()); // Potential contribution
rhsCd_x -= dRM_x.transpose() * lambdaCd.segment(machIdx, dRM_M.cols()); // Mach number contribution
rhsCd_x -= dRrho_x.transpose() * lambdaCd.segment(rhoIdx, dRrho_rho.cols()); // Density contribution
@@ -293,6 +339,8 @@ void CoupledAdjoint::run()
tdCd_x[n->row](m) = lambdaCd_x[isol->pbl->nDim * (n->row) + m];
}
}
+ tms["12-solveCxMesh"].stop();
+ //std::cout << tms << std::endl;
}
/**
@@ -363,6 +411,57 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
// Remove zeros
dRM_phi.prune(0.); dRrho_phi.prune(0.); dRv_phi.prune(0.);
+ // // Analytical
+ // std::vector<Eigen::Triplet<double>> tripletsdM_an;
+
+ // size_t jj = 0;
+ // int iReg = 0;
+ // for (auto sec : vsol->sections[0])
+ // {
+ // if (sec->name == "wake")
+ // iReg = 1;
+ // else if (sec->name == "upper" || sec->name == "lower")
+ // iReg = 0;
+ // else
+ // throw std::runtime_error("gradientswrtInviscidFlow: Unknown section name");
+ // for (size_t iNode = 0; iNode < sec->nNodes; ++iNode)
+ // {
+ // tbox::Node* srfNode = pbl->bBCs[iReg]->nodes[sec->rows[iNode]];
+ // for (auto e: pbl->fluid->neMap[srfNode])
+ // {
+ // // Get pre-computed values
+ // tbox::Cache &cache = e->getVCache();
+ // tbox::Gauss &gauss = cache.getVGauss();
+ // Eigen::VectorXd grad = Eigen::VectorXd::Zero(e->nodes.size());
+ // for (size_t k = 0; k < gauss.getN(); ++k)
+ // {
+ // // Gauss point and determinant
+ // double wdj = gauss.getW(k) * e->getDetJ(k);
+ // // Shape functions and solution gradient
+ // Eigen::MatrixXd const &dSf = e->getJinv(k) * cache.getDsf(k);
+ // Eigen::VectorXd dPhi = e->computeGradient(isol->phi, k);
+ // grad += isol->pbl->fluid->mach->evalGrad(*e, isol->phi, k) * dSf.transpose() * dPhi;
+ // }
+ // for (size_t k = 0; k < e->nodes.size(); ++k)
+ // {
+ // tripletsdM_an.push_back(Eigen::Triplet<double>(jj, e->nodes[k]->row, -grad(k)/e->nodes.size()));
+ // }
+ // }
+ // ++jj;
+ // }
+ // }
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi_an(dRM_phi.rows(), dRM_phi.cols());
+ // dRM_phi_an.setFromTriplets(tripletsdM_an.begin(), tripletsdM_an.end());
+ // dRM_phi_an.prune(0.);
+ // writeMatrixToFile(dRM_phi_an, "dRM_phi_an.txt");
+ // writeMatrixToFile(dRM_phi, "dRM_phi.txt");
+ // writeMatrixToFile((dRM_phi_an - dRM_phi), "dif.txt");
+ // std::cout << "written" << std::endl;
+ // std::cout << dRM_phi_an - dRM_phi << std::endl;
+
+ // throw;
+
+
//**** partials dCl_phi, dCd_phi ****//
std::vector<std::vector<double>> dCx_phi = iadjoint->computeGradientCoefficientsFlow();
dCl_phi = Eigen::VectorXd::Map(dCx_phi[0].data(), dCx_phi[0].size());
@@ -385,6 +484,7 @@ void CoupledAdjoint::gradientswrtMachNumber()
size_t i = 0;
double saveM = 0.;
std::vector<Eigen::VectorXd> ddStar(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ std::vector<double> dCdf(2, 0.);
for (auto sec: vsol->sections[0])
for (size_t j = 0; j < sec->nNodes; ++j)
{
@@ -392,13 +492,17 @@ void CoupledAdjoint::gradientswrtMachNumber()
sec->Me[j] = saveM + delta;
ddStar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
sec->Me[j] = saveM - delta;
ddStar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
sec->Me[j] = saveM;
for (int k = 0; k < ddStar[0].size(); ++k)
T.push_back(Eigen::Triplet<double>(k, i, -(ddStar[0][k] - ddStar[1][k])/(2.*delta)));
+ // Collect Cd gradients
+ dCd_M(i) = (dCdf[0] - dCdf[1])/(2.*delta);
++i;
}
dRdStar_M.setFromTriplets(T.begin(), T.end());
@@ -501,6 +605,8 @@ void CoupledAdjoint::gradientswrtVelocity()
i = 0;
double saveM = 0.;
std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ std::vector<double> dCdf(2, 0.);
+ Eigen::VectorXd dCd_vt = Eigen::VectorXd::Zero(dRM_M.cols());
for (auto sec: vsol->sections[0])
for (size_t j = 0; j < sec->nNodes; ++j)
{
@@ -508,13 +614,16 @@ void CoupledAdjoint::gradientswrtVelocity()
sec->vt[j] = saveM + delta;
dvt[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
sec->vt[j] = saveM - delta;
dvt[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
sec->vt[j] = saveM;
for (int k = 0; k < dvt[0].size(); ++k)
T.push_back(Eigen::Triplet<double>(k, i, -(dvt[0][k] - dvt[1][k])/(2.*delta)));
+ dCd_vt(i) = (dCdf[0] - dCdf[1])/(2.*delta);
++i;
}
Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(), dRdStar_M.cols());
@@ -525,6 +634,8 @@ void CoupledAdjoint::gradientswrtVelocity()
dRdStar_v = dRdStar_vt * dVt_v;
dRdStar_v.prune(0.);
+ //**** dRCdf_v ****//
+ dCd_v = dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim]
//**** dRblw_vt (analytical) ****//
std::deque<Eigen::Triplet<double>> T_blw;
@@ -561,6 +672,7 @@ void CoupledAdjoint::gradientswrtDeltaStar()
double delta = 1e-6;
double saveDStar = 0.;
std::vector<Eigen::VectorXd> ddstar(2, Eigen::VectorXd::Zero(dRdStar_dStar.cols()));
+ std::vector<double> dCdf(2, 0.);
size_t i = 0;
for (auto sec: vsol->sections[0])
for (size_t j = 0; j < sec->nNodes; ++j)
@@ -569,13 +681,16 @@ void CoupledAdjoint::gradientswrtDeltaStar()
sec->deltaStarExt[j] = saveDStar + delta;
ddstar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
sec->deltaStarExt[j] = saveDStar - delta;
ddstar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
sec->deltaStarExt[j] = saveDStar;
for (int k = 0; k < ddstar[0].size(); ++k)
T.push_back(Eigen::Triplet<double>(k, i, (ddstar[0][k] - ddstar[1][k])/(2.*delta)));
+ dCd_dStar(i) = (dCdf[0] - dCdf[1])/(2.*delta);
++i;
}
Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum(dRdStar_dStar.rows(), dRdStar_dStar.cols());
@@ -663,13 +778,101 @@ void CoupledAdjoint::gradientwrtMesh(){
/**** dCk_dx ****/
std::vector<std::vector<double>> dCx_x = iadjoint->computeGradientCoefficientsMesh();
dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size());
- dCd_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
+ dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
/**** dRphi_x ****/
dRphi_x = iadjoint->getRu_X();
/**** dRx_x ****/
dRx_x = iadjoint->getRx_X();
+ /**** dRdStar_x ****/
+ std::deque<Eigen::Triplet<double>> T_dStar_x;
+ double delta = 1e-8;
+ size_t i = 0;
+ int iReg = 0;
+ double saveX = 0.;
+ std::vector<Eigen::VectorXd> ddStar(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ std::vector<double> dCdf(2, 0.);
+ for (auto sec: vsol->sections[0])
+ {
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row;
+ // x
+ saveX = sec->x[j];
+
+ sec->x[j] = saveX + delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
+ sec->x[j] = saveX - delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
+
+ sec->x[j] = saveX;
+
+ for (int k = 0; k < ddStar[0].size(); ++k)
+ T_dStar_x.push_back(Eigen::Triplet<double>(k, rowj*isol->pbl->nDim+0, -(ddStar[0][k] - ddStar[1][k])/(2.*delta)));
+ dCdf_x(rowj*isol->pbl->nDim+0) += (dCdf[0] - dCdf[1])/(2.*delta);
+ // y
+ saveX = sec->y[j];
+
+ sec->y[j] = saveX + delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[0] = this->runViscous();
+ dCdf[0] = vsol->Cdf;
+ sec->y[j] = saveX - delta;
+ sec->computeLocalCoordinates();
+ sec->xxExt = sec->loc;
+ ddStar[1] = this->runViscous();
+ dCdf[1] = vsol->Cdf;
+
+ sec->y[j] = saveX;
+
+ for (int k = 0; k < ddStar[0].size(); ++k)
+ T_dStar_x.push_back(Eigen::Triplet<double>(k, rowj*isol->pbl->nDim+1, -(ddStar[0][k] - ddStar[1][k])/(2.*delta)));
+ dCdf_x(rowj*isol->pbl->nDim+1) += (dCdf[0] - dCdf[1])/(2.*delta);
+ ++i;
+ }
+ }
+ dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end());
+ dRdStar_x.prune(0.);
+
+ /**** dRblw_x ****/
+ std::deque<Eigen::Triplet<double>> T_blw_x;
+ int offSetElms = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ if (sec->name == "wake")
+ iReg = 1;
+ else if (sec->name == "upper" || sec->name == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ std::vector<std::vector<double>> gradX = sec->evalGradwrtX();
+ std::vector<std::vector<double>> gradY = sec->evalGradwrtY();
+ for (size_t i = 0; i < gradX.size(); ++i)
+ for (size_t j = 0; j < gradX[i].size(); ++j)
+ {
+ int rowj = isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row;
+ T_blw_x.push_back(Eigen::Triplet<double>(i + offSetElms, rowj*isol->pbl->nDim+0, -gradX[i][j]));
+ T_blw_x.push_back(Eigen::Triplet<double>(i + offSetElms, rowj*isol->pbl->nDim+1, -gradY[i][j]));
+ }
+ offSetElms += sec->nElms;
+ }
+ dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end());
+ dRblw_x.prune(0.);
+
/**** dRM_x, dRrho_x, dRv_x ****/
double delta_x = 1e-8;
auto pbl = isol->pbl;
@@ -831,8 +1034,6 @@ void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_x){
Eigen::VectorXd CoupledAdjoint::runViscous()
{
int exitCode = vsol->run();
- if (exitCode != 0)
- throw std::runtime_error("CoupledAdjoint::runViscous: Solver did not converge");
// Extract deltaStar
std::vector<Eigen::VectorXd> dStar_p;
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 09936987b47e07fdb33074654f432ce15536bd29..086a70d31cd853d8508bf4d59e111bd4d8dac390 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -20,6 +20,7 @@
#include "blast.h"
#include "wObject.h"
#include "dart.h"
+#include "wTimers.h"
#include <Eigen/Sparse>
#include <iostream>
@@ -99,16 +100,23 @@ private:
// Objective functions
Eigen::VectorXd dCl_phi; ///< Derivative of the lift coefficient with respect to the inviscid flow
Eigen::VectorXd dCd_phi; ///< Derivative of the drag coefficient with respect to the inviscid flow
+ Eigen::VectorXd dCd_M; ///< Derivative of the drag coefficient with respect to the Mach number
+ Eigen::VectorXd dCd_rho; ///< Derivative of the drag coefficient with respect to the density
+ Eigen::VectorXd dCd_v; ///< Derivative of the drag coefficient with respect to the velocity
+ Eigen::VectorXd dCd_dStar; ///< Derivative of the drag coefficient with respect to delta*
// Angle of attack
Eigen::VectorXd dRphi_AoA; ///< Derivative of the inviscid flow residual with respect to the angle of attack
double dCl_AoA; ///< Partial derivative of the lift coefficient with respect to the angle of attack
- double dCd_AoA; ///< Partial erivative of the drag coefficient with respect to the angle of attack
+ double dCd_AoA; ///< Partial derivative of the total drag coefficient with respect to the angle of attack
// Mesh
Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_x; ///< Derivative of the mesh residual with respect to the mesh coordinates
Eigen::VectorXd dCl_x; ///< Partial derivative of the lift coefficient with respect to the mesh coordinates
- Eigen::VectorXd dCd_x; ///< Partial derivative of the drag coefficient with respect to the mesh coordinates
+ Eigen::VectorXd dCdp_x; ///< Partial derivative of the pressure drag coefficient with respect to the mesh coordinates
+ Eigen::VectorXd dCdf_x; ///< Partial derivative of the friction drag coefficient with respect to the mesh coordinates
+
+ fwk::Timers tms; //< internal timers
public:
CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
diff --git a/blast/tests/dart/adjointVII_2D.py b/blast/tests/dart/adjointVII_2D.py
index 44401961d30404d840627b175d422ca6e2cf1c9a..ea6db8211ea1d332f2aaf6157edf445efe66570b 100644
--- a/blast/tests/dart/adjointVII_2D.py
+++ b/blast/tests/dart/adjointVII_2D.py
@@ -18,13 +18,11 @@
# @author Paul Dechamps
# @date 2024
-# Test the blast adjoint implementation.
+# Test the blaster adjoint implementation.
# Imports.
import blast.utils as viscUtils
-from blast.interfaces.DAdjointInterface import AdjointInterface
-import modules.dartflo.dart.default as floD
import blast
import numpy as np
@@ -38,9 +36,6 @@ from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
-import os
-import sys
-
def cfgInviscid(nthrds, verb):
import os.path
# Parameters
@@ -110,24 +105,18 @@ def main():
alpha = icfg['AoA']*np.pi/180
dalpha = 1e-6 * np.pi/180
daoa = [alpha - dalpha, alpha + dalpha]
- dcl = []
- dcd = []
coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
morpher = isol.iobj['mrf']
# Adjoint objects
cAdj = blast.CoupledAdjoint(isol.adjointSolver, vsol)
- pyAdjoint = AdjointInterface(cAdj, isol, vsol)
# Run coupled case
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['forward'].start()
aeroCoeffs = coupler.run()
- tms['forward'].stop()
-
- # Print with precision of 1e-12
- print(f"{isol.solver.Cl:.12f}", f"{isol.solver.Cd:.12f}")
+ tms['forward'].stop()
print(ccolors.ANSI_BLUE + 'PyDartAdjoint...' + ccolors.ANSI_RESET)
tms['adj_dart'].start()
@@ -135,9 +124,6 @@ def main():
tms['adj_dart'].stop()
print(ccolors.ANSI_BLUE + 'PyBlasterAdjoint...' + ccolors.ANSI_RESET)
- tms['adj_py'].start()
- pyAdjoint.build()
- tms['adj_py'].stop()
tms['adj_blast'].start()
cAdj.run()
tms['adj_blast'].stop()
@@ -161,9 +147,9 @@ def main():
coupler.reset()
aeroCoeffs = coupler.run(write=False)
dcl.append(isol.getCl())
- dcd.append(isol.getCd())
- dCl_AoA_FD = (dcl[-1] - dcl[-2])/(2*dalpha)
- dCd_AoA_FD = (dcd[-1] - dcd[-2])/(2*dalpha)
+ dcd.append(isol.getCd()+vsol.Cdf)
+ dCl_AoA_FD = (dcl[-1] - dcl[-2])/(2.*dalpha)
+ dCd_AoA_FD = (dcd[-1] - dcd[-2])/(2.*dalpha)
tms['fd_blaster_aoa'].stop()
isol.setAoA(icfg['AoA']*np.pi/180) # Reset AoA after FD
@@ -190,7 +176,6 @@ def main():
import time
start_time = time.time()
for inod, n in enumerate(isol.blw[0].nodes):
- vidx = np.where(isol.vBnd[0][0].nodesCoord[:,3] == n.row)[0][0]
# Safety check
if isol.iobj['bnd'].nodes[inod].row != n.row:
raise RuntimeError('Nodes do not match', n.row, isol.iobj['bnd'].nodes[inod].row)
@@ -201,40 +186,35 @@ def main():
morpher.savePos()
n.pos[idim] = saveNodes[n.row, idim] + sign * delta
morpher.deform()
- isol.vBnd[0][0].nodesCoord[vidx, idim] = saveNodes[n.row, idim] + sign * delta
- setvMesh(isol.vBnd, isol.sec[0])
+ isol.updateMesh()
coupler.reset()
# prevent print
from contextlib import redirect_stdout
with redirect_stdout(None):
- aeroCoeffs = coupler.run(write=False)
+ coupler.run(write=False)
cmp += 1
if cmp % 20 == 0:
print('F-D opers', cmp, '/', isol.blw[0].nodes.size() * isol.getnDim() * 2, '(node '+str(n.row)+')', 'time', str(round(time.time() - start_time, 3))+'s')
cl[isgn] = isol.getCl()
- cd[isgn] = isol.getCd()
+ cd[isgn] = isol.getCd() + vsol.Cdf
# Reset mesh
for node in isol.iobj['sol'].pbl.msh.nodes:
for d in range(isol.getnDim()):
node.pos[d] = saveNodes[node.row, d]
- for jj, jrow in enumerate(isol.vBnd[0][0].nodesCoord[:,3]):
- for d in range(isol.getnDim()):
- isol.vBnd[0][0].nodesCoord[jj, d] = saveNodes[int(jrow), d]
- dCl_x_blaster_fd[inod, idim] = (cl[1] - cl[0]) / (2 * delta)
- dCd_x_blaster_fd[inod, idim] = (cd[1] - cd[0]) / (2 * delta)
+ dCl_x_blaster_fd[inod, idim] = (cl[1] - cl[0]) / (2. * delta)
+ dCd_x_blaster_fd[inod, idim] = (cd[1] - cd[0]) / (2. * delta)
tms['fd_blaster_msh'].stop()
print('Statistics')
print(tms)
-
print(ccolors.ANSI_BLUE + 'PyTesting' + ccolors.ANSI_RESET)
tests = CTests()
tests.add(CTest('dCl_dAoA', cAdj.tdCl_AoA, dCl_AoA_FD, 1e-5))
tests.add(CTest('dCd_dAoA', cAdj.tdCd_AoA, dCd_AoA_FD, 1e-4))
- tests.add(CTest('norm dCl_dx (adj - fd)', np.linalg.norm(dCl_x), np.linalg.norm(dCl_x_blaster_fd), 1e-4))
- tests.add(CTest('norm dCd_dx (adj - fd)', np.linalg.norm(dCd_x), np.linalg.norm(dCd_x_blaster_fd), 1e-3))
+ tests.add(CTest('norm dCl_dx (adj - fd)', np.linalg.norm(dCl_x), np.linalg.norm(dCl_x_blaster_fd), 1e-6))
+ tests.add(CTest('norm dCd_dx (adj - fd)', np.linalg.norm(dCd_x), np.linalg.norm(dCd_x_blaster_fd), 1e-4))
tests.add(CTest('dCl_dAoA (msh)', dCl_AoA_fromX, cAdj.tdCl_AoA, 1e-2)) # The value gets closer if the mesh if refined
tests.add(CTest('dCd_dAoA (msh)', dCd_AoA_fromX, cAdj.tdCd_AoA, 1e-1)) # The value gets closer if the mesh if refined
tests.run()
@@ -242,48 +222,5 @@ def main():
# eof
print('')
-def setvMesh(vBnd, secs):
- import tbox
- # Compute section chord
- xMin = np.min(vBnd[0][0].nodesCoord[:,0])
- xMax = np.max(vBnd[0][0].nodesCoord[:,0])
- chord = xMax - xMin
- for reg in secs:
- if reg.name == 'wake':
- iReg = 1
- elif reg.name == 'lower' or reg.name == 'upper':
- iReg = 0
- else:
- raise RuntimeError('Invalid region', reg.name)
-
- nodeRows = tbox.std_vector_int()
- for val in vBnd[0][iReg].getNodesRow(reg.name):
- nodeRows.push_back(int(val))
- connectElems = tbox.std_vector_int()
- for val in vBnd[0][iReg].getConnectList(reg.name):
- connectElems.push_back(int(val))
- # x, y, z
- xv = tbox.std_vector_double()
- for val in vBnd[0][iReg].getNodesCoord(reg.name, 0):
- xv.push_back(val)
- yv = tbox.std_vector_double()
- for val in vBnd[0][iReg].getNodesCoord(reg.name, 1):
- yv.push_back(val)
- zv = tbox.std_vector_double()
- for val in vBnd[0][iReg].getNodesCoord(reg.name, 2):
- zv.push_back(val)
-
- # Set mesh
- reg.setMesh(xv,
- yv,
- zv,
- chord,
- xMin,
- nodeRows,
- connectElems)
-
- reg.xxExt = reg.loc
-
-
if __name__ == '__main__':
main()
\ No newline at end of file