diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index 8ce769401e43e382e873dbf70704ce457c363d07..bf1cd1a4d5ede687cbb146be8df51fdd3c967911 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -72,7 +72,7 @@ class SolversInterface:
reg.getYUpper(),
reg.getZUpper())
self.vSolver.setGlobMesh(iSec, 1,
- abs(reg.getXLower()),
+ reg.getXLower(),
reg.getYLower(),
reg.getZLower())
# External variables
@@ -90,7 +90,30 @@ class SolversInterface:
struture but name was', reg.name)
## Inviscid variables
for reg in self.vBnd[iSec]:
+ # External variables
reg.updatePrev()
+ if 'vAirfoil' in reg.name:
+ self.xxExt[iSec][0] = np.zeros(reg.stagPoint+1)
+ self.deltaStarExt[iSec][0] = np.zeros(reg.stagPoint+1)
+ #self.vtOld[iSec][0] = np.zeros(reg.stagPoint+1)
+ self.xxExt[iSec][1] = np.zeros(reg.nPoints
+ - reg.stagPoint)
+ self.deltaStarExt[iSec][1] = np.zeros(reg.nPoints
+ - reg.stagPoint)
+ #self.vtOld[iSec][1] = np.zeros(reg.nPoints
+ # - reg.stagPoint)
+ elif 'vWake' in reg.name:
+ iReg = 2
+ self.vSolver.setGlobMesh(iSec, iReg,
+ reg.nodesCoord[:,0],
+ reg.nodesCoord[:,1],
+ reg.nodesCoord[:,2])
+ # External variables
+ self.xxExt[iSec][2] = np.zeros(reg.nPoints)
+ self.deltaStarExt[iSec][2] = np.zeros(reg.nPoints)
+ #self.vtOld[iSec][2] = np.zeros(reg.nPoints)
+ else:
+ raise RuntimeError('Wrong region')
if 'vWake' in reg.name:
iReg = 2
self.vSolver.setVelocities(iSec, iReg,
diff --git a/blast/interfaces/dart/DartAdjointInterface.py b/blast/interfaces/dart/DartAdjointInterface.py
index e9e62d749b6ed0ab810215c8a7be822421c15f00..2bea93d2c80e7ebdf499029f7d1d45f473923c30 100644
--- a/blast/interfaces/dart/DartAdjointInterface.py
+++ b/blast/interfaces/dart/DartAdjointInterface.py
@@ -29,55 +29,76 @@ class DartAdjointInterface():
for i in range(nEs):
saveBlw[i] = self.iSolverAPI.blw[0].getU(i)
- saveMach = np.zeros(nNs)
- saverho = np.zeros(nNs)
- savev = np.zeros(nDim*nNs)
- savex = np.zeros((nNs, nDim))
- for iNo in range(nNs):
- irow = self.iSolverAPI.blw[0].nodes[iNo].row
- saveMach[iNo] = self.iSolverAPI.solver.M[irow]
- saverho[iNo] = self.iSolverAPI.solver.rho[irow]
- for jj in range(nDim):
- savev[iNo*nDim + jj] = self.iSolverAPI.solver.U[irow][jj]
- 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
- blowing = self.__runViscous()
- dRblw_M[:,iNo] = (blowing-saveBlw)/delta
- self.iSolverAPI.solver.M[irow] = saveMach[iNo]
+ for i, n in enumerate(self.iSolverAPI.blw[0].nodes):
+ # vidx = np.where(self.iSolverAPI.iBnd[0].nodesCoord[:,3] == n.row)[0][0]
+ # if self.iSolverAPI.iBnd[0].nodesCoord[vidx,0] != n.pos[0] or self.iSolverAPI.iBnd[0].nodesCoord[vidx,1] != n.pos[1]:
+ # print('WRONG NODE MACH')
+ # quit()
+ savemach = self.iSolverAPI.solver.M[n.row]
+ self.iSolverAPI.solver.M[n.row] = savemach + delta
+ blowing1 = self.__runViscous()
+
+ self.iSolverAPI.solver.M[n.row] = savemach - delta
+ blowing2 = self.__runViscous()
+ self.iSolverAPI.solver.M[n.row] = savemach
+ dRblw_M[:,i] = -(blowing1 - blowing2)/(2*delta)
# Rho
print('Density')
- for iNo in range(nNs):
- irow = self.iSolverAPI.blw[0].nodes[iNo].row
- self.iSolverAPI.solver.rho[irow] += delta
- blowing = self.__runViscous()
- dRblw_rho[:,iNo] = (blowing-saveBlw)/delta
- self.iSolverAPI.solver.rho[irow] = saverho[iNo]
+ for i, n in enumerate(self.iSolverAPI.blw[0].nodes):
+ # vidx = np.where(self.iSolverAPI.iBnd[0].nodesCoord[:,3] == n.row)[0][0]
+ # if self.iSolverAPI.iBnd[0].nodesCoord[vidx,0] != n.pos[0] or self.iSolverAPI.iBnd[0].nodesCoord[vidx,1] != n.pos[1]:
+ # print('WRONG NODE DENSITY')
+ # quit()
+ saverho = self.iSolverAPI.solver.rho[n.row]
+ self.iSolverAPI.solver.rho[n.row] = saverho + delta
+ blowing1 = self.__runViscous()
+
+ self.iSolverAPI.solver.rho[n.row] = saverho - delta
+ blowing2 = self.__runViscous()
+ self.iSolverAPI.solver.rho[n.row] = saverho
+ dRblw_rho[:,i] = -(blowing1 - blowing2)/(2*delta)
# # V
print('Velocity')
- for iNo in range(nNs):
- irow = self.iSolverAPI.blw[0].nodes[iNo].row
- for iDim in range(nDim):
- self.iSolverAPI.solver.U[irow][iDim] += delta
- blowing = self.__runViscous()
- dRblw_v[:,iNo*nDim + iDim] = (blowing-saveBlw)/delta
- self.iSolverAPI.solver.U[irow][iDim] = savev[iNo*nDim + iDim]
+ 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 VELOCITY')
+ # quit()
+ savev = self.iSolverAPI.solver.U[n.row][jj]
+ self.iSolverAPI.solver.U[n.row][jj] = savev + delta
+ blowing1 = self.__runViscous()
+
+ self.iSolverAPI.solver.U[n.row][jj] = savev - delta
+ blowing2 = self.__runViscous()
+ self.iSolverAPI.solver.U[n.row][jj] = savev
+ dRblw_v[:,i*nDim + jj] = -(blowing1 - blowing2)/(2*delta)
# # Mesh coordinates
print('Mesh')
- for iNo in range(nNs):
- 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()
- dRblw_x[:, irow*nDim+iDim] = (blowing-saveBlw)/delta
- self.iSolverAPI.iBnd[0].nodesCoord[iNo, iDim] = savex[iNo, iDim]
+ 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)
+
+ # print('dRblw_M',dRblw_M[dRblw_M!=0])
+ # print('dRblw_rho', dRblw_rho[dRblw_rho!=0])
+ # print('dRblw_v', dRblw_v[dRblw_v!=0])
+ # print('dRblw_x', dRblw_x[dRblw_x!=0])
self.coupledAdjointSolver.setdRblw_M(dRblw_M)
self.coupledAdjointSolver.setdRblw_rho(dRblw_rho)
diff --git a/blast/interfaces/dart/DartInterface.py b/blast/interfaces/dart/DartInterface.py
index e2cc4f7a5af8197f31e1ee881fd269647cd35009..68656e12751e8d4c9b43d7fdf40f72e95b8b8559 100644
--- a/blast/interfaces/dart/DartInterface.py
+++ b/blast/interfaces/dart/DartInterface.py
@@ -152,6 +152,9 @@ class DartInterface(SolversInterface):
if self.cfg['nDim'] == 2:
self.iBnd[n].connectBlowingVel()
for i, blw in enumerate(self.iBnd[n].blowingVel):
+ if n == 1 and blw != 0:
+ print('blw not 0 in wake')
+ quit()
self.blw[n].setU(i, blw)
def setMesh(self):
diff --git a/blast/src/DBoundaryLayer.cpp b/blast/src/DBoundaryLayer.cpp
index dae9606aef41cdda02a72105fe6fbe2159c2025b..ed33ed251c9c7cc72a86947e447f04df2e6e96e6 100644
--- a/blast/src/DBoundaryLayer.cpp
+++ b/blast/src/DBoundaryLayer.cpp
@@ -46,21 +46,21 @@ void BoundaryLayer::setMesh(std::vector<double> x,
mesh->setGlob(x, y, z);
// Resize and reset all variables if needed.
- if (mesh->getDiffnElm() != 0)
- {
- u.resize(nVar * nMarkers, 0.);
- Ret.resize(nMarkers, 0.);
- cd.resize(nMarkers, 0.);
- cf.resize(nMarkers, 0.);
- Hstar.resize(nMarkers, 0.);
- Hstar2.resize(nMarkers, 0.);
- Hk.resize(nMarkers, 0.);
- ctEq.resize(nMarkers, 0.);
- us.resize(nMarkers, 0.);
- deltaStar.resize(nMarkers, 0.);
- delta.resize(nMarkers, 0.);
- regime.resize(nMarkers, 0);
- }
+ // if (mesh->getDiffnElm() != 0)
+ // {
+ u.resize(nVar * nMarkers, 0.);
+ Ret.resize(nMarkers, 0.);
+ cd.resize(nMarkers, 0.);
+ cf.resize(nMarkers, 0.);
+ Hstar.resize(nMarkers, 0.);
+ Hstar2.resize(nMarkers, 0.);
+ Hk.resize(nMarkers, 0.);
+ ctEq.resize(nMarkers, 0.);
+ us.resize(nMarkers, 0.);
+ deltaStar.resize(nMarkers, 0.);
+ delta.resize(nMarkers, 0.);
+ regime.resize(nMarkers, 0);
+ // }
}
/**
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 072c5893c5a7b348ff4c550e6607e2ecc4829422..793d5207304e58025e7678dc4b6631c98c5399c7 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -126,6 +126,22 @@ void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMat
// Create a list of triplets for the larger matrix
std::vector<Eigen::Triplet<double>> triplets;
+ // Sanity check
+ for (size_t irow = 0; irow < matrices.size(); ++irow)
+ {
+ if (irow > 0)
+ if (matrices[irow].size() != matrices[irow-1].size()){
+ std::cout << "WRONG ROW SIZE " << std::endl;
+ throw;
+ }
+ for (size_t icol = 0; icol < matrices[irow].size(); ++icol)
+ if (icol > 0)
+ if (matrices[irow][icol]->rows() != matrices[irow][icol-1]->rows()){
+ std::cout << "WRONG COL SIZE " << std::endl;
+ throw;
+ }
+ }
+
// Iterate over the rows and columns of the vector
int rowOffset = 0;
for (const auto& row : matrices) {
@@ -162,8 +178,8 @@ void CoupledAdjoint::run()
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_rho, &dRM_rho, &dRv_rho, &dRrho_rho, &dRblw_rho},
{&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}
};
@@ -187,20 +203,31 @@ void CoupledAdjoint::run()
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());
+ // std::vector<std::vector<double>> dCf_U = iadjoint->computeGradientCoefficientsFlow(); // {dCl, dCd}/dU
+ // std::vector<double> lamClFlodum= iadjoint->solveFlow(dCf_U[0]);
+ // Eigen::VectorXd lamClFlodum2 = Eigen::Map<Eigen::VectorXd> (lamClFlodum.data(), lamClFlodum.size());
+ // Eigen::VectorXd diff = lambdaClPhi - lamClFlodum2;
+ // for (auto i=0; i<diff.size(); ++i)
+ // std::cout << "dart " << lamClFlodum2[i] << ", vii " << lambdaClPhi[i] << ", diff " << diff[i] << std::endl;
+ // throw;
+
+ // std::vector<std::vector<double>> dCf_U = iadjoint->computeGradientCoefficientsFlow(); // {dCl, dCd}/dU
+ // std::vector<double> lamClFlodum= iadjoint->solveFlow(dCf_U[0]); // lambdaClU = dRu/dU^-T * dCl/dU
+ // lambdaClPhi = Eigen::Map<Eigen::VectorXd>(lamClFlodum.data(), lamClFlodum.size());
tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaClPhi; // Total gradient
// 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);
+ // 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
+ // tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCdPhi; // Total gradient
//***** Gradient wrt mesh coordinates *****//
//*****************************************//
@@ -210,7 +237,14 @@ void CoupledAdjoint::run()
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());
+ // std::vector<std::vector<double>> dCf_U = iadjoint->computeGradientCoefficientsFlow(); // {dCl, dCd}/dU
+ // std::vector<double> lambdaClPhidum = iadjoint->solveFlow(dCf_U[0]); // lambdaClU = dRu/dU^-T * dCl/dU
+ // lambdaClPhi.resize(lambdaClPhidum.size());
+ // for (auto ii = 0; ii < lambdaClPhi.size(); ii++)
+ // lambdaClPhi[ii] = lambdaClPhidum[ii];
+
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_);
for (auto n : isol->pbl->msh->nodes)
@@ -222,6 +256,54 @@ void CoupledAdjoint::run()
tdCd_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m];
}
}
+
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_phi = Eigen::Sparsetrix<double, Eigen::RowMajor> (isol->pbl->nDim * isol->pbl->msh->nodes.size(), isol->pbl->msh->nodes.size());
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_M = Eigen::SparseMatrix<double, Eigen::RowMajor> (isol->pbl->nDim * isol->pbl->msh->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_rho = Eigen::SparseMatrix<double, Eigen::RowMajor> (isol->pbl->nDim * isol->pbl->msh->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_v = Eigen::SparseMatrix<double, Eigen::RowMajor> (isol->pbl->nDim * isol->pbl->msh->nodes.size(), isol->pbl->nDim * isol->pbl->bBCs[0]->nodes.size());
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_blw = Eigen::SparseMatrix<double, Eigen::RowMajor> (isol->pbl->nDim * isol->pbl->msh->nodes.size(), isol->pbl->bBCs[0]->uE.size());
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_x = Eigen::SparseMatrix<double, Eigen::RowMajor> (isol->pbl->bBCs[0]->nodes.size(), isol->pbl->nDim * isol->pbl->msh->nodes.size());
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_x = Eigen::SparseMatrix<double, Eigen::RowMajor> (isol->pbl->bBCs[0]->nodes.size(), isol->pbl->nDim * isol->pbl->msh->nodes.size());
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_x = Eigen::SparseMatrix<double, Eigen::RowMajor> (isol->pbl->bBCs[0]->nodes.size(), isol->pbl->nDim * isol->pbl->msh->nodes.size());
+ // dRx_phi = dRx_phi.transpose();
+ // dRx_M = dRx_M.transpose();
+ // dRx_rho = dRx_rho.transpose();
+ // dRx_v = dRx_v.transpose();
+ // dRx_blw = dRx_blw.transpose();
+ // dRx_x = dRx_x.transpose();
+ // dRphi_x = dRphi_x.transpose();
+ // dRM_x = dRM_x.transpose();
+ // dRv_x = dRv_x.transpose();
+ // dRrho_x = dRrho_x.transpose();
+ // dRblw_x = dRblw_x.transpose();
+ // std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> matricesTest = {
+ // {&dRphi_phi, &dRM_phi, &dRv_phi, &dRrho_phi, &dRblw_phi, &dRx_phi},
+ // {&dRphi_M, &dRM_M, &dRv_M, &dRrho_M, &dRblw_M, &dRx_M},
+ // {&dRphi_rho, &dRM_rho, &dRv_rho, &dRrho_rho, &dRblw_rho, &dRx_rho},
+ // {&dRphi_v, &dRM_v, &dRv_v, &dRrho_v, &dRblw_v, &dRx_v},
+ // {&dRphi_blw, &dRM_blw, &dRv_blw, &dRrho_blw, &dRblw_blw, &dRx_blw},
+ // {&dRphi_x, &dRM_x, &dRv_x, &dRrho_x, &dRblw_x, &dRx_x}
+ // };
+
+ // rows = 0; cols = 0;
+ // for (const auto& row : matricesTest) rows += row[0]->rows(); // All matrices in a row have the same number of rows
+ // for (const auto& mat1 : matricesTest[0]) cols += mat1->cols(); // All matrices in a column have the same number of columns
+ // Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint_test(rows, cols);
+ // buildAdjointMatrix(matricesTest, A_adjoint_test);
+ // std::cout << "MATRIX BUILT " << std::endl;
+
+ // std::vector<double> rhsCl_test(A_adjoint_test.cols(), 0.0);
+ // for (auto i = 0; i<dCl_phi.size(); ++i)
+ // rhsCl_test[i] = dCl_phi[i];
+ // for (auto ii = rhsCl_test.size() - isol->pbl->nDim * isol->pbl->msh->nodes.size(); ii < rhsCl_test.size(); ++ii)
+ // rhsCl_test[ii] = dCl_x[ii];
+ // Eigen::VectorXd rhsCl_test_ = Eigen::Map<const Eigen::VectorXd>(rhsCl_test.data(), rhsCl_test.size());
+
+ // Eigen::VectorXd lambdaCl_test(A_adjoint_test.cols()); // Solution vector for lambdasCl
+ // Eigen::Map<Eigen::VectorXd> lambdaCl_test_(lambdaCl_test.data(), lambdaCl_test.size());
+ // alinsol->compute(A_adjoint_test, Eigen::Map<Eigen::VectorXd>(rhsCl_test_.data(), rhsCl_test_.size()), lambdaCl_test_);
+ // Eigen::VectorXd lambdaCl_x = lambdaCl_test.segment(lambdaCl_test.size() - isol->pbl->nDim * isol->pbl->msh->nodes.size(), isol->pbl->nDim * isol->pbl->msh->nodes.size());
+
}
void CoupledAdjoint::gradientswrtInviscidFlow()
@@ -232,8 +314,8 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
//**** dRM_phi, dRrho_phi, dRv_phi ****//
auto pbl = isol->pbl;
// Finite difference
- std::vector<double> Phi_save = isol->phi; // Differential of the independant variable (phi)
- double deltaPhi = 1e-8; // perturbation
+ std::vector<double> Phi_save = isol->phi; // Differential of the independent variable (phi)
+ double deltaPhi = 1e-6; // perturbation
std::vector<Eigen::Triplet<double>> tripletsdM;
std::vector<Eigen::Triplet<double>> tripletsdrho;
@@ -243,7 +325,7 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
Phi_save[n->row] = isol->phi[n->row] + deltaPhi;
// Recompute the quantities on surface nodes.
- size_t jj = 0;
+ int jj = 0;
for (auto srfNode : pbl->bBCs[0]->nodes){
double dM = 0.; double drho = 0.;
@@ -261,10 +343,11 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
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 (size_t 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));
+ 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
@@ -337,7 +420,7 @@ void CoupledAdjoint::gradientwrtMesh(){
}
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());
+ //dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dM.size(), _dRblw_dM[0].size());
// Convert std::vector<double> to Eigen::Triplets
std::vector<Eigen::Triplet<double>> triplets;
for (size_t i = 0; i < _dRblw_dM.size(); i++){
@@ -350,7 +433,7 @@ void CoupledAdjoint::setdRblw_M(std::vector<std::vector<double>> _dRblw_dM){
dRblw_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());
+ //dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dv.size(), _dRblw_dv[0].size());
// Convert std::vector<double> to Eigen::Triplets
std::vector<Eigen::Triplet<double>> triplets;
for (size_t i = 0; i < _dRblw_dv.size(); i++){
@@ -363,7 +446,7 @@ void CoupledAdjoint::setdRblw_v(std::vector<std::vector<double>> _dRblw_dv){
dRblw_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());
+ //dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_drho.size(), _dRblw_drho[0].size());
// Convert std::vector<double> to Eigen::Triplets
std::vector<Eigen::Triplet<double>> triplets;
for (size_t i = 0; i < _dRblw_drho.size(); i++){
@@ -376,7 +459,7 @@ void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_drho){
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());
+ //dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dx.size(), _dRblw_dx[0].size());
// Convert std::vector<double> to Eigen::Triplets
std::vector<Eigen::Triplet<double>> triplets;
for (size_t i = 0; i < _dRblw_dx.size(); i++){
@@ -546,7 +629,13 @@ void CoupledAdjoint::transposeMatrices(){
// std::cout << "dRblw_M " << dRblw_M.rows() << " " << dRblw_M.cols() << std::endl;
// std::cout << "dRblw_rho " << dRblw_rho.rows() << " " << dRblw_rho.cols() << std::endl;
// std::cout << "dRblw_v " << dRblw_v.rows() << " " << dRblw_v.cols() << std::endl;
- // std::cout << "dRblw_blw " << dRblw_blw.rows() << " " << dRblw_blw.cols() << std::endl;
+ // std::cout << "dRphi_x " << dRphi_x.rows() << " " << dRphi_x.cols() << std::endl;
+ // std::cout << "dRM_x " << dRM_x.rows() << " " << dRM_x.cols() << std::endl;
+ // std::cout << "dRv_x " << dRv_x.rows() << " " << dRv_x.cols() << std::endl;
+ // std::cout << "dRrho_x " << dRrho_x.rows() << " " << dRrho_x.cols() << std::endl;
+ // std::cout << "dRblw_x " << dRblw_x.rows() << " " << dRblw_x.cols() << std::endl;
+ // std::cout << "dRphi_x " << dRphi_x.rows() << " " << dRphi_x.cols() << std::endl;
+ // std::cout << "dRx_x " << dRx_x.rows() << " " << dRx_x.cols() << std::endl;
}
void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename) {
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 9ae18616fde6d7071861850b7cf83caea47b385b..a551acb9b12e88ef2d0719a58736e186f5994d32 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -121,7 +121,7 @@ int Driver::run(unsigned int const couplIter)
tms["0-CheckIn"].stop();
// Loop over the regions (upper, lower and wake).
- for (size_t iRegion = 0; iRegion < sections[iSec].size(); ++iRegion)
+ for (size_t iRegion = 0; iRegion < 2; ++iRegion)
{
convergenceStatus[iSec][iRegion].resize(0);
@@ -132,6 +132,11 @@ int Driver::run(unsigned int const couplIter)
if (verbose > 1)
std::cout << sections[iSec][iRegion]->name << " ";
+ std::vector<double> xxExtCurr(sections[iSec][iRegion]->mesh->getnMarkers(), 0.);
+ for (size_t iPoint = 0; iPoint < xxExtCurr.size(); ++iPoint)
+ xxExtCurr[iPoint] = sections[iSec][iRegion]->mesh->getLoc(iPoint);
+ sections[iSec][iRegion]->mesh->setExt(xxExtCurr);
+
// Check if safe mode is required.
if (couplIter == 0 || sections[iSec][iRegion]->mesh->getDiffnElm() != 0 || stagPtMvmt[iSec] == true || solverExitCode != 0)
{
@@ -243,7 +248,7 @@ int Driver::run(unsigned int const couplIter)
if (verbose > 1)
std::cout << std::endl;
}
-
+ sections[0][2]->computeBlwVel();
for (size_t i = 0; i < sections.size(); ++i)
computeSectionalDrag(sections[i], i);
computeTotalDrag();
diff --git a/blast/tests/dart/adjointVII.py b/blast/tests/dart/adjointVII.py
index 065e06a797465019027b8a4d9614f60fb1a00a7f..9f2cbdbab4e4d12f1f984874c1a1dbd9554b186a 100644
--- a/blast/tests/dart/adjointVII.py
+++ b/blast/tests/dart/adjointVII.py
@@ -93,8 +93,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': 75, # Maximum number of iterations
- 'couplTol' : 1e-4, # Tolerance of the VII methodology
+ 'couplIter': 200, # 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.
'sections' : [0],
diff --git a/blast/tests/dart/adjointVIImesh.py b/blast/tests/dart/adjointVIImesh.py
index 27fde997fee5d04d44e9c2c551f5162111f78a0e..5617a2fd02e42349052ec4c48a88a78e3fae953b 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' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.05, 'msLe' : 0.05}, # parameters for input file model
+ 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.1, 'msLe' : 0.1}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -136,6 +136,7 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
aeroCoeffs = coupler.run()
+ iSolverAPI.run()
tms['solver'].stop()
Clref = iSolverAPI.getCl()
@@ -168,9 +169,37 @@ def main():
for i in range(3):
dCl_x_mat[n.row, i] = coupledAdjointSolver.tdCl_x[n.row][i]
+ ### DART ###
+ delta = 1e-6
+ dCl_x_FD_DART = np.zeros((len(iSolverAPI.argOut['pbl'].msh.nodes), 3))
+ for n in iSolverAPI.argOut['bnd'].nodes:
+ for i in range(iSolverAPI.argOut['pbl'].nDim):
+ #print('BEFORE', iSolverAPI.solver.U[n.row][i])
+ pos = n.pos[i] # eviter les erreurs d'arrondi
+ n.pos[i] = pos + delta # deforme la surface
+ #morpher.deform() # deforme le maillage et mets à jour les elements
+ #solver.reset()
+ iSolverAPI.reset()
+ iSolverAPI.run()
+ Cl1 = iSolverAPI.getCl()
+ #print('+', iSolverAPI.solver.U[n.row][i])
+
+ n.pos[i] = pos - delta # deforme la surface
+ #morpher.deform() # deforme le maillage et mets à jour les elements
+ iSolverAPI.reset()
+ iSolverAPI.run()
+ Cl2 = iSolverAPI.getCl()
+ #print('-', iSolverAPI.solver.U[n.row])
+ dCl_x_FD_DART[n.row, i] = (Cl1 - Cl2)/(2*delta)
+ n.pos[i] = pos
+ #quit()
+ # for n in iSolverAPI.argOut['bnd'].nodes:
+ # print((dCl_x_FD_DART[n.row, :] - dCl_x_mat_DART[n.row, :]) / dCl_x_mat_DART[n.row, :])
+ # quit()
+
### Finite difference ####
print('Starting FD')
- delta = 1e-6
+ delta = 1e-8
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]
@@ -186,21 +215,24 @@ def main():
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()
+ # for i in range(len(iSolverAPI.iBnd[0].elemsCoord[:,0])):
+ # iSolverAPI.blw[0].setU(i, 0.)
+ # iSolverAPI.imposeInvBC(0, adj=1)
+ iSolverAPI.imposeInvBC(0, adj=0)
+ for it in range(100):
+ eci = iSolverAPI.run()
iSolverAPI.imposeInvBC(it, adj=0)
- vSolver.run(it)
+ ec = vSolver.run(it)
+ if ec != 0:
+ raise RuntimeError('Not converged in 1')
iSolverAPI.imposeBlwVel(adj=0)
clCurr = iSolverAPI.getCl()
- if abs((clCurr - clPrev)/clCurr) < 1e-12:
- print('1conv', n.row, it, abs((clCurr - clPrev)/clCurr))
+ if abs((clCurr - clPrev)/clCurr) < 1e-13:
+ #print('1conv', n.row, it, abs((clCurr - clPrev)/clCurr))
break
clPrev = clCurr
- if (it > 49):
- print('Not converged')
+ if (it >= 90):
+ print('Not converged', it)
quit()
Cl1 = iSolverAPI.getCl()
@@ -208,42 +240,160 @@ def main():
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):
+ # for i in range(len(iSolverAPI.iBnd[0].elemsCoord[:,0])):
+ # iSolverAPI.blw[0].setU(i, 0.)
+ # iSolverAPI.imposeInvBC(0, adj=1)
+ iSolverAPI.imposeInvBC(0, adj=0)
+ for it in range(100):
iSolverAPI.run()
iSolverAPI.imposeInvBC(it, adj=0)
- vSolver.run(it)
+ ec = vSolver.run(it)
+ if ec != 0:
+ raise RuntimeError('Not converged in 2')
iSolverAPI.imposeBlwVel(adj=0)
clCurr = iSolverAPI.getCl()
- if abs((clCurr - clPrev)/clCurr) < 1e-12:
- print('2conv', n.row, it, abs((clCurr - clPrev)/clCurr))
+ if abs((clCurr - clPrev)/clCurr) < 1e-13:
+ #print('2conv', n.row, it, abs((clCurr - clPrev)/clCurr))
break
clPrev = clCurr
- if (it > 49):
- print('Not converged')
+ if (it >= 90):
+ print('Not converged', it)
quit()
Cl2 = iSolverAPI.getCl()
n.pos[idim] = pos
iSolverAPI.iBnd[0].nodesCoord[idxV,idim] = pos
- print(n.row, 'fd ', (Cl1 - Cl2)/(2*delta))
+ #print(n.row, 'fd ', (Cl1 - Cl2)/(2*delta))
dCl_x_FD_coupled[n.row, idim] = (Cl1 - Cl2)/(2*delta)
+ #absDiff = (dCl_x_FD_coupled[n.row, :] - dCl_x_mat[n.row, :])
+ #print(f"row {n.row:<5} pos: [{n.pos[0]:<5.3}, {n.pos[1]:<5.3}, {n.pos[2]:<5.3}] Cl1 {Cl1:<5.7} Cl2 {Cl2:<5.7} absDiff {absDiff}")
+
+ ### FD AoA ###
+ delta = 0.000001*np.pi/180
+ iSolverAPI.reset()
+ iSolverAPI.argOut['pbl'].update(alpha+delta)
+ iSolverAPI.imposeInvBC(0, adj=0)
+ for it in range(100):
+ eci = iSolverAPI.run()
+ iSolverAPI.imposeInvBC(it, adj=0)
+ ec = vSolver.run(it)
+ if ec != 0:
+ raise RuntimeError('Not converged in 1')
+ iSolverAPI.imposeBlwVel(adj=0)
+ clCurr = iSolverAPI.getCl()
+ if abs((clCurr - clPrev)/clCurr) < 1e-13:
+ print('1conv', n.row, it, abs((clCurr - clPrev)/clCurr))
+ break
+ clPrev = clCurr
+ if (it >= 90):
+ print('Not converged')
+ quit()
+ Cl1 = iSolverAPI.getCl()
+ print('1', it, Cl1)
+ iSolverAPI.reset()
+ iSolverAPI.argOut['pbl'].update(alpha-delta)
+ for it in range(100):
+ eci = iSolverAPI.run()
+ iSolverAPI.imposeInvBC(it, adj=0)
+ ec = vSolver.run(it)
+ if ec != 0:
+ raise RuntimeError('Not converged in 1')
+ iSolverAPI.imposeBlwVel(adj=0)
+ clCurr = iSolverAPI.getCl()
+ if abs((clCurr - clPrev)/clCurr) < 1e-13:
+ print('1conv', n.row, it, abs((clCurr - clPrev)/clCurr))
+ break
+ clPrev = clCurr
+ if (it > 90):
+ print('Not converged')
+ quit()
+ Cl2 = iSolverAPI.getCl()
+ print('2', it, Cl2)
+ dCl_aoa_coupled_FD = (Cl1 - Cl2)/(2*delta)
+
+ ### FD AoA DART ###
+
+ ### DART ###
+ delta = 1e-6 * np.pi/180
+ #solver.reset()
+ iSolverAPI.argOut['pbl'].update(alpha + delta)
+ iSolverAPI.reset()
+ iSolverAPI.run()
+ Cl1 = iSolverAPI.getCl()
+ #print('+', iSolverAPI.solver.U[n.row][i])
+
+ iSolverAPI.argOut['pbl'].update(alpha - delta)
+ #morpher.deform() # deforme le maillage et mets à jour les elements
+ iSolverAPI.reset()
+ iSolverAPI.run()
+ Cl2 = iSolverAPI.getCl()
- print(coupledAdjointSolver.tdCl_AoA)
+ iSolverAPI.argOut['pbl'].update(alpha)
+
+ dCl_AoA_FD_DART = (Cl1 - Cl2) / (2*delta)
+ print(iSolverAPI.adjointSolver.dClAoa, dCl_AoA_FD_DART, (iSolverAPI.adjointSolver.dClAoa - dCl_AoA_FD_DART)/iSolverAPI.adjointSolver.dClAoa)
+
+ # 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_FROMMESH_FD = 0
+ #dCdAoA = 0
+ for n in iSolverAPI.argOut['bnd'].nodes:
+ dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1]]))
+ dClAoA_FROMMESH_FD += np.array([dCl_x_FD_coupled[n.row,0], dCl_x_FD_coupled[n.row,1]]).dot(dx)
+ #dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]).dot(dx)
+
+ # 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_FROMMESH_ADJ = 0
+ #dCdAoA = 0
for n in iSolverAPI.argOut['bnd'].nodes:
- print('adj node', n.row, dCl_x_mat[n.row, :])
+ dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1]]))
+ dClAoA_FROMMESH_ADJ += np.array([dCl_x_mat[n.row,0], dCl_x_mat[n.row,1]]).dot(dx)
+ #dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]).dot(dx)
+ print(dClAoA_FROMMESH_ADJ)
+
+ print('from fd', dClAoA_FROMMESH_FD)
+ print('from adj', dClAoA_FROMMESH_ADJ)
+
+
+
+
+ #quit()
+ # 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('absDiff node', n.row, (dCl_x_FD_coupled[n.row, :] - dCl_x_mat[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)
+ mat = np.zeros((iSolverAPI.argOut['bnd'].nodes.size(), 6))
+ ii = 0
for n in iSolverAPI.argOut['bnd'].nodes:
- print('fd node', n.row, dCl_x_FD_coupled[n.row, :])
+ mat[ii, :] = [n.row, n.pos[0], dCl_x_mat[n.row,0], dCl_x_FD_coupled[n.row, 0], dCl_x_mat_DART[n.row,0], dCl_x_FD_DART[n.row, 0]]
+ ii +=1
+ mat = mat[mat[:,0].argsort()]
+ print('x')
+ for row in mat:
+ print(' '.join('{:>15.10f}'.format(val) for val in row))
+ ii = 0
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)
+ mat[ii, :] = [n.row, n.pos[1], dCl_x_mat[n.row,1], dCl_x_FD_coupled[n.row, 1], dCl_x_mat_DART[n.row,1], dCl_x_FD_DART[n.row, 1]]
+ ii +=1
+ #print('nodes {:8.8f} x {:8.8f} AdjCoupl {:8.8f} FDCoupl {:8.8f} AdjDart {:8.8f} FDDart {:8.8f}'.format(n.row, n.pos[0], dCl_x_mat[n.row,0], dCl_x_FD_coupled[n.row, 0], dCl_x_mat_DART[n.row,0], dCl_x_FD_DART[n.row, 0]))
+ mat = mat[mat[:,0].argsort()]
+ print('y')
+ for row in mat:
+ print(' '.join('{:>15.10f}'.format(val) for val in row))
+ print('dClaoa_AdjCoupl', coupledAdjointSolver.tdCl_AoA, 'dClaoa_FDCoupl', dCl_aoa_coupled_FD, 'dClaoa_AdjDart', iSolverAPI.adjointSolver.dClAoa, 'dClaoa_FDDart', dCl_AoA_FD_DART)
+ print('err viscous', (coupledAdjointSolver.tdCl_AoA - dCl_aoa_coupled_FD)/coupledAdjointSolver.tdCl_AoA)
+ print('err inviscid', (iSolverAPI.adjointSolver.dClAoa - dCl_AoA_FD_DART)/iSolverAPI.adjointSolver.dClAoa)
quit()
-
# np.savetxt('dCl_x_FD_coupled.dat',dCl_x_FD_coupled)
# print('COUPLED', dCl_x_mat - dCl_x_FD_coupled)
@@ -267,8 +417,8 @@ def main():
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- 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'])
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.05, 'msLe' : 0.05}
+ msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
# create mesh deformation handler
morpher = floD.morpher(msh, dim, 'airfoil')
tms['msh'].stop()