From 0346d90b9a00d84255319cce31925f1ca838d261 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Sun, 29 Sep 2024 10:44:00 +0200
Subject: [PATCH] (feat) Added viscous drag contribution to adjoint
Friction drag is now part of the implementation. The rhs of the adjoint system is updated and the test uses total drag. The adjoint test now uses the mesh update capability.
---
blast/interfaces/DAdjointInterface.py | 566 --------------------------
blast/src/DCoupledAdjoint.cpp | 229 ++++++++++-
blast/src/DCoupledAdjoint.h | 12 +-
blast/tests/dart/adjointVII_2D.py | 87 +---
4 files changed, 237 insertions(+), 657 deletions(-)
delete mode 100644 blast/interfaces/DAdjointInterface.py
diff --git a/blast/interfaces/DAdjointInterface.py b/blast/interfaces/DAdjointInterface.py
deleted file mode 100644
index 7eee2fe..0000000
--- 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 e712037..8ea9318 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 0993698..086a70d 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 4440196..ea6db82 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
--
GitLab