From 499014582b134f00d3be06272def7605dc95da54 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Mon, 23 Sep 2024 15:02:37 +0200
Subject: [PATCH] (feat) Adjoint 2D
Added adjoint for AoA and mesh in 2D. Nodes indices are now passed to the BoundaryLayer object
---
blast/coupler.py | 27 +-
blast/interfaces/DAdjointInterface.py | 612 ++++++++++++++++++---
blast/interfaces/DDataStructure.py | 34 ++
blast/interfaces/DSolversInterface.py | 65 ++-
blast/interfaces/dart/DartInterface.py | 36 +-
blast/src/DBoundaryLayer.cpp | 41 +-
blast/src/DBoundaryLayer.h | 26 +-
blast/src/DClosures.cpp | 4 +-
blast/src/DCoupledAdjoint.cpp | 706 +++++++++++++++++++++----
blast/src/DCoupledAdjoint.h | 29 +-
blast/src/DDriver.cpp | 2 +-
blast/tests/dart/adjointVII_2D.py | 181 +++++--
blast/utils.py | 2 +-
blast/validation/raeValidation.py | 8 +-
14 files changed, 1482 insertions(+), 291 deletions(-)
diff --git a/blast/coupler.py b/blast/coupler.py
index b7213a4..9c3d903 100644
--- a/blast/coupler.py
+++ b/blast/coupler.py
@@ -142,6 +142,20 @@ class Coupler:
# Inviscid solver
self.isol.reset()
+ # Viscous solver
+ for s in self.isol.sec:
+ for reg in s:
+ reg.reset()
+
+ # Blowing velocity
+ for b in self.isol.iBnd:
+ for i in range(len(b.blowingVel)):
+ b.blowingVel[i] = 0.
+ for sec in self.isol.vBnd:
+ for side in sec:
+ for i in range(len(side.blowingVel)):
+ side.blowingVel[i] = 0.
+ self.isol.setBlowingVelocity()
for iSec in range(self.isol.cfg['nSections']):
for i, reg in enumerate(self.isol.sec[iSec]):
@@ -156,13 +170,8 @@ class Coupler:
loc[iPoint] = reg.loc[iPoint]
self.isol.xxExt[iSec][i] = loc
self.isol.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
+ self.isol.vtExt[iSec][i] = np.zeros(reg.getnNodes())
+ reg.setxxExt(self.isol.xxExt[iSec][i])
+ reg.setVtExt(self.isol.vtExt[iSec][i])
+ reg.setDeltaStarExt(self.isol.deltaStarExt[iSec][i])
- # Blowing velocity
- for b in self.isol.iBnd:
- for i in range(len(b.blowingVel)):
- b.blowingVel[i] = 0.
- for sec in self.isol.vBnd:
- for side in sec:
- for i in range(len(side.blowingVel)):
- side.blowingVel[i] = 0.
- self.isol.setBlowingVelocity()
diff --git a/blast/interfaces/DAdjointInterface.py b/blast/interfaces/DAdjointInterface.py
index cf3f60f..7eee2fe 100644
--- a/blast/interfaces/DAdjointInterface.py
+++ b/blast/interfaces/DAdjointInterface.py
@@ -10,93 +10,557 @@ class AdjointInterface():
def build(self):
nDim = self.isol.solver.pbl.nDim
nNd = self.isol.solver.pbl.msh.nodes.size()
- nEs = len(self.isol.iBnd[0].elemsCoord[:,0])
- nNs = len(self.isol.blw[0].nodes)
- dRblw_M = np.zeros((nEs, nNs))
- dRblw_rho = np.zeros((nEs, nNs))
- dRblw_v = np.zeros((nEs, nDim*nNs))
- dRblw_x = np.zeros((nEs, nDim*nNd))
+
+ 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
- saveBlw = np.zeros(nEs)
- for i in range(nEs):
- saveBlw[i] = self.isol.blw[0].getU(i)
-
- # Mach
- print('Mach', end='...')
- for i, n in enumerate(self.isol.blw[0].nodes):
- savemach = self.isol.solver.M[n.row]
- self.isol.solver.M[n.row] = savemach + delta
- blowing1 = self.__runViscous()
-
- self.isol.solver.M[n.row] = savemach - delta
- blowing2 = self.__runViscous()
- self.isol.solver.M[n.row] = savemach
- dRblw_M[:,i] = -(blowing1 - blowing2)/(2*delta)
- print('done')
+ # 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)
- # Rho
- print('Density', end='...')
- for i, n in enumerate(self.isol.blw[0].nodes):
- saverho = self.isol.solver.rho[n.row]
- self.isol.solver.rho[n.row] = saverho + delta
- blowing1 = self.__runViscous()
-
- self.isol.solver.rho[n.row] = saverho - delta
- blowing2 = self.__runViscous()
- self.isol.solver.rho[n.row] = saverho
- dRblw_rho[:,i] = -(blowing1 - blowing2)/(2*delta)
- print('done')
+ 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
- # # V
- print('Velocity', end='...')
- for i, n in enumerate(self.isol.blw[0].nodes):
- for jj in range(nDim):
- savev = self.isol.solver.U[n.row][jj]
- self.isol.solver.U[n.row][jj] = savev + delta
- blowing1 = self.__runViscous()
-
- self.isol.solver.U[n.row][jj] = savev - delta
- blowing2 = self.__runViscous()
- self.isol.solver.U[n.row][jj] = savev
- dRblw_v[:,i*nDim + jj] = -(blowing1 - blowing2)/(2*delta)
+ 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')
- """# # Mesh coordinates
- print('Mesh')
- for i, n in enumerate(self.iSolverAPI.blw[0].nodes):
- vidx = np.where(self.iSolverAPI.iBnd[0].nodesCoord[:,3] == n.row)[0][0]
- for jj in range(nDim):
- if self.iSolverAPI.iBnd[0].nodesCoord[vidx,jj] != n.pos[jj]:
- print('WRONG NODE MESH')
- quit()
- savemesh = n.pos[jj]
- self.iSolverAPI.iBnd[0].nodesCoord[vidx, jj] = savemesh + delta
- blowing1 = self.__runViscous()
-
- self.iSolverAPI.iBnd[0].nodesCoord[vidx, jj] = savemesh - delta
- blowing2 = self.__runViscous()
- self.iSolverAPI.iBnd[0].nodesCoord[vidx, jj] = savemesh
- dRblw_x[:, n.row*nDim+jj] = -(blowing1 - blowing2)/(2*delta)"""
- self.coupledAdjointSolver.setdRblw_M(dRblw_M)
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):
- self.isol.getInviscidBC()
- self.isol.interpolate('i2v')
- self.isol.setViscousSolver(adjoint=True)
ext = self.vsol.run()
if ext != 0:
raise ValueError('Viscous solver did not converge')
- self.isol.getBlowingBoundary(1, adjoint=True)
- self.isol.interpolate('v2i')
- self.isol.setBlowingVelocity()
-
- blowing = np.zeros(len(self.isol.iBnd[0].elemsCoord[:,0]))
- for i in range(len(self.isol.iBnd[0].elemsCoord[:,0])):
- blowing[i] = self.isol.blw[0].getU(i)
- return blowing
+ deltaStar = self.__getDeltaStar()
+ return deltaStar
diff --git a/blast/interfaces/DDataStructure.py b/blast/interfaces/DDataStructure.py
index 4a63dc6..a5049f6 100644
--- a/blast/interfaces/DDataStructure.py
+++ b/blast/interfaces/DDataStructure.py
@@ -69,6 +69,23 @@ class Group:
else:
raise RuntimeError('Invalid region', reg)
+ def getNodesRow(self, reg:str):
+ """Return the row of the nodes
+
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ """
+ if reg == "wake":
+ return self.connectListNodes
+ elif reg == "upper":
+ return np.flip(self.connectListNodes[:self.stagPoint+1])
+ elif reg == "lower":
+ return self.connectListNodes[self.stagPoint:]
+ else:
+ raise RuntimeError('Invalid region', reg)
+
def getVelocity(self, reg:str, dim:int):
"""Return the velocity in the direction dim
@@ -140,6 +157,23 @@ class Group:
return self.Rho[self.stagPoint:]
else:
raise RuntimeError('Invalid region', reg)
+
+ def getConnectList(self, reg:str):
+ """Return the connectivity list of the elements
+
+ args:
+ ----
+ reg: str
+ Region of the boundary layer (upper, lower, wake)
+ """
+ if reg == "wake":
+ return self.connectListElems
+ elif reg == "upper":
+ return np.flip(self.connectListElems[:self.stagPoint])
+ elif reg == "lower":
+ return self.connectListElems[self.stagPoint:]
+ else:
+ raise RuntimeError('Invalid region', reg)
def setConnectList(self, _connectListNodes, _connectListElems):
"""Set connectivity lists for viscous structures.
diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index dc85a47..50eae4b 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -1,5 +1,6 @@
import numpy as np
import blast
+import tbox
from blast.interfaces.DDataStructure import Group
import time
@@ -26,6 +27,8 @@ class SolversInterface:
reg.initStructures(self.iBnd[iReg].nPoints)
self.vBnd[iSec][iReg].nodesCoord = self.iBnd[iReg].nodesCoord
self.vBnd[iSec][iReg].elemsCoord = self.iBnd[iReg].elemsCoord
+ self.vBnd[iSec][iReg].connectListNodes = self.iBnd[iReg].connectListNodes
+ self.vBnd[iSec][iReg].connectListElems = self.iBnd[iReg].connectListElems
elif self.cfg['interpolator'] == 'Rbf':
from blast.interfaces.interpolators.DRbfInterpolator import RbfInterpolator as interp
@@ -47,7 +50,7 @@ class SolversInterface:
for iSec in range(self.cfg['nSections'])]
self.vinit = False
- def setViscousSolver(self, adjoint=False):
+ def setViscousSolver(self):
"""Interpolate solution to viscous mesh and sets the inviscid boundary in the viscous solver.
"""
for iSec in range(self.cfg['nSections']):
@@ -63,11 +66,11 @@ class SolversInterface:
reg.setRhoe(self.vBnd[iSec][iReg].getRho(reg.name))
reg.setVt(self.vBnd[iSec][iReg].getVt(reg.name))
- if adjoint == False:
- reg.setDeltaStarExt(self.deltaStarExt[iSec][i])
- reg.setxxExt(self.xxExt[iSec][i])
+ reg.setDeltaStarExt(self.deltaStarExt[iSec][i])
+ reg.setxxExt(self.xxExt[iSec][i])
+ reg.setVtExt(self.vtExt[iSec][i])
- def getBlowingBoundary(self, couplIter, adjoint=False):
+ def getBlowingBoundary(self, couplIter):
"""Get blowing boundary condition from the viscous solver.
args:
----
@@ -82,11 +85,11 @@ class SolversInterface:
elif 'vAirfoil' in reg.name:
reg.blowingVel = np.concatenate((np.flip(np.asarray(self.sec[iSec][0].blowingVelocity)),
np.asarray(self.sec[iSec][1].blowingVelocity)))
-
- if adjoint == False:
- for i in range(len(self.sec[iSec])):
- self.xxExt[iSec][i] = np.asarray(self.sec[iSec][i].loc)
- self.deltaStarExt[iSec][i] = np.asarray(self.sec[iSec][i].deltaStar)
+
+ for i in range(len(self.sec[iSec])):
+ self.xxExt[iSec][i] = np.asarray(self.sec[iSec][i].loc)
+ self.deltaStarExt[iSec][i] = np.asarray(self.sec[iSec][i].deltaStar)
+ self.vtExt[iSec][i] = np.asarray(self.sec[iSec][i].vt)
def interpolate(self, dir):
if dir == 'i2v':
@@ -117,7 +120,6 @@ class SolversInterface:
name = "upper"
xtrF = self.cfg['xtrF'][j]
-
self.sec[i].append(blast.BoundaryLayer(xtrF, name))
self.vSolver.addSection(i, self.sec[i][j])
@@ -153,12 +155,41 @@ class SolversInterface:
else:
raise RuntimeError('Incorrect number of dimensions', self.getnDim())
- reg.setMesh(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx),
- self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx),
- self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx),
- chord,
- xMin)
+ nodeRows = tbox.std_vector_int()
+ if self.cfg['interpolator'] == 'Matching':
+ for val in self.vBnd[iSec][iReg].getNodesRow(reg.name):
+ nodeRows.push_back(int(val))
+ else:
+ for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))):
+ nodeRows.push_back(int(0))
+ connectElems = tbox.std_vector_int()
+ if self.cfg['interpolator'] == 'Matching':
+ for val in self.vBnd[iSec][iReg].getConnectList(reg.name):
+ connectElems.push_back(int(val))
+ else:
+ for i in range(len(self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx))-1):
+ connectElems.push_back(int(0))
+
+ # x, y, z
+ xv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, xIdx):
+ xv.push_back(val)
+ yv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, yIdx):
+ yv.push_back(val)
+ zv = tbox.std_vector_double()
+ for val in self.vBnd[iSec][iReg].getNodesCoord(reg.name, zIdx):
+ zv.push_back(val)
+
+ # Set mesh
+ reg.setMesh(xv,
+ yv,
+ zv,
+ chord,
+ xMin,
+ nodeRows,
+ connectElems)
# External variables
for i, reg in enumerate(self.sec[iSec]):
if reg.name == 'wake':
@@ -172,7 +203,9 @@ class SolversInterface:
loc[iPoint] = reg.loc[iPoint]
self.xxExt[iSec][i] = loc
self.deltaStarExt[iSec][i] = np.zeros(reg.getnNodes())
+ self.vtExt[iSec][i] = np.zeros(reg.getnNodes())
self.vinit = True
+ return 0
def getVWing(self):
"""Obtain the nodes' location and row and cg of all elements.
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index 09e06be..0334541 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -28,11 +28,11 @@ class DartInterface(SolversInterface):
def __init__(self, dartCfg, vSolver, _cfg):
try:
from modules.dartflo.dart.api.core import init_dart
- self.inviscidObjects = init_dart(dartCfg, task=dartCfg['task'], viscous=True)
- self.solver = self.inviscidObjects.get('sol') # Solver
- self.blw = [self.inviscidObjects.get('blwb'), self.inviscidObjects.get('blww')]
+ self.iobj = init_dart(dartCfg, task=dartCfg['task'], viscous=True)
+ self.solver = self.iobj.get('sol') # Solver
+ self.blw = [self.iobj.get('blwb'), self.iobj.get('blww')]
if dartCfg['task'] == 'optimization':
- self.adjointSolver = self.inviscidObjects.get('adj')
+ self.adjointSolver = self.iobj.get('adj')
print('Dart successfully loaded for optimization.')
else:
print('Dart successfully loaded for analysis.')
@@ -81,20 +81,20 @@ class DartInterface(SolversInterface):
elif self.getnDim() == 3:
# Save surface cp
- data = np.zeros((self.inviscidObjects['bnd'].nodes.size(), 4))
- for i, n in enumerate(self.inviscidObjects['bnd'].nodes):
+ data = np.zeros((self.iobj['bnd'].nodes.size(), 4))
+ for i, n in enumerate(self.iobj['bnd'].nodes):
data[i,0] = n.pos[0]
data[i,1] = n.pos[1]
data[i,2] = n.pos[2]
data[i,3] = self.solver.Cp[n.row]
- np.savetxt(self.inviscidObjects['msh'].name+'_Cp_surface'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
+ np.savetxt(self.iobj['msh'].name+'_Cp_surface'+sfx+'.dat', data, fmt='%1.6e', delimiter=', ', header='x, y, z, Cp', comments='')
import modules.dartflo.dart.utils as invUtils
if self.cfg['Format'] == 'vtk':
- print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.inviscidObjects['msh'].name, self.cfg['saveTag']))
- invUtils.write_slices(self.inviscidObjects['msh'].name, self.cfg['writeSections'], self.cfg['saveTag'], sfx=sfx)
+ print('Saving Cp files in vtk format. Msh {:s}, Tag {:.0f}'.format(self.iobj['msh'].name, self.cfg['saveTag']))
+ invUtils.write_slices(self.iobj['msh'].name, self.cfg['writeSections'], self.cfg['saveTag'], sfx=sfx)
def save(self, sfx='inviscid'):
- self.solver.save(self.inviscidObjects['wrt'], sfx)
+ self.solver.save(self.iobj['wrt'], sfx)
def getAoA(self):
return self.solver.pbl.alpha
@@ -118,7 +118,7 @@ class DartInterface(SolversInterface):
return self.solver.verbose
def getnDim(self):
- return self.inviscidObjects['pbl'].nDim
+ return self.iobj['pbl'].nDim
def getnNodesDomain(self):
return self.solver.pbl.msh.nodes.size()
@@ -266,9 +266,14 @@ class DartInterface(SolversInterface):
data[:,2] = data[connectListNodes,2]
data[:,3] = data[connectListNodes,3]
data[:,4] = data[connectListNodes,4]
+
+ connect2 = np.zeros((len(connectListElems)), dtype=int)
+ for i, e in enumerate(self.blw[0].tag.elems):
+ connect2[connectListElems[i]] = i
+
self.iBnd[0].nodesCoord = np.column_stack((data[:,1], data[:,2],\
data[:,3], data[:,4]))
- self.iBnd[0].setConnectList(connectListNodes, connectListElems)
+ self.iBnd[0].setConnectList(connectListNodes, connect2)
elemCoord = np.zeros((len(self.iBnd[0].nodesCoord)-1, 4))
for i in range(len(elemCoord[:,0])):
elemCoord[i,0] = 0.5 * (self.iBnd[0].nodesCoord[i,0] + self.iBnd[0].nodesCoord[i+1,0])
@@ -277,7 +282,6 @@ class DartInterface(SolversInterface):
elemCoord[i,3] = i
self.iBnd[0].elemsCoord = elemCoord
-
# Wake
self.iBnd[1].initStructures(self.blw[1].nodes.size())
# Node number
@@ -314,7 +318,11 @@ class DartInterface(SolversInterface):
dataW[:,4] = dataW[connectListNodes,4]
self.iBnd[1].nodesCoord = np.column_stack((dataW[:,1], dataW[:,2], \
dataW[:,3], dataW[:,4]))
- self.iBnd[1].setConnectList(connectListNodes, connectListElems)
+
+ connect2_w = np.zeros((len(connectListElems)), dtype=int)
+ for i, e in enumerate(self.blw[1].tag.elems):
+ connect2_w[connectListElems[i]] = i
+ self.iBnd[1].setConnectList(connectListNodes, connect2_w)
elemCoordW = np.zeros((len(self.iBnd[1].nodesCoord)-1, 4))
for i in range(len(elemCoordW[:,0])):
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index ad18c99..5f84a69 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -25,18 +25,51 @@ BoundaryLayer::~BoundaryLayer()
delete closSolver;
}
-void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z, double _chord, double xMin)
+void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std::vector<double> _z, double _chord, double _xMin, std::vector<int> const _rows, std::vector<int> const _no)
{
if (_x.size() != _y.size() || _x.size() != _z.size())
throw std::runtime_error("Mesh coordinates are not consistent.");
nNodes = _x.size();
nElms = nNodes - 1;
+ if (_rows.size() == 0)
+ {
+ rows.resize(0);
+ for (size_t i = 0; i < nNodes; ++i)
+ rows.push_back(static_cast<int>(i));
+ }
+ else if (_rows.size() != nNodes)
+ throw std::runtime_error("Row vector is not consistent.");
+ else
+ {
+ rows = _rows;
+ }
+ // Elements
+ if (_no.size() == 0)
+ {
+ no.resize(0);
+ for (size_t i = 0; i < nElms; ++i)
+ no.push_back(static_cast<int>(i));
+ }
+ else if (_no.size() != nElms)
+ throw std::runtime_error("No vector is not consistent.");
+ else
+ {
+ no = _no;
+ }
+
x = _x;
y = _y;
z = _z;
chord = _chord;
+ xMin = _xMin;
// Compute the local coordinate.
+ this->computeLocalCoordinates();
+ this->reset();
+}
+
+void BoundaryLayer::computeLocalCoordinates()
+{
loc.resize(nNodes, 0.);
alpha.resize(nElms, 0.);
dx.resize(nElms, 0.);
@@ -55,8 +88,6 @@ void BoundaryLayer::setMesh(std::vector<double> _x, std::vector<double> _y, std:
xoc.resize(nNodes, 0.);
for (size_t iPoint = 0; iPoint < nNodes; ++iPoint)
xoc[iPoint] = (x[iPoint] - xMin) / chord;
-
- this->reset();
}
void BoundaryLayer::reset()
@@ -94,7 +125,7 @@ void BoundaryLayer::reset()
*/
void BoundaryLayer::setStagnationBC(double Re)
{
- double dv0 = (vt[1] - vt[0]) / (loc[1] - loc[0]);
+ double dv0 = std::abs((vt[1] - vt[0]) / (loc[1] - loc[0]));
if (dv0 > 0)
u[0] = sqrt(0.075 / (Re * dv0)); // Theta
else
@@ -173,11 +204,9 @@ void BoundaryLayer::setWakeBC(double Re, std::vector<double> UpperCond, std::vec
*/
void BoundaryLayer::computeBlowingVelocity()
{
- //deltaStar[0] = u[0] * u[1];
blowingVelocity.resize(nNodes - 1, 0.);
for (size_t i = 1; i < nNodes; ++i)
{
- //deltaStar[i] = u[i * nVar + 0] * u[i * nVar + 1];
blowingVelocity[i - 1] = (rhoe[i] * vt[i] * deltaStar[i] - rhoe[i-1] * vt[i-1] * deltaStar[i - 1]) / (rhoe[i] * (loc[i] - loc[i-1]));
if (std::isnan(blowingVelocity[i - 1]))
{
diff --git a/blast/src/DBoundaryLayer.h b/blast/src/DBoundaryLayer.h
index 030bb1d..5f7adb2 100644
--- a/blast/src/DBoundaryLayer.h
+++ b/blast/src/DBoundaryLayer.h
@@ -26,16 +26,19 @@ public:
std::string name; ///< Name of the region.
// Mesh
- size_t nNodes; ///< Number of points in the domain.
- size_t nElms; ///< Number of cells in the domain.
- std::vector<double> x; ///< x coordinate in the global frame of reference.
- std::vector<double> y; ///< y coordinate in the global frame of reference.
- std::vector<double> z; ///< z coordinate in the global frame of reference.
- std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
- std::vector<double> loc; ///< xi coordinate in the local frame of reference.
- std::vector<double> dx; ///< Cell size.
- std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
- double chord; ///< Chord of the section.
+ size_t nNodes; ///< Number of points in the domain.
+ size_t nElms; ///< Number of cells in the domain.
+ std::vector<double> x; ///< x coordinate in the global frame of reference.
+ std::vector<double> y; ///< y coordinate in the global frame of reference.
+ std::vector<double> z; ///< z coordinate in the global frame of reference.
+ std::vector<double> xoc; ///< x/c coordinate in the global frame of reference.
+ std::vector<double> loc; ///< xi coordinate in the local frame of reference.
+ std::vector<int> rows; ///< Reference numbers of the nodes.
+ std::vector<int> no; ///< Reference numbers of the nodes.
+ std::vector<double> dx; ///< Cell size.
+ std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
+ double chord; ///< Chord of the section.
+ double xMin; ///< Minimum x coordinate of the mesh (for local coordinates computation).
// Boundary layer variables.
std::vector<double> u; ///< Solution vector (theta, H, N, ue, Ct).
@@ -81,7 +84,8 @@ public:
double getMaxMach() const { return *std::max_element(Me.begin(), Me.end()); };
// Setters
- void setMesh(std::vector<double> const _x, std::vector<double> const y, std::vector<double> const z, double const _chord, double const xMin);
+ void setMesh(std::vector<double> const _x, std::vector<double> const y, std::vector<double> const z, double const _chord, double const _xMin, std::vector<int> const _rows=std::vector<int>(0), std::vector<int> const _no=std::vector<int>(0));
+ void computeLocalCoordinates();
void setMe(std::vector<double> const _Me) {if (_Me.size() != nNodes) throw std::runtime_error("Mach number vector is not consistent."); Me = _Me;};
void setVt(std::vector<double> const _vt) {if (_vt.size() != nNodes) throw std::runtime_error("Velocity vector is not consistent."); vt = _vt;};
void setRhoe(std::vector<double> const _rhoe) {if (_rhoe.size() != nNodes) throw std::runtime_error("Density vector is not consistent."); rhoe = _rhoe;};
diff --git a/blast/src/DClosures.cpp b/blast/src/DClosures.cpp
index c411d1c..2cf9f44 100644
--- a/blast/src/DClosures.cpp
+++ b/blast/src/DClosures.cpp
@@ -169,8 +169,8 @@ void Closures::turbulentClosures(std::vector<double> &closureVars, double theta,
if (us > 0.95)
us = 0.98;
n = 1.0;
- cf = (1 / Fc) * (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + 0.00011 * (std::tanh(4 - (Hk / 0.875)) - 1));
- Hkc = std::max(Hk - 1 - 18 / Ret, 0.01);
+ cf = (1 / Fc) * (0.3 * std::exp(arg) * std::pow(logRt / 2.3026, (-1.74 - 0.31 * Hk)) + 0.00011 * (std::tanh(4. - (Hk / 0.875)) - 1.));
+ Hkc = std::max(Hk - 1. - 18. / Ret, 0.01);
// Dissipation coefficient.
Hmin = 1 + 2.1 / std::log(Ret);
Fl = (Hk - 1) / (Hmin - 1);
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 8139d32..e712037 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -15,6 +15,7 @@
*/
#include "DCoupledAdjoint.h"
+#include "DDriver.h"
#include "wAdjoint.h"
#include "wNewton.h"
#include "wSolver.h"
@@ -23,9 +24,12 @@
#include "wFluid.h"
#include "wMshData.h"
#include "wNode.h"
+#include "wBlowingResidual.h"
#include "wLinearSolver.h"
#include "wGmres.h"
+#include "wCache.h"
+#include "wGauss.h"
#include <Eigen/Sparse>
#include<Eigen/SparseLU>
@@ -35,6 +39,7 @@
#include <deque>
#include <algorithm>
#include <iomanip>
+#include <chrono>
using namespace blast;
@@ -58,49 +63,71 @@ CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::sh
void CoupledAdjoint::reset()
{
size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
- size_t nNs = isol->pbl->bBCs[0]->nodes.size(); // Number of surface nodes
size_t nNd = isol->pbl->msh->nodes.size(); // Number of domain nodes
- size_t nEs = isol->pbl->bBCs[0]->uE.size(); // Number of elements on the surface
+ size_t nNs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nNs += bBC->nodes.size(); // Number of surface nodes
+ nNs += 1; // Duplicate stagnation point
+ size_t nEs = 0;
+ 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);
dRv_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNd);
+ dRdStar_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
dRblw_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNd);
dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
dRM_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRdStar_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs);
dRphi_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
dRM_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
dRrho_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
dRv_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRdStar_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs);
dRphi_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim*nNs);
dRM_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
dRrho_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
dRv_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nDim*nNs);
+ dRdStar_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim*nNs);
+ dRphi_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
+ dRM_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRrho_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRv_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRdStar_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRblw_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs);
+
dRphi_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nEs);
dRM_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
dRrho_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nEs);
+ dRdStar_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nEs);
dRphi_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim*nNd);
dRM_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNd);
dRrho_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNd);
dRv_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nDim*nNd);
+ dRdStar_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNd);
dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim*nNd);
- dRM_M.setIdentity();
- dRrho_rho.setIdentity();
- dRv_v.setIdentity();
- dRblw_blw.setIdentity();
+
+ // Gradient wrt objective function
+ dRphi_AoA = Eigen::VectorXd::Zero(nNd);
// Objective functions
dCl_phi = Eigen::VectorXd::Zero(nNd);
@@ -166,21 +193,24 @@ void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMat
void CoupledAdjoint::run()
{
gradientswrtInviscidFlow();
+ gradientswrtMachNumber();
+ gradientswrtDensity();
+ gradientswrtVelocity();
+ gradientswrtDeltaStar();
gradientswrtBlowingVelocity();
gradientwrtAoA();
gradientwrtMesh();
transposeMatrices();
- // Adjoint matrix
-
// Store pointers to the matrices in a 2D vector
std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> matrices = {
- {&dRphi_phi, &dRM_phi, &dRv_phi, &dRrho_phi, &dRblw_phi},
- {&dRphi_M, &dRM_M, &dRv_M, &dRrho_M, &dRblw_M},
- {&dRphi_v, &dRM_v, &dRv_v, &dRrho_v, &dRblw_v},
- {&dRphi_rho, &dRM_rho, &dRv_rho, &dRrho_rho, &dRblw_rho},
- {&dRphi_blw, &dRM_blw, &dRv_blw, &dRrho_blw, &dRblw_blw}
+ {&dRphi_phi, &dRM_phi, &dRrho_phi, &dRv_phi, &dRdStar_phi, &dRblw_phi},
+ {&dRphi_M, &dRM_M, &dRrho_M, &dRv_M, &dRdStar_M, &dRblw_M},
+ {&dRphi_rho, &dRM_rho, &dRrho_rho, &dRv_rho, &dRdStar_rho, &dRblw_rho},
+ {&dRphi_v, &dRM_v, &dRrho_v, &dRv_v, &dRdStar_v, &dRblw_v},
+ {&dRphi_dStar, &dRM_dStar, &dRrho_dStar, &dRv_dStar, &dRdStar_dStar, &dRblw_dStar},
+ {&dRphi_blw, &dRM_blw, &dRrho_blw, &dRv_blw, &dRdStar_blw, &dRblw_blw}
};
int rows = 0; int cols = 0;
@@ -189,59 +219,87 @@ void CoupledAdjoint::run()
Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols);
buildAdjointMatrix(matrices, A_adjoint);
+ size_t phiIdx = 0;
+ size_t machIdx = dRphi_phi.cols();
+ size_t rhoIdx = machIdx + dRM_M.cols();
+ size_t vIdx = rhoIdx + dRrho_rho.cols();
+ size_t dStarIdx = vIdx + dRv_v.cols();
+ size_t blwIdx = dStarIdx + dRdStar_dStar.cols();
+
//***** Gradient wrt angle of attack *****//
//****************************************//
// CL
std::vector<double> rhsCl(A_adjoint.cols(), 0.0);
for (auto i = 0; i<dCl_phi.size(); ++i)
- rhsCl[i] = dCl_phi[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_);
- Eigen::VectorXd lambdaClPhi = lambdaCl.segment(0, isol->pbl->msh->nodes.size());
- Eigen::VectorXd lambdaClBlw = lambdaCl.segment(lambdaCl.size() - isol->pbl->bBCs[0]->uE.size(), isol->pbl->bBCs[0]->uE.size());
// Total gradient
- tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaClPhi;
+ tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaCl.segment(phiIdx, dRphi_phi.cols());
// CD
- // std::vector<double> rhsCd(A_adjoint.cols(), 0.0);
- // for (auto i = 0; i<dCd_phi.size(); ++i)
- // rhsCd[i] = dCd_phi[i];
- // 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());
- // alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCd_.data(), rhsCd_.size()), lambdaCd_);
- // Eigen::VectorXd lambdaCdPhi = lambdaCd.block(0, 0, isol->pbl->msh->nodes.size(), 1);
-
- // tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCdPhi; // Total gradient
+ std::vector<double> rhsCd(A_adjoint.cols(), 0.0);
+ for (auto i = 0; i<dCd_phi.size(); ++i)
+ rhsCd[i] = dCd_phi[i];
+ 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
//***** Gradient wrt mesh coordinates *****//
//*****************************************//
- // CL
- // Solution of (from augmented Lagrangian, eqn d/dx )
- // "dCl_x + dRphi_x * lambdaCl_phi + dRblw_x * lambdaCl_blw + dRx_x * lambdaCl_x"
+ // Ck (cL, cd)
+ // 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"
+
+ 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
+ rhsCl_x -= dRrho_x.transpose() * lambdaCl.segment(rhoIdx, dRrho_rho.cols()); // Density contribution
+ rhsCl_x -= dRv_x.transpose() * lambdaCl.segment(vIdx, dRv_v.cols()); // Velocity contribution
+ rhsCl_x -= dRdStar_x.transpose() * lambdaCl.segment(dStarIdx, dRdStar_dStar.cols()); // Displacement thickness contribution
+ rhsCl_x -= dRblw_x.transpose() * lambdaCl.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution
+
Eigen::VectorXd lambdaCl_x = Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size());
Eigen::Map<Eigen::VectorXd> lambdaCl_x_(lambdaCl_x.data(), lambdaCl_x.size());
-
- Eigen::VectorXd rhsCl_x = dCl_x - dRphi_x.transpose() * lambdaClPhi - dRblw_x.transpose() * lambdaClBlw;
- dRx_x.setIdentity();
alinsol->compute(dRx_x.transpose(), Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size()), lambdaCl_x_);
+ Eigen::VectorXd rhsCd_x = dCd_x;
+ 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
+ rhsCd_x -= dRv_x.transpose() * lambdaCd.segment(vIdx, dRv_v.cols()); // Velocity contribution
+ rhsCd_x -= dRdStar_x.transpose() * lambdaCd.segment(dStarIdx, dRdStar_dStar.cols()); // Displacement thickness contribution
+ rhsCd_x -= dRblw_x.transpose() * lambdaCd.segment(blwIdx, dRblw_blw.cols()); // Blowing contribution
+
+ Eigen::VectorXd lambdaCd_x = Eigen::VectorXd::Zero(isol->pbl->nDim * isol->pbl->msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> lambdaCd_x_(lambdaCd_x.data(), lambdaCd_x.size());
+ alinsol->compute(dRx_x.transpose(), Eigen::Map<Eigen::VectorXd>(rhsCd_x.data(), rhsCd_x.size()), lambdaCd_x_);
+
for (auto n : isol->pbl->msh->nodes)
{
for (int m = 0; m < isol->pbl->nDim; m++)
{
tdCl_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m];
- // TODO: CHANGE FOR CD
- tdCd_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m];
+ tdCd_x[n->row](m) = lambdaCd_x[isol->pbl->nDim * (n->row) + m];
}
}
}
+/**
+ * @brief Compute all the gradients wrt the inviscid flow (phi) in the domain
+ * Non-zero are: dRphi_phi, dRM_phi, dRrho_phi, dRv_phi
+ * Zeros are: dRdStar_phi, dRblw_phi, dRx_phi
+ */
void CoupledAdjoint::gradientswrtInviscidFlow()
{
//**** dRphi_phi ****//
@@ -251,7 +309,7 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
auto pbl = isol->pbl;
// Finite difference
std::vector<double> Phi_save = isol->phi; // Differential of the independent variable (phi)
- double deltaPhi = 1e-6; // perturbation
+ double deltaPhi = 1.e-8; // perturbation
std::vector<Eigen::Triplet<double>> tripletsdM;
std::vector<Eigen::Triplet<double>> tripletsdrho;
@@ -261,30 +319,39 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
Phi_save[n->row] = isol->phi[n->row] + deltaPhi;
// Recompute the quantities on surface nodes.
- int jj = 0;
- for (auto srfNode : pbl->bBCs[0]->nodes){
-
- double dM = 0.; double drho = 0.;
- Eigen::Vector3d dv = Eigen::Vector3d::Zero();
- // Average of the quantity on the elements adjacent to the considered node.
- for (auto e : pbl->fluid->neMap[srfNode]){
- dM += pbl->fluid->mach->eval(*e, Phi_save, 0);
- drho += pbl->fluid->rho->eval(*e, Phi_save, 0);
- Eigen::VectorXd grad = e->computeGradient(Phi_save, 0);
- for (int i = 0; i < grad.size(); ++i)
- dv(i) += grad(i);
+ size_t jj = 0;
+ for (auto sec: vsol->sections[0]){
+ int iReg = 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]];
+ double dM = 0.; double drho = 0.;
+ Eigen::Vector3d dv = Eigen::Vector3d::Zero();
+ // Average of the quantity on the elements adjacent to the considered node.
+ for (auto e : pbl->fluid->neMap[srfNode]){
+ dM += pbl->fluid->mach->eval(*e, Phi_save, 0);
+ drho += pbl->fluid->rho->eval(*e, Phi_save, 0);
+ Eigen::VectorXd grad = e->computeGradient(Phi_save, 0);
+ for (int i = 0; i < grad.size(); ++i)
+ dv(i) += grad(i);
+ }
+ dM /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ drho /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ dv /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+
+ // Form triplets
+ // Minus because M = f(phi) and therefore RM = M - f(phi) = 0
+ tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row, -(dM - isol->M[srfNode->row])/deltaPhi));
+ tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row, -(drho - isol->rho[srfNode->row])/deltaPhi));
+ for (int i = 0; i < pbl->nDim; ++i)
+ tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row, -(dv(i) - isol->U[srfNode->row](i))/deltaPhi));
+ jj++;
}
- dM /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
- drho /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
- dv /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
-
- // Form triplets
- tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row, -(dM - isol->M[srfNode->row])/deltaPhi));
- tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row, -(drho - isol->rho[srfNode->row])/deltaPhi));
- for (int i = 0; i < pbl->nDim; ++i)
- tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row, -(dv(i) - isol->U[srfNode->row](i))/deltaPhi));
-
- jj++;
}
// Reset phi
Phi_save[n->row] = isol->phi[n->row];
@@ -296,111 +363,464 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
// Remove zeros
dRM_phi.prune(0.); dRrho_phi.prune(0.); dRv_phi.prune(0.);
- //**** dCl_phi, dCd_phi ****//
+ //**** 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());
dCd_phi = Eigen::VectorXd::Map(dCx_phi[1].data(), dCx_phi[1].size());
}
-void CoupledAdjoint::gradientswrtMachNumber(){
- // dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->msh->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
- // dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
- // dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->nDim*isol->pbl->bBCs[0]->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
- // dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->uE.size(), isol->pbl->bBCs[0]->nodes.size());
+/**
+ * @brief Compute all the gradients wrt the Mach number on the surface
+ * Non-zero are: dRM_M, dRdStar_M
+ * Zeros are: dRphi_M, dRrho_M, dRv_M, dRblw_M, dRx_M
+ */
+void CoupledAdjoint::gradientswrtMachNumber()
+{
+ /**** dRM_M ****/
+ dRM_M.setIdentity();
+
+ //**** dRdStar_M ****//
+ std::deque<Eigen::Triplet<double>> T;
+ double delta = 1e-6;
+ size_t i = 0;
+ double saveM = 0.;
+ std::vector<Eigen::VectorXd> ddStar(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ for (auto sec: vsol->sections[0])
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ saveM = sec->Me[j];
+
+ sec->Me[j] = saveM + delta;
+ ddStar[0] = this->runViscous();
+ sec->Me[j] = saveM - delta;
+ ddStar[1] = this->runViscous();
+
+ 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)));
+ ++i;
+ }
+ dRdStar_M.setFromTriplets(T.begin(), T.end());
+ dRdStar_M.prune(0.);
}
-void CoupledAdjoint::gradientswrtDensity(){
- // dRphi_dRho = 0
- // dRM_dRho = 0
- // dRrho_dRho =
- // dRv_dRho = 0
- // dRblw_dRho =
+/**
+ * @brief Compute all the gradients wrt the density on the surface
+ * Non-zero are: dRrho_rho, dRblw_rho
+ * Zeros are: dRphi_rho, dRM_rho, dRv_rho, dRdStar_rho, dRx_rho
+ */
+void CoupledAdjoint::gradientswrtDensity()
+{
+ size_t nNs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nNs += bBC->nodes.size(); // Number of surface nodes
+ nNs += 1; // Duplicate stagnation point
+ size_t nEs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nEs += bBC->uE.size(); // Number of blowing elements
+
+ //**** dRrho_rho ****//
+ dRrho_rho.setIdentity();
- // dCl_dRho =
- // dCd_dRho =
+ //**** dRblw_rho ****//
+ std::deque<Eigen::Triplet<double>> T_blw;
+ int offSetElms = 0;
+ int offSetNodes = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ std::vector<std::vector<double>> grad = sec->evalGradwrtRho();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, -grad[i][j]));
+ offSetElms += sec->nElms;
+ offSetNodes += sec->nNodes;
+ }
+ dRblw_rho.setFromTriplets(T_blw.begin(), T_blw.end());
+ dRblw_rho.prune(0.);
}
-void CoupledAdjoint::gradientswrtVelocity(){
- // dRphi_dV = 0
- // dRM_dV = 0
- // dRrho_dV = 0
- // dRv_dV =
- // dRblw_dV =
+/**
+ * @brief Compute all the gradients wrt the velocity on the surface
+ * Non-zeros are: dRv_v, dRdStar_v, dRblw_v
+ * Zeros are: dRphi_v, dRM_v, dRrho_v, dRx_v
+ */
+void CoupledAdjoint::gradientswrtVelocity()
+{
+ size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
+ size_t nNs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nNs += bBC->nodes.size(); // Number of surface nodes
+ nNs += 1; // Duplicate stagnation point
+ size_t nEs = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nEs += bBC->uE.size(); // Number of blowing elements
+
+ /**** dRv_v ****/
+ dRv_v.setIdentity();
+
+ //**** Get velocity****/
+ Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nNs, nDim);
+ size_t i = 0;
+ int iReg = -1;
+ 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");
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ for (size_t idim = 0; idim < nDim; ++idim)
+ v(i, idim) = isol->U[isol->pbl->bBCs[iReg]->nodes[sec->rows[j]]->row](idim);
+ ++i;
+ }
+ }
+
+ //**** dvt_v (analytical) ****//
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim);
+ std::deque<Eigen::Triplet<double>> T_vt;
+
+ i = 0;
+ for (const auto& sec : vsol->sections[0]) {
+ for (size_t j = 0; j < sec->nNodes; ++j) {
+ for (size_t idim = 0; idim < nDim; ++idim) {
+ T_vt.push_back(Eigen::Triplet<double>(i, i * nDim + idim, 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim)));
+ }
+ ++i;
+ }
+ }
+ dVt_v.setFromTriplets(T_vt.begin(), T_vt.end());
+ dVt_v.prune(0.);
+
+ //**** dRdStar_vt (finite-difference) ****//
+ std::deque<Eigen::Triplet<double>> T;
+ double delta = 1e-6;
+ i = 0;
+ double saveM = 0.;
+ std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
+ for (auto sec: vsol->sections[0])
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ saveM = sec->vt[j];
+
+ sec->vt[j] = saveM + delta;
+ dvt[0] = this->runViscous();
+ sec->vt[j] = saveM - delta;
+ dvt[1] = this->runViscous();
- // dCl_dV =
- // dCd_dV =
+ 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)));
+ ++i;
+ }
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(), dRdStar_M.cols());
+ dRdStar_vt.setFromTriplets(T.begin(), T.end());
+ dRdStar_vt.prune(0.);
+
+ //**** dRdStar_v ****//
+ dRdStar_v = dRdStar_vt * dVt_v;
+ dRdStar_v.prune(0.);
+
+
+ //**** dRblw_vt (analytical) ****//
+ std::deque<Eigen::Triplet<double>> T_blw;
+ int offSetElms = 0;
+ int offSetNodes = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ std::vector<std::vector<double>> grad = sec->evalGradwrtVt();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, -grad[i][j]));
+ offSetElms += sec->nElms;
+ offSetNodes += sec->nNodes;
+ }
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt(nEs, nNs);
+ dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end());
+ dRblw_vt.prune(0.);
+
+ //**** dRblw_v ****//
+ dRblw_v = dRblw_vt * dVt_v;
+ dRblw_v.prune(0.);
}
-void CoupledAdjoint::gradientswrtBlowingVelocity(){
- dRphi_blw = iadjoint->getRu_Blw();
+/**
+ * @brief Compute all the gradients wrt the displacement thickness on the surface
+ * Non-zero are: dRdStar_dStar, dRblw_dStar
+ * Zeros are: dRphi_dStar, dRM_dStar, dRrho_dStar, dRv_dStar, dRx_dStar
+ */
+void CoupledAdjoint::gradientswrtDeltaStar()
+{
+ //**** dRdStar_dStar (finite-difference) ****//
+ std::deque<Eigen::Triplet<double>> T;
+ double delta = 1e-6;
+ double saveDStar = 0.;
+ std::vector<Eigen::VectorXd> ddstar(2, Eigen::VectorXd::Zero(dRdStar_dStar.cols()));
+ size_t i = 0;
+ for (auto sec: vsol->sections[0])
+ for (size_t j = 0; j < sec->nNodes; ++j)
+ {
+ saveDStar = sec->deltaStarExt[j];
- // All others are zero.
+ sec->deltaStarExt[j] = saveDStar + delta;
+ ddstar[0] = this->runViscous();
+ sec->deltaStarExt[j] = saveDStar - delta;
+ ddstar[1] = this->runViscous();
+
+ 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)));
+ ++i;
+ }
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum(dRdStar_dStar.rows(), dRdStar_dStar.cols());
+ // RdStar = dStar - f(dStar) = 0
+ // dRdStar_dStar = 1 - df_dStar (minus is 3 lines below and not in the triplets)
+ dRdStar_dStar_dum.setFromTriplets(T.begin(), T.end());
+ dRdStar_dStar.setIdentity();
+ dRdStar_dStar -= dRdStar_dStar_dum;
+ dRdStar_dStar.prune(0.);
+
+ //**** dRblw_dStar (analytical) ****//
+ std::deque<Eigen::Triplet<double>> T_blw;
+ int offSetElms = 0;
+ int offSetNodes = 0;
+ for (const auto& sec: vsol->sections[0])
+ {
+ std::vector<std::vector<double>> grad = sec->evalGradwrtDeltaStar();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes, -grad[i][j]));
+ offSetElms += sec->nElms;
+ offSetNodes += sec->nNodes;
+ }
+ dRblw_dStar.setFromTriplets(T_blw.begin(), T_blw.end());
+ dRblw_dStar.prune(0.);
}
+/**
+ * @brief Compute all the gradients wrt the blowing velocity
+ * Non-zero are: dRphi_blw, dRblw_blw
+ * Zeros are: dRM_blw, dRrho_blw, dRv_blw, dRdStar_blw, dRx_blw
+ */
+void CoupledAdjoint::gradientswrtBlowingVelocity()
+{
+ /**** dRphi_blw ****/
+ std::deque<Eigen::Triplet<double>> T;
+
+ 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("Unknown section name");
+ for (size_t iElm = 0; iElm < sec->nElms; ++iElm){
+ auto blowingElement = isol->pbl->bBCs[iReg]->uE[sec->no[iElm]].first;
+ Eigen::VectorXd be = dart::BlowingResidual::buildGradientBlowing(*blowingElement);
+ for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++)
+ {
+ tbox::Node *nodi = blowingElement->nodes[ii];
+ T.push_back(Eigen::Triplet<double>(isol->rows[nodi->row], jj, be(ii)));
+ }
+ ++jj;
+ }
+ }
+ dRphi_blw.setFromTriplets(T.begin(), T.end());
+ dRphi_blw.prune(0.);
+ dRphi_blw.makeCompressed();
+
+ /**** dRblw_blw ****/
+ dRblw_blw.setIdentity();
+}
+
+/**
+ * @brief Compute all the gradients wrt the angle of attack
+ * Non-zero are: dRphi_AoA
+ * Zeros are: dRM_AoA, dRrho_AoA, dRv_AoA, dRdStar_AoA, dRblw_AoA, dRx_AoA
+ */
void CoupledAdjoint::gradientwrtAoA(){
std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa();
dCl_AoA = dCx_AoA[0]; dCd_AoA = dCx_AoA[1];
dRphi_AoA = iadjoint->getRu_A();
-
- // All others are zero.
}
+/**
+ * @brief Compute all the gradients wrt the mesh coordinates
+ * Non-zero are: dRphi_x, dRM_x, dRrho_x, dRv_x, dRdStar_x, dRblw_x, dRx_x
+ * Zeros are: -
+ */
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());
+ /**** dRphi_x ****/
dRphi_x = iadjoint->getRu_X();
+ /**** dRx_x ****/
dRx_x = iadjoint->getRx_X();
+
+ /**** dRM_x, dRrho_x, dRv_x ****/
+ double delta_x = 1e-8;
+ auto pbl = isol->pbl;
+ std::vector<Eigen::Triplet<double>> tripletsdM;
+ std::vector<Eigen::Triplet<double>> tripletsdrho;
+ std::vector<Eigen::Triplet<double>> tripletsdv;
+
+ for (auto n : pbl->msh->nodes){
+ for (int idim = 0; idim < isol->pbl->nDim; ++idim){
+ // Modify node position
+ double save = n->pos[idim];
+ n->pos[idim] = save + delta_x;
+ for (auto e: pbl->msh->elems)
+ if (e->type() == tbox::ElType::LINE2 || e->type() == tbox::ElType::TRI3) e->update();
+
+ // Recompute the quantities on surface nodes.
+ size_t jj = 0;
+ for (auto sec: vsol->sections[0]){
+ int iReg = 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]];
+ double dM = 0.; double drho = 0.;
+ Eigen::Vector3d dv = Eigen::Vector3d::Zero();
+ // Average of the quantity on the elements adjacent to the considered node.
+ for (auto e : pbl->fluid->neMap[srfNode]){
+ dM += pbl->fluid->mach->eval(*e, isol->phi, 0);
+ drho += pbl->fluid->rho->eval(*e, isol->phi, 0);
+ Eigen::VectorXd grad = e->computeGradient(isol->phi, 0);
+ for (int i = 0; i < grad.size(); ++i)
+ dv(i) += grad(i);
+ }
+ dM /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ drho /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ dv /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+
+ // Form triplets
+ // Minus because M = f(phi) and therefore RM = M - f(phi) = 0
+ tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row*pbl->nDim + idim, -(dM - isol->M[srfNode->row])/delta_x));
+ tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row*pbl->nDim + idim, -(drho - isol->rho[srfNode->row])/delta_x));
+ for (int i = 0; i < pbl->nDim; ++i)
+ tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row*pbl->nDim + idim, -(dv(i) - isol->U[srfNode->row](i))/delta_x));
+ jj++;
+ }
+ }
+ // Reset node pos
+ n->pos[idim] = save;
+ }
+ }
+ // Fill matrices with gradients
+ dRM_x.setFromTriplets(tripletsdM.begin(), tripletsdM.end());
+ dRrho_x.setFromTriplets(tripletsdrho.begin(), tripletsdrho.end());
+ dRv_x.setFromTriplets(tripletsdv.begin(), tripletsdv.end());
+ // Remove zeros
+ dRM_x.prune(0.); dRrho_x.prune(0.); dRv_x.prune(0.);
}
-void CoupledAdjoint::setdRblw_M(std::vector<std::vector<double>> _dRblw_dM){
- //dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dM.size(), _dRblw_dM[0].size());
- // Convert std::vector<double> to Eigen::Triplets
+void CoupledAdjoint::setdRdStar_M(std::vector<std::vector<double>> _dRdStar_dM){
std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_dM.size(); i++){
- for (size_t j = 0; j < _dRblw_dM[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dM[i][j]));
+ for (size_t i = 0; i < _dRdStar_dM.size(); i++){
+ for (size_t j = 0; j < _dRdStar_dM[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dM[i][j]));
}
}
// Set matrix
- dRblw_M.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_M.prune(0.);
+ dRdStar_M.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_M.prune(0.);
}
-void CoupledAdjoint::setdRblw_v(std::vector<std::vector<double>> _dRblw_dv){
- //dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dv.size(), _dRblw_dv[0].size());
- // Convert std::vector<double> to Eigen::Triplets
+
+void CoupledAdjoint::setdRdStar_v(std::vector<std::vector<double>> _dRdStar_dv){
std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_dv.size(); i++){
- for (size_t j = 0; j < _dRblw_dv[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dv[i][j]));
+ for (size_t i = 0; i < _dRdStar_dv.size(); i++){
+ for (size_t j = 0; j < _dRdStar_dv[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dv[i][j]));
}
}
// Set matrix
- dRblw_v.setFromTriplets(triplets.begin(), triplets.end());
- dRblw_v.prune(0.);
+ dRdStar_v.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_v.prune(0.);
}
-void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_drho){
- //dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_drho.size(), _dRblw_drho[0].size());
- // Convert std::vector<double> to Eigen::Triplets
+
+void CoupledAdjoint::setdRdStar_dStar(std::vector<std::vector<double>> _dRdStar_dStar){
std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_drho.size(); i++){
- for (size_t j = 0; j < _dRblw_drho[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_drho[i][j]));
+ for (size_t i = 0; i < _dRdStar_dStar.size(); i++){
+ for (size_t j = 0; j < _dRdStar_dStar[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_dStar[i][j]));
+ }
+ }
+ // Set matrix
+ dRdStar_dStar.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_dStar.prune(0.);
+}
+
+void CoupledAdjoint::setdRdStar_x(std::vector<std::vector<double>> _dRdStar_x){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRdStar_x.size(); i++){
+ for (size_t j = 0; j < _dRdStar_x[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRdStar_x[i][j]));
+ }
+ }
+ // Set matrix
+ dRdStar_x.setFromTriplets(triplets.begin(), triplets.end());
+ dRdStar_x.prune(0.);
+}
+
+void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_rho){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_rho.size(); i++){
+ for (size_t j = 0; j < _dRblw_rho[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_rho[i][j]));
}
}
// Set matrix
dRblw_rho.setFromTriplets(triplets.begin(), triplets.end());
dRblw_rho.prune(0.);
}
-void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_dx){
- //dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dx.size(), _dRblw_dx[0].size());
- // Convert std::vector<double> to Eigen::Triplets
+
+void CoupledAdjoint::setdRblw_v(std::vector<std::vector<double>> _dRblw_v){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_v.size(); i++){
+ for (size_t j = 0; j < _dRblw_v[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_v[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_v.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_v.prune(0.);
+}
+
+void CoupledAdjoint::setdRblw_dStar(std::vector<std::vector<double>> _dRblw_dStar){
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_dStar.size(); i++){
+ for (size_t j = 0; j < _dRblw_dStar[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dStar[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_dStar.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_dStar.prune(0.);
+}
+
+void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_x){
std::vector<Eigen::Triplet<double>> triplets;
- for (size_t i = 0; i < _dRblw_dx.size(); i++){
- for (size_t j = 0; j < _dRblw_dx[i].size(); j++){
- triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dx[i][j]));
+ for (size_t i = 0; i < _dRblw_x.size(); i++){
+ for (size_t j = 0; j < _dRblw_x[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_x[i][j]));
}
}
// Set matrix
@@ -408,47 +828,101 @@ void CoupledAdjoint::setdRblw_x(std::vector<std::vector<double>> _dRblw_dx){
dRblw_x.prune(0.);
}
+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;
+
+ // Extract parts in a loop
+ for (size_t i = 0; i < vsol->sections[0].size(); ++i) {
+ std::vector<double> deltaStar_vec = this->vsol->sections[0][i]->deltaStar;
+ Eigen::VectorXd deltaStar_eigen = Eigen::Map<Eigen::VectorXd>(deltaStar_vec.data(), deltaStar_vec.size());
+ dStar_p.push_back(deltaStar_eigen);
+ }
+ // Calculate the total size
+ int nNs = 0;
+ for (const auto& p : dStar_p)
+ nNs += p.size();
+
+ // Concatenate the vectors
+ Eigen::VectorXd deltaStar(nNs);
+ int idx = 0;
+ for (const auto& p : dStar_p) {
+ deltaStar.segment(idx, p.size()) = p;
+ idx += p.size();
+ }
+ return deltaStar;
+}
+
/**
* @brief Transpose all matrices of the adjoint problem included in the system. Not the others.
*/
void CoupledAdjoint::transposeMatrices(){
+
dRphi_phi = dRphi_phi.transpose();
dRphi_M = dRphi_M.transpose();
dRphi_rho = dRphi_rho.transpose();
dRphi_v = dRphi_v.transpose();
+ dRphi_dStar = dRphi_dStar.transpose();
dRphi_blw = dRphi_blw.transpose();
dRM_phi = dRM_phi.transpose();
dRM_M = dRM_M.transpose();
dRM_rho = dRM_rho.transpose();
dRM_v = dRM_v.transpose();
+ dRM_dStar = dRM_dStar.transpose();
dRM_blw = dRM_blw.transpose();
dRrho_phi = dRrho_phi.transpose();
dRrho_M = dRrho_M.transpose();
dRrho_rho = dRrho_rho.transpose();
dRrho_v = dRrho_v.transpose();
+ dRrho_dStar = dRrho_dStar.transpose();
dRrho_blw = dRrho_blw.transpose();
dRv_phi = dRv_phi.transpose();
dRv_M = dRv_M.transpose();
dRv_rho = dRv_rho.transpose();
dRv_v = dRv_v.transpose();
+ dRv_dStar = dRv_dStar.transpose();
dRv_blw = dRv_blw.transpose();
+ dRdStar_phi = dRdStar_phi.transpose();
+ dRdStar_M = dRdStar_M.transpose();
+ dRdStar_rho = dRdStar_rho.transpose();
+ dRdStar_v = dRdStar_v.transpose();
+ dRdStar_dStar = dRdStar_dStar.transpose();
+ dRdStar_blw = dRdStar_blw.transpose();
+
dRblw_phi = dRblw_phi.transpose();
dRblw_M = dRblw_M.transpose();
dRblw_rho = dRblw_rho.transpose();
dRblw_v = dRblw_v.transpose();
+ dRblw_dStar = dRblw_dStar.transpose();
dRblw_blw = dRblw_blw.transpose();
+
}
void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename) {
std::ofstream file(filename);
if (file.is_open()) {
+ // Write the column headers
+ for (int j = 0; j < matrix.cols(); ++j) {
+ file << std::setw(15) << j;
+ if (j != matrix.cols() - 1) {
+ file << " ";
+ }
+ }
+ file << "\n";
+
+ // Write the matrix data
for (int i = 0; i < matrix.rows(); ++i) {
for (int j = 0; j < matrix.cols(); ++j) {
- file << std::setprecision(20) << matrix(i, j);
+ file << std::fixed << std::setprecision(12) << std::setw(15) << matrix(i, j);
if (j != matrix.cols() - 1) {
file << " ";
}
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 3ae603d..0993698 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -52,6 +52,7 @@ private:
Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_M; ///< Derivative of the inviscid flow residual with respect to the Mach number
Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_v; ///< Derivative of the inviscid flow residual with respect to the velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_rho; ///< Derivative of the inviscid flow residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_dStar; ///< Derivative of the inviscid flow residual with respect to delta*
Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_blw; ///< Derivative of the inviscid flow residual with respect to the blowing velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_x; ///< Derivative of the potential residual with respect to the mesh coordinates
@@ -59,6 +60,7 @@ private:
Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_M; ///< Derivative of the Mach number residual with respect to the Mach number
Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_v; ///< Derivative of the Mach number residual with respect to the velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_rho; ///< Derivative of the Mach number residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_dStar; ///< Derivative of the Mach number residual with respect to delta*
Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_blw; ///< Derivative of the Mach number residual with respect to the blowing velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_x; ///< Derivative of the Mach number residual with respect to the mesh coordinates
@@ -66,6 +68,7 @@ private:
Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_M; ///< Derivative of the density residual with respect to the Mach number
Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_v; ///< Derivative of the density residual with respect to the velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_rho; ///< Derivative of the density residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_dStar; ///< Derivative of the density residual with respect to delta*
Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_blw; ///< Derivative of the density residual with respect to the blowing velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_x; ///< Derivative of the density residual with respect to the mesh coordinates
@@ -73,13 +76,23 @@ private:
Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_M; ///< Derivative of the velocity residual with respect to the Mach number
Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_v; ///< Derivative of the velocity residual with respect to the velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_rho; ///< Derivative of the velocity residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_dStar; ///< Derivative of the velocity residual with respect to delta*
Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_blw; ///< Derivative of the velocity residual with respect to the blowing velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_x; ///< Derivative of the velocity residual with respect to the mesh coordinates
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_phi; ///< Derivative of delta* residual with respect to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_M; ///< Derivative of delta* residual with respect to the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_v; ///< Derivative of delta* residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_rho; ///< Derivative of delta* residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar; ///< Derivative of delta* residual with respect to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_blw; ///< Derivative of delta* residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_x; ///< Derivative of delta* residual with respect to the mesh coordinates
+
Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_phi; ///< Derivative of the blowing velocity residual with respect to the inviscid flow
Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_M; ///< Derivative of the blowing velocity residual with respect to the Mach number
Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_v; ///< Derivative of the blowing velocity residual with respect to the velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_rho; ///< Derivative of the blowing velocity residual with respect to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_dStar; ///< Derivative of the blowing velocity residual with respect to delta*
Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_blw; ///< Derivative of the blowing velocity residual with respect to the blowing velocity
Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_x; ///< Derivative of the blowing velocity residual with respect to the mesh coordinates
@@ -104,18 +117,26 @@ public:
void run();
void reset();
- void setdRblw_M(std::vector<std::vector<double>> dRblw_dM);
- void setdRblw_v(std::vector<std::vector<double>> dRblw_dv);
- void setdRblw_rho(std::vector<std::vector<double>> dRblw_drho);
- void setdRblw_x(std::vector<std::vector<double>> dRblw_dx);
+ void setdRdStar_M(std::vector<std::vector<double>> _dRdStar_M);
+ void setdRdStar_v(std::vector<std::vector<double>> _dRdStar_v);
+ void setdRdStar_dStar(std::vector<std::vector<double>> _dRdStar_dStar);
+ void setdRdStar_x(std::vector<std::vector<double>> _dRdStar_x);
+
+ void setdRblw_rho(std::vector<std::vector<double>> _dRblw_rho);
+ void setdRblw_v(std::vector<std::vector<double>> _dRblw_v);
+ void setdRblw_dStar(std::vector<std::vector<double>> _dRblw_ddStar);
+ void setdRblw_x(std::vector<std::vector<double>> _dRblw_dx);
private:
void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
+ Eigen::VectorXd runViscous();
+
void gradientswrtInviscidFlow();
void gradientswrtMachNumber();
void gradientswrtDensity();
void gradientswrtVelocity();
+ void gradientswrtDeltaStar();
void gradientswrtBlowingVelocity();
void gradientwrtAoA();
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 2fee439..8a4fb03 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -266,7 +266,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, bool forced)
double avTurb = (bl->x[iPoint] - bl->xtr) / (bl->x[iPoint] - bl->x[iPoint-1]);
double avLam = 1. - avTurb;
- //std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
+ // std::cout << "avTurb " << avTurb << " avLam " << avLam << std::endl;
if (avTurb < 0. || avTurb > 1. || avLam < 0. || avLam > 1.)
{
bl->printSolution(iPoint);
diff --git a/blast/tests/dart/adjointVII_2D.py b/blast/tests/dart/adjointVII_2D.py
index 7341e98..4440196 100644
--- a/blast/tests/dart/adjointVII_2D.py
+++ b/blast/tests/dart/adjointVII_2D.py
@@ -24,6 +24,7 @@
import blast.utils as viscUtils
from blast.interfaces.DAdjointInterface import AdjointInterface
+import modules.dartflo.dart.default as floD
import blast
import numpy as np
@@ -37,7 +38,8 @@ from fwk.testing import *
from fwk.coloring import ccolors
import numpy as np
-import math
+import os
+import sys
def cfgInviscid(nthrds, verb):
import os.path
@@ -50,7 +52,7 @@ def cfgInviscid(nthrds, verb):
'Verb' : verb, # verbosity
# Model (geometry or mesh)
'File' : os.path.dirname(os.path.abspath(__file__)) + '/../../models/dart/n0012.geo', # Input file containing the model
- 'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 20, 'msTe' : 0.01, 'msLe' : 0.01}, # parameters for input file model
+ 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 3, 'msTe' : 0.1, 'msLe' : 0.1}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -72,7 +74,7 @@ def cfgInviscid(nthrds, verb):
'y_ref' : 0.0, # reference point for moment computation (y)
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
- 'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+ 'LSolver' : 'SparseLu', # inner solver (Pardiso, MUMPS or GMRES)
'G_fill' : 2, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
@@ -87,8 +89,8 @@ def cfgBlast(verb):
'Minf' : 0.2, # Freestream Mach number (used for the computation of the time step only)
'CFL0' : 1, # Inital CFL number of the calculation
'Verb': verb, # Verbosity level of the solver
- 'couplIter': 200, # Maximum number of iterations
- 'couplTol' : 1e-8, # Tolerance of the VII methodology
+ 'couplIter': 100, # Maximum number of iterations
+ 'couplTol' : 1e-10, # Tolerance of the VII methodology
'iterPrint': 20, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0], # List of sections for boundary layer calculation
@@ -109,20 +111,10 @@ def main():
dalpha = 1e-6 * np.pi/180
daoa = [alpha - dalpha, alpha + dalpha]
dcl = []
+ dcd = []
coupler, isol, vsol = viscUtils.initBlast(icfg, vcfg)
-
- print(ccolors.ANSI_BLUE + 'PyInviscid...' + ccolors.ANSI_RESET)
- isol.run()
- # Inviscid finite difference
- for aa in daoa:
- isol.setAoA(aa)
- isol.reset()
- isol.run()
- dcl.append(isol.getCl())
- dcl_daoa_dart_fd = (dcl[-1] - dcl[-2])/(2*dalpha)
- isol.adjointSolver.run()
- dcl_daoa_dart_ad = isol.adjointSolver.dClAoa
+ morpher = isol.iobj['mrf']
# Adjoint objects
cAdj = blast.CoupledAdjoint(isol.adjointSolver, vsol)
@@ -132,7 +124,10 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['forward'].start()
aeroCoeffs = coupler.run()
- tms['forward'].stop()
+ tms['forward'].stop()
+
+ # Print with precision of 1e-12
+ print(f"{isol.solver.Cl:.12f}", f"{isol.solver.Cd:.12f}")
print(ccolors.ANSI_BLUE + 'PyDartAdjoint...' + ccolors.ANSI_RESET)
tms['adj_dart'].start()
@@ -147,28 +142,148 @@ def main():
cAdj.run()
tms['adj_blast'].stop()
- # Coupled finite difference
- tms['fd_blast'].start()
+ # Get mesh matrix
+ dCl_x = np.zeros((isol.iobj['bnd'].nodes.size(), isol.getnDim()))
+ dCd_x = np.zeros((isol.iobj['bnd'].nodes.size(), isol.getnDim()))
+ for inod, n in enumerate(isol.blw[0].nodes):
+ for i in range(isol.getnDim()):
+ dCl_x[inod, i] = cAdj.tdCl_x[n.row][i]
+ dCd_x[inod, i] = cAdj.tdCd_x[n.row][i]
+
+
+ print(ccolors.ANSI_BLUE + 'PyFiniteDifference...' + ccolors.ANSI_RESET)
+ # Finite difference on AoA
+ tms['fd_blaster_aoa'].start()
dcl = []
+ dcd = []
for aa in daoa:
isol.setAoA(aa)
coupler.reset()
aeroCoeffs = coupler.run(write=False)
dcl.append(isol.getCl())
- dcl_daoa_blast_fd = (dcl[-1] - dcl[-2])/(2*dalpha)
- tms['fd_blast'].stop()
-
- print(ccolors.ANSI_BLUE + 'PyResults' + ccolors.ANSI_RESET)
- print('--------- Dart ---------')
- print('AD: ' + str(dcl_daoa_dart_ad))
- print('FD: ' + str(dcl_daoa_dart_fd))
- print('Error: ' + str(abs(dcl_daoa_dart_ad - dcl_daoa_dart_fd)/dcl_daoa_dart_ad))
- print('--------- Blaster ---------')
- print('AD: ' + str(cAdj.tdCl_AoA))
- print('FD: ' + str(dcl_daoa_blast_fd))
- print('Error: ' + str(abs(cAdj.tdCl_AoA - dcl_daoa_blast_fd)/cAdj.tdCl_AoA))
- tms['total'].stop()
+ dcd.append(isol.getCd())
+ 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
+
+ # Recover grad AoA from grad x
+ drot = np.array([[-np.sin(alpha), np.cos(alpha)], [-np.cos(alpha), -np.sin(alpha)]])
+
+ dCl_AoA_fromX = 0.
+ dCd_AoA_fromX = 0.
+ for inod, n in enumerate(isol.iobj['bnd'].nodes):
+ dx = drot.dot([n.pos[0] - 0.25, n.pos[1]])
+ dCl_AoA_fromX += np.array([cAdj.tdCl_x[n.row][0], cAdj.tdCl_x[n.row][1]]).dot(dx)
+ dCd_AoA_fromX += np.array([cAdj.tdCd_x[n.row][0], cAdj.tdCd_x[n.row][1]]).dot(dx)
+
+ ### BLASTER FD MESH ###
+ tms['fd_blaster_msh'].start()
+ saveNodes = np.zeros((len(isol.iobj['pbl'].msh.nodes), isol.getnDim()))
+ for n in isol.iobj['pbl'].msh.nodes:
+ for i in range(isol.getnDim()):
+ saveNodes[n.row, i] = n.pos[i]
+ dCl_x_blaster_fd = np.zeros((len(isol.blw[0].nodes), isol.getnDim()))
+ dCd_x_blaster_fd = np.zeros((len(isol.blw[0].nodes), isol.getnDim()))
+ delta = 1e-6
+ cmp = 0
+ 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)
+ for idim in range(isol.getnDim()):
+ cl = [0, 0]
+ cd = [0, 0]
+ for isgn, sign in enumerate([-1, 1]):
+ 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])
+
+ coupler.reset()
+ # prevent print
+ from contextlib import redirect_stdout
+ with redirect_stdout(None):
+ aeroCoeffs = 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()
+
+ # 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)
+ 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('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()
+
+ # 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
diff --git a/blast/utils.py b/blast/utils.py
index 93efffc..4da3ce6 100644
--- a/blast/utils.py
+++ b/blast/utils.py
@@ -123,7 +123,7 @@ def getSolution(sections, write=False, toW=['x', 'xelm', 'theta', 'H', 'deltaSta
varKeys = ['theta', 'H', 'N', 'ue', 'Ct']
attrKeys = ['x', 'y', 'z', 'xoc', 'deltaStar', \
'cd', 'cf', 'Hstar', 'Hstar2', 'Hk', 'ctEq', 'us', 'delta',\
- 'vt', 'rhoe', 'Ret']
+ 'vt', 'vtExt', 'deltaStarExt', 'rhoe', 'Ret', 'regime']
elemsKeys = ['xelm', 'yelm', 'zelm']
for iSec, sec in enumerate(sections):
diff --git a/blast/validation/raeValidation.py b/blast/validation/raeValidation.py
index 2ee1444..ef462b8 100644
--- a/blast/validation/raeValidation.py
+++ b/blast/validation/raeValidation.py
@@ -132,12 +132,12 @@ def main():
# Test solution
print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
tests = CTests()
- tests.add(CTest('Cl', isol.getCl(), 0.759, 5e-2))
- tests.add(CTest('Cd wake', vsol.Cdt, 0.0093, 1e-3, forceabs=True))
- tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0134, 1e-3, forceabs=True))
+ tests.add(CTest('Cl', isol.getCl(), 0.730, 5e-2))
+ tests.add(CTest('Cd wake', vsol.Cdt, 0.0095, 1e-3, forceabs=True))
+ tests.add(CTest('Cd integral', isol.getCd() + vsol.Cdf, 0.0132, 1e-3, forceabs=True))
tests.add(CTest('Cdf', vsol.Cdf, 0.0068, 1e-3, forceabs=True))
if icfg['LSolver'] == 'SparseLu':
- tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 29, 0, forceabs=True))
+ tests.add(CTest('Iterations', len(aeroCoeffs['Cl']), 21, 0, forceabs=True))
tests.run()
expResults = np.loadtxt(os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/references/rae2822_AR138_case6.dat')
--
GitLab