From 3619cea4bfe0c2a25c08520423bdd0b1eece4a0b Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Tue, 26 Mar 2024 23:41:37 +0100
Subject: [PATCH] (fix) Update adjoint
Validating mesh w/ finite difference. Removed morpher from computation of fd
---
blast/coupler.py | 1 +
blast/interfaces/dart/DartAdjointInterface.py | 9 +-
blast/src/DCoupledAdjoint.cpp | 7 +-
blast/tests/dart/adjointVIImesh.py | 336 +++++++++++++++---
4 files changed, 299 insertions(+), 54 deletions(-)
diff --git a/blast/coupler.py b/blast/coupler.py
index ba45552..64fb433 100644
--- a/blast/coupler.py
+++ b/blast/coupler.py
@@ -69,6 +69,7 @@ class Coupler:
self.iSolverAPI.imposeBlwVel()
self.tms['processing'].stop()
+ # TODO Change aeroCoeffs to dict
aeroCoeffs = np.vstack((aeroCoeffs, [self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.Cdf + self.iSolverAPI.getCd(), self.vSolver.Cdf]))
# Check convergence status.
diff --git a/blast/interfaces/dart/DartAdjointInterface.py b/blast/interfaces/dart/DartAdjointInterface.py
index 527e4eb..e9e62d7 100644
--- a/blast/interfaces/dart/DartAdjointInterface.py
+++ b/blast/interfaces/dart/DartAdjointInterface.py
@@ -23,7 +23,7 @@ class DartAdjointInterface():
# for iRegv in range(2):
# self.deltaStarExt[0][iRegv] = self.vSolver.extractDeltaStar(0, iRegv)
- delta = 1e-5
+ delta = 1e-6
saveBlw = np.zeros(nEs)
for i in range(nEs):
@@ -42,6 +42,7 @@ class DartAdjointInterface():
savex[iNo, jj] = self.iSolverAPI.iBnd[0].nodesCoord[iNo, jj]
# Mach
+ print('Mach')
for iNo in range(nNs):
irow = self.iSolverAPI.blw[0].nodes[iNo].row
self.iSolverAPI.solver.M[irow] += delta
@@ -50,6 +51,7 @@ class DartAdjointInterface():
self.iSolverAPI.solver.M[irow] = saveMach[iNo]
# Rho
+ print('Density')
for iNo in range(nNs):
irow = self.iSolverAPI.blw[0].nodes[iNo].row
self.iSolverAPI.solver.rho[irow] += delta
@@ -58,6 +60,7 @@ class DartAdjointInterface():
self.iSolverAPI.solver.rho[irow] = saverho[iNo]
# # V
+ print('Velocity')
for iNo in range(nNs):
irow = self.iSolverAPI.blw[0].nodes[iNo].row
for iDim in range(nDim):
@@ -67,9 +70,9 @@ class DartAdjointInterface():
self.iSolverAPI.solver.U[irow][iDim] = savev[iNo*nDim + iDim]
# # Mesh coordinates
+ print('Mesh')
for iNo in range(nNs):
- PAS BON CA. IL FAUT ALLER RECHERCHER LA BONNE ROW DANS nodesCoord (colonne 4 normalement).
- irow = self.iSolverAPI.blw[0].nodes[iNo].row
+ irow = int(self.iSolverAPI.iBnd[0].nodesCoord[iNo,3])
for iDim in range(nDim):
self.iSolverAPI.iBnd[0].nodesCoord[iNo, iDim] += delta
blowing = self.__runViscous()
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 64a530d..072c589 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -213,16 +213,12 @@ void CoupledAdjoint::run()
Eigen::VectorXd rhsCl_x = dCl_x - dRphi_x.transpose() * lambdaClPhi - dRblw_x.transpose() * lambdaClBlw;
alinsol->compute(dRx_x.transpose(), Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size()), lambdaCl_x_);
- std::cout << "dRblw_x.transpose().norm() " << std::setprecision(15) << dRblw_x.transpose().norm() << std::endl;
- std::cout << " lambdaClBlw.norm() " << std::setprecision(15) << lambdaClBlw.norm() << std::endl;
- std::cout << "normClAoA " << std::setprecision(15) << tdCl_AoA << " Norm lambdaCl_x " << lambdaCl_x.norm() << std::endl;
- throw;
-
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: CHANGER POUR CD
tdCd_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m];
}
}
@@ -337,7 +333,6 @@ void CoupledAdjoint::gradientwrtMesh(){
dCd_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
dRphi_x = iadjoint->getRu_X();
- std::cout << "AT input dRphi_x " << dRphi_x.rows() << " " << dRphi_x.cols() << " " << dRphi_x.norm() << std::endl;
dRx_x = iadjoint->getRx_X();
}
diff --git a/blast/tests/dart/adjointVIImesh.py b/blast/tests/dart/adjointVIImesh.py
index ec40bb9..27fde99 100644
--- a/blast/tests/dart/adjointVIImesh.py
+++ b/blast/tests/dart/adjointVIImesh.py
@@ -56,7 +56,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.0075}, # parameters for input file model
+ 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.05, 'msLe' : 0.05}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -66,7 +66,7 @@ def cfgInviscid(nthrds, verb):
'Wake' : 'wake', # LIST of names of physical group containing the wake
'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
'Te' : 'te', # LIST of names of physical group containing the trailing edge
- 'dbc' : True,
+ 'dbc' : False,
'Upstream' : 'upstream',
# Freestream
'M_inf' : 0.2, # freestream Mach number
@@ -93,10 +93,10 @@ 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': 75, # Maximum number of iterations
- 'couplTol' : 1e-4, # Tolerance of the VII methodology
+ 'couplIter': 150, # Maximum number of iterations
+ 'couplTol' : 1e-12, # Tolerance of the VII methodology
'iterPrint': 5, # int, number of iterations between outputs
- 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
+ 'resetInv' : False, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0],
'xtrF' : [0.2, 0.2], # Forced transition location
'interpolator' : 'Matching',
@@ -128,6 +128,8 @@ def main():
tms['pre'].stop()
coupledAdjointSolver = blast.CoupledAdjoint(iSolverAPI.adjointSolver, vSolver)
+ import modules.dartflo.dart.default as floD
+ morpher = floD.morpher(iSolverAPI.argOut['msh'], iSolverAPI.argOut['pbl'].nDim, 'airfoil')
from blast.interfaces.dart.DartAdjointInterface import DartAdjointInterface as DartAdjointinterface
pyAdjoint = DartAdjointinterface(coupledAdjointSolver, iSolverAPI, vSolver)
@@ -136,6 +138,8 @@ def main():
aeroCoeffs = coupler.run()
tms['solver'].stop()
+ Clref = iSolverAPI.getCl()
+
tms['TotalAdjoint'].start()
tms['dartAdjoint'].start()
"""for iReg in range(2):
@@ -151,32 +155,108 @@ def main():
coupledAdjointSolver.run()
tms['CoupledAdjoint'].stop()
tms['TotalAdjoint'].stop()
- print(coupledAdjointSolver.tdCl_AoA)
- # Get matrix
+ # Get inviscid matrix
+ dCl_x_mat_DART = np.zeros((iSolverAPI.argOut['msh'].nodes.size(), 3))
+ for n in iSolverAPI.argOut['bnd'].nodes:
+ for i in range(3):
+ dCl_x_mat_DART[n.row, i] = iSolverAPI.adjointSolver.dClMsh[n.row][i]
+
+ # Get coupled matrix
dCl_x_mat = np.zeros((iSolverAPI.argOut['msh'].nodes.size(), 3))
- for n in iSolverAPI.argOut['msh'].nodes:
+ for n in iSolverAPI.argOut['bnd'].nodes:
for i in range(3):
dCl_x_mat[n.row, i] = coupledAdjointSolver.tdCl_x[n.row][i]
-
- print(dCl_x_mat)
- print(np.shape(dCl_x_mat))
- print(coupledAdjointSolver.tdCl_AoA, np.linalg.norm(dCl_x_mat))
-
- # Get matrix
- dCl_x_mat_DART = np.zeros((iSolverAPI.argOut['msh'].nodes.size(), 3))
- for n in iSolverAPI.argOut['msh'].nodes:
- for i in range(3):
- dCl_x_mat_DART[n.row, i] = iSolverAPI.adjointSolver.dClMsh[n.row][i]
+ ### Finite difference ####
+ print('Starting FD')
+ delta = 1e-6
+ dCl_x_FD_coupled = np.zeros((iSolverAPI.argOut['msh'].nodes.size(),3))
+ for n in iSolverAPI.argOut['bnd'].nodes:
+ idxV = np.where(iSolverAPI.iBnd[0].nodesCoord[:,3] == n.row)[0][0]
+ for idim in range(iSolverAPI.argOut['pbl'].nDim):
+ if iSolverAPI.iBnd[0].nodesCoord[idxV,idim] != n.pos[idim]:
+ print('WRONG NODE ', n.pos[i], n.row)
+ quit()
+ #morpher.savePos()
+ pos = n.pos[idim]
+ n.pos[idim] = pos + delta
+ iSolverAPI.iBnd[0].nodesCoord[idxV,idim] = pos + delta
+ #morpher.deform() # deforme le maillage et mets à jour les elements
+
+ clPrev = 0
+ iSolverAPI.reset()
+ for i in range(len(iSolverAPI.iBnd[0].elemsCoord[:,0])):
+ iSolverAPI.blw[0].setU(i, 0.)
+ iSolverAPI.imposeInvBC(0, adj=1)
+ for it in range(50):
+ iSolverAPI.run()
+ iSolverAPI.imposeInvBC(it, adj=0)
+ vSolver.run(it)
+ iSolverAPI.imposeBlwVel(adj=0)
+ clCurr = iSolverAPI.getCl()
+ if abs((clCurr - clPrev)/clCurr) < 1e-12:
+ print('1conv', n.row, it, abs((clCurr - clPrev)/clCurr))
+ break
+ clPrev = clCurr
+ if (it > 49):
+ print('Not converged')
+ quit()
+ Cl1 = iSolverAPI.getCl()
+
+ n.pos[idim] = pos - delta
+ iSolverAPI.iBnd[0].nodesCoord[idxV,idim] = pos - delta
+ clPrev = 0
+ iSolverAPI.reset()
+ for i in range(len(iSolverAPI.iBnd[0].elemsCoord[:,0])):
+ iSolverAPI.blw[0].setU(i, 0.)
+ iSolverAPI.imposeInvBC(0, adj=1)
+ for it in range(50):
+ iSolverAPI.run()
+ iSolverAPI.imposeInvBC(it, adj=0)
+ vSolver.run(it)
+ iSolverAPI.imposeBlwVel(adj=0)
+ clCurr = iSolverAPI.getCl()
+ if abs((clCurr - clPrev)/clCurr) < 1e-12:
+ print('2conv', n.row, it, abs((clCurr - clPrev)/clCurr))
+ break
+ clPrev = clCurr
+ if (it > 49):
+ print('Not converged')
+ quit()
+ Cl2 = iSolverAPI.getCl()
+
+ n.pos[idim] = pos
+ iSolverAPI.iBnd[0].nodesCoord[idxV,idim] = pos
+ print(n.row, 'fd ', (Cl1 - Cl2)/(2*delta))
+ dCl_x_FD_coupled[n.row, idim] = (Cl1 - Cl2)/(2*delta)
+ print(coupledAdjointSolver.tdCl_AoA)
+ for n in iSolverAPI.argOut['bnd'].nodes:
+ print('adj node', n.row, dCl_x_mat[n.row, :])
+ for n in iSolverAPI.argOut['bnd'].nodes:
+ print('fd node', n.row, dCl_x_FD_coupled[n.row, :])
+ for n in iSolverAPI.argOut['bnd'].nodes:
+ print('relDiff node', n.row, (dCl_x_FD_coupled[n.row, :] - dCl_x_mat[n.row, :]) / dCl_x_mat[n.row, :])
+ print('norm dCl_x_mat', np.linalg.norm(dCl_x_mat))
+ print('norm dCl_x_FD_coupled', np.linalg.norm(dCl_x_FD_coupled))
+ print('norm dCl_x_adj_DART', np.linalg.norm(dCl_x_mat_DART))
+ print('dCl_AoA', coupledAdjointSolver.tdCl_AoA)
quit()
+
+ # np.savetxt('dCl_x_FD_coupled.dat',dCl_x_FD_coupled)
- ################# ONLY DART #################
+ # print('COUPLED', dCl_x_mat - dCl_x_FD_coupled)
+ # print('DART', dCl_x_mat_DART - dCl_x_FD_dart)
+ # print(coupledAdjointSolver.tdCl_AoA)
- args = parseargs()
- icfg = cfgInviscid(args.k, args.verb)
- vcfg = cfgBlast(args.verb)
+ # quit()
+
+ # # # # ################# ONLY DART #################
+
+ # # args = parseargs()
+ # # icfg = cfgInviscid(args.k, args.verb)
+ # # vcfg = cfgBlast(args.verb)
# define flow variables
alpha = 2*np.pi/180
@@ -187,8 +267,8 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 2., 'msTe' : 0.01, 'msLe' : 0.01}
+ msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
# create mesh deformation handler
morpher = floD.morpher(msh, dim, 'airfoil')
tms['msh'].stop()
@@ -198,10 +278,9 @@ def main():
pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', vsc=True)
tms['pre'].stop()
- baseBlw = np.loadtxt('/home/paulzer/lab/blaster/blast/tests/dart/blwVelAF.dat')
- for i in range(len(blw[0].nodes)-1):
- blw[0].setU(int(baseBlw[i,0]), baseBlw[i,1])
-
+ # baseBlw = np.loadtxt('/home/paulzer/lab/blaster/blast/tests/dart/blwVelAF.dat')
+ # for i in range(len(blw[0].nodes)-1):
+ # blw[0].setU(int(baseBlw[i,0]), baseBlw[i,1])
# solve problem
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -219,25 +298,192 @@ def main():
adjoint.save(gmshWriter)
tms['adjoint'].stop()
- # compute norm of gradient wrt to mesh
- dClX = 0
- dCdX = 0
- for n in msh.nodes:
- dClX += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1]]).dot(np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1]]))
- dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]))
- dClX = np.sqrt(dClX)
- dCdX = np.sqrt(dCdX)
-
- # recover gradient wrt to AoA from gradient wrt mesh
- drot = np.array([[-np.sin(alpha), np.cos(alpha)], [-np.cos(alpha), -np.sin(alpha)]])
- dClAoA = 0
- dCdAoA = 0
+ # Get matrix
+ dCl_x_mat_DART = np.zeros((msh.nodes.size(), 3))
+ for n in bnd.nodes:
+ for i in range(3):
+ dCl_x_mat_DART[n.row, i] = adjoint.dClMsh[n.row][i]
+ print(dCl_x_mat_DART)
+
+ dCl_x_FD = np.zeros((len(pbl.msh.nodes), 3))
+ Clref = solver.Cl
+ delta = 1e-8
+
+ nodesPos = np.zeros((msh.nodes.size(), 3))
+
+ morpher.savePos() # garder en memoire la position initiale
+ for n in bnd.nodes:
+ for i in range(pbl.nDim):
+ pos = n.pos[i] # eviter les erreurs d'arrondi
+ nodesPos[n.row, i] = pos
+ n.pos[i] = pos + delta # deforme la surface
+ #morpher.deform() # deforme le maillage et mets à jour les elements
+ #solver.reset()
+ solver.run()
+ Cl1 = solver.Cl
+ n.pos[i] = pos - delta # deforme la surface
+ #morpher.deform() # deforme le maillage et mets à jour les elements
+ solver.run()
+ Cl2 = solver.Cl
+ dCl_x_FD[n.row, i] = (Cl1 - Cl2)/(2*delta)
+ n.pos[i] = pos
for n in bnd.nodes:
- dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1]]))
- dClAoA += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1]]).dot(dx)
- dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]).dot(dx)
+ print((dCl_x_FD[n.row, :] - dCl_x_mat_DART[n.row, :]) / dCl_x_mat_DART[n.row, :])
+ print('norm adj', np.linalg.norm(dCl_x_mat_DART))
+ print('norm FD', np.linalg.norm(dCl_x_FD))
+
+ quit()
+ # print(nodesPos)
+ # idxNodeInterest = np.where(abs(nodesPos[:,0] - 8.19144e-2)<1e-4)[0]
+ # print('nodePos', nodesPos[idxNodeInterest,:])
+ # print('adjoint',dCl_x_mat_DART[idxNodeInterest,:])
+ # print('fd', dCl_x_FD[idxNodeInterest,:])
+ # diffabs = dCl_x_FD - dCl_x_mat_DART
+ # print('diffabs', diffabs[idxNodeInterest,:])
+ # diff = (dCl_x_FD - dCl_x_mat_DART)/dCl_x_mat_DART
+ # print('diffrel', diff[idxNodeInterest,:])
+
+
+ # ################# ONLY DART AOA #################
+
+ # args = parseargs()
+ # icfg = cfgInviscid(args.k, args.verb)
+ # vcfg = cfgBlast(args.verb)
+
+ # # define flow variables
+ # alpha = 2*np.pi/180
+ # M_inf = 0.2
+ # c_ref = 1
+ # dim = 2
+
+ # # mesh the geometry
+ # print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+ # tms['msh'].start()
+ # pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 2., 'msTe' : 0.1, 'msLe' : 0.1}
+ # msh, gmshWriter = floD.mesh(dim, 'models/rae2822.geo', pars, ['field', 'airfoil', 'downstream'])
+ # # create mesh deformation handler
+ # morpher = floD.morpher(msh, dim, 'airfoil')
+ # tms['msh'].stop()
+
+ # # set the problem and add fluid, initial/boundary conditions
+ # tms['pre'].start()
+ # pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', vsc=True)
+ # tms['pre'].stop()
+
+ # # baseBlw = np.loadtxt('/home/paulzer/lab/blaster/blast/tests/dart/blwVelAF.dat')
+ # # for i in range(len(blw[0].nodes)-1):
+ # # blw[0].setU(int(baseBlw[i,0]), baseBlw[i,1])
+
+ # # solve problem
+ # print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ # tms['solver'].start()
+ # solver = floD.newton(pbl)
+ # solver.run()
+ # solver.save(gmshWriter)
+ # tms['solver'].stop()
+
+ # # solve adjoint problem
+ # print(ccolors.ANSI_BLUE + 'PyGradient...' + ccolors.ANSI_RESET)
+ # tms['adjoint'].start()
+ # adjoint = floD.adjoint(solver, morpher)
+ # adjoint.run()
+ # adjoint.save(gmshWriter)
+ # tms['adjoint'].stop()
+
+ # # Get matrix
+ # dClAoaAdj = adjoint.dClAoa
+ # delta = 1e-6
+ # clv = []
+ # pbl.update(alpha + delta)
+ # solver.run()
+ # cl1 = solver.Cl
+
+ # pbl.update(alpha - delta)
+ # solver.run()
+ # cl2 = solver.Cl
+
+ # dClAoaFD = (cl1 - cl2) / (2*delta)
+
+ # print('adj:', dClAoaAdj*np.pi/180)
+ # print('fd:', dClAoaFD*np.pi/180)
+
+ ##### ONLY DART THROUGH API ####
+
+ # args = parseargs()
+ # icfg = cfgInviscid(args.k, args.verb)
+ # vcfg = cfgBlast(args.verb)
+
+ # coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+
+ # # baseBlw = np.loadtxt('/home/paulzer/lab/blaster/blast/tests/dart/blwVelAF.dat')
+ # # for i in range(len(blw[0].nodes)-1):
+ # # blw[0].setU(int(baseBlw[i,0]), baseBlw[i,1])
+
+ # iSolverAPI.solver.run()
+
+ # # solve adjoint problem
+ # print(ccolors.ANSI_BLUE + 'PyGradient...' + ccolors.ANSI_RESET)
+ # tms['adjoint'].start()
+ # iSolverAPI.adjointSolver.run()
+ # tms['adjoint'].stop()
+
+ # # Get matrix
+ # dCl_x_mat_DART = np.zeros((iSolverAPI.argOut['msh'].nodes.size(), 3))
+ # for n in iSolverAPI.argOut['bnd'].nodes:
+ # for i in range(3):
+ # dCl_x_mat_DART[n.row, i] = iSolverAPI.adjointSolver.dClMsh[n.row][i]
+ # print(dCl_x_mat_DART)
+
+ # dCl_x_FD = np.zeros((len(iSolverAPI.argOut['pbl'].msh.nodes), 3))
+ # Clref = iSolverAPI.solver.Cl
+ # delta = 1e-6
+
+ # nodesPos = np.zeros((iSolverAPI.argOut['pbl'].msh.nodes.size(), 3))
+
+ # for n in iSolverAPI.argOut['bnd'].nodes:
+ # for i in range(iSolverAPI.argOut['pbl'].nDim):
+ # iSolverAPI.argOut['mrf'].savePos() # garder en memoire la position initiale
+ # pos = n.pos[i] # eviter les erreurs d'arrondi
+ # nodesPos[n.row, i] = pos
+ # n.pos[i] = pos + delta # deforme la surface
+ # iSolverAPI.argOut['mrf'].deform() # deforme le maillage et mets à jour les elements
+ # iSolverAPI.solver.reset()
+ # iSolverAPI.solver.run()
+ # Cl1 = iSolverAPI.solver.Cl
+ # n.pos[i] = pos - delta # deforme la surface
+ # iSolverAPI.argOut['mrf'].deform() # deforme le maillage et mets à jour les elements
+ # iSolverAPI.solver.run()
+ # Cl2 = iSolverAPI.solver.Cl
+ # dCl_x_FD[n.row, i] = (Cl1 - Cl2)/(2*delta)
+ # n.pos[i] = pos
+ # print('pos', nodesPos)
+ # print('adjoint',dCl_x_mat_DART)
+ # print('fd', dCl_x_FD)
+ # diff = (dCl_x_FD - dCl_x_mat_DART)/dCl_x_mat_DART
+ # import math
+ # diff[np.isnan(diff)] = 0.
+ # print(diff[diff!=0])
+
+
+ # compute norm of gradient wrt to mesh
+ # dClX = 0
+ # dCdX = 0
+ # for n in msh.nodes:
+ # dClX += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1]]).dot(np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1]]))
+ # dCdX += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]).dot(np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]))
+ # dClX = np.sqrt(dClX)
+ # dCdX = np.sqrt(dCdX)
+
+ # # recover gradient wrt to AoA from gradient wrt mesh
+ # drot = np.array([[-np.sin(alpha), np.cos(alpha)], [-np.cos(alpha), -np.sin(alpha)]])
+ # dClAoA = 0
+ # dCdAoA = 0
+ # for n in bnd.nodes:
+ # dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1]]))
+ # dClAoA += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1]]).dot(dx)
+ # dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]).dot(dx)
- print(dClAoA, adjoint.dClAoa)
+ # print(dClAoA, adjoint.dClAoa)
quit()
--
GitLab