diff --git a/blast/interfaces/DSolversInterface.py b/blast/interfaces/DSolversInterface.py
index fcf09e9e8e54a4c27e3e8b88bd9b69fc32e2dc80..8ce769401e43e382e873dbf70704ce457c363d07 100644
--- a/blast/interfaces/DSolversInterface.py
+++ b/blast/interfaces/DSolversInterface.py
@@ -54,7 +54,7 @@ class SolversInterface:
self.interpolator.inviscidToViscous(self.iBnd, self.vBnd, couplIter)
for iSec in range(self.cfg['nSections']):
## Mesh
- if (self.vBnd[iSec][0].isMeshChanged()):
+ if (self.vBnd[iSec][0].isMeshChanged()) or _adj == 1:
for reg in self.vBnd[iSec]:
if 'vWake' in reg.name:
iReg = 2
diff --git a/blast/interfaces/dart/DartAdjointInterface.py b/blast/interfaces/dart/DartAdjointInterface.py
index c1aa816b7a96f5e838de7996f931217057577ef9..527e4eb48eea0d2512c21cb1e59908f14e11d7c7 100644
--- a/blast/interfaces/dart/DartAdjointInterface.py
+++ b/blast/interfaces/dart/DartAdjointInterface.py
@@ -15,6 +15,7 @@ class DartAdjointInterface():
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))
#self.__runViscous()
@@ -31,12 +32,14 @@ class DartAdjointInterface():
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
for iNo in range(nNs):
@@ -63,9 +66,20 @@ class DartAdjointInterface():
dRblw_v[:,iNo*nDim + iDim] = (blowing-saveBlw)/delta
self.iSolverAPI.solver.U[irow][iDim] = savev[iNo*nDim + iDim]
+ # # Mesh coordinates
+ 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
+ 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]
+
self.coupledAdjointSolver.setdRblw_M(dRblw_M)
self.coupledAdjointSolver.setdRblw_rho(dRblw_rho)
self.coupledAdjointSolver.setdRblw_v(dRblw_v)
+ self.coupledAdjointSolver.setdRblw_x(dRblw_x)
def __runViscous(self):
self.iSolverAPI.imposeInvBC(1, adj=1)
diff --git a/blast/src/DBoundaryLayer.h b/blast/src/DBoundaryLayer.h
index d56795eb194d8bd9a9d5769553359cb190bcde3f..b2a948126b437c1cda0b7d4f37362d01bd3a8865 100644
--- a/blast/src/DBoundaryLayer.h
+++ b/blast/src/DBoundaryLayer.h
@@ -56,6 +56,7 @@ public:
void setMesh(std::vector<double> x,
std::vector<double> y,
std::vector<double> z);
+ void setxtrF(const double _xtrF) {std::cout << "Setting xtr reg " << name << " " << _xtrF << std::endl; xtrF = _xtrF;};
// Getters.
size_t getnVar() const { return nVar; };
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 96f1a49c397873c1520808f93c1109a5776c64c5..64a530d526389611be26828cdc34c9271641d8b1 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -51,8 +51,6 @@ CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::sh
else
alinsol = isol->linsol;
- sizeBlowing = iadjoint->getSizeBlowing();
-
// Initalize all derivatives
this->reset();
}
@@ -94,18 +92,33 @@ void CoupledAdjoint::reset()
dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*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);
+ dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim*nNd);
dRM_M.setIdentity();
dRrho_rho.setIdentity();
dRv_v.setIdentity();
dRblw_blw.setIdentity();
+ // Objective functions
dCl_phi = Eigen::VectorXd::Zero(nNd);
dCd_phi = Eigen::VectorXd::Zero(nNd);
+ // Angle of attack
dCl_AoA = 0.0;
dCd_AoA = 0.0;
tdCl_AoA = 0.0;
tdCd_AoA = 0.0;
+
+ // Mesh
+ dRx_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNd, nDim*nNd);
+ dCl_x = Eigen::VectorXd::Zero(nDim*nNd);
+ dCd_x = Eigen::VectorXd::Zero(nDim*nNd);
+
+ tdCl_x.resize(nNd, Eigen::Vector3d::Zero());
+ tdCd_x.resize(nNd, Eigen::Vector3d::Zero());
}
void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint)
@@ -131,15 +144,17 @@ void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMat
// Create a new matrix from the deque of triplets
matrixAdjoint.setFromTriplets(triplets.begin(), triplets.end());
matrixAdjoint.prune(0.);
- //matrixAdjoint.makeCompressed();
+ matrixAdjoint.makeCompressed();
}
void CoupledAdjoint::run()
{
gradientswrtInviscidFlow();
gradientswrtBlowingVelocity();
- transposeMatrices();
gradientwrtAoA();
+ gradientwrtMesh();
+
+ transposeMatrices();
// Adjoint matrix
@@ -158,89 +173,59 @@ void CoupledAdjoint::run()
Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols);
buildAdjointMatrix(matrices, A_adjoint);
- /******** CL ********/
- // Build right hand side for Cl
+ //***** 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];
- // std::copy(dCl_phi.data(), dCl_phi.data() + dCl_phi.size(), rhsCl.begin());
Eigen::VectorXd rhsCl_ = Eigen::Map<const Eigen::VectorXd>(rhsCl.data(), rhsCl.size());
- // Solution vector for lambdasCl
- Eigen::VectorXd lambdaCl(A_adjoint.cols());
+ Eigen::VectorXd lambdaCl(A_adjoint.cols()); // Solution vector for lambdasCl
Eigen::Map<Eigen::VectorXd> lambdaCl_(lambdaCl.data(), lambdaCl.size());
-
- // Solve coupled system
- // std::cout << "A_adjoint " << A_adjoint.rows() << " " << A_adjoint.cols() << " " << std::setprecision(12) << A_adjoint.norm() << std::endl;
- // std::cout << "rhsCl " << std::setprecision(12) << Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()).norm() << std::endl;
- // std::cout << "lambdasCl_ " << std::setprecision(12) << lambdasCl_.norm() << std::endl;
-
- // std::cout << "Supposed size " << isol->pbl->msh->nodes.size() + isol->pbl->bBCs[0]->nodes.size() + 2*isol->pbl->bBCs[0]->nodes.size() + isol->pbl->bBCs[0]->nodes.size() + isol->pbl->bBCs[0]->uE.size() << std::endl;
- std::cout << "A_adjoint " << std::setprecision(20) << A_adjoint.norm() << std::endl;
- // // std::cout << "A_adjoint " << std::setprecision(3) << A_adjoint << std::endl;
- std::cout << "rhsCl " << std::setprecision(20) << Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()).norm() << std::endl;
- // for (auto i=0; i<rhsCl.size(); i++)
- // {
- // std::cout << i << " " << rhsCl[i] << "\n";
- // }
-
- // solverTest.compute(A_adjoint);
- // lambdaCl_ = solverTest.solve(Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()));
-
- //std::cout << *alinsol << std::endl;
alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl_.data(), rhsCl_.size()), lambdaCl_);
- // Eigen::MatrixXd denseMat = Eigen::MatrixXd(A_adjoint);
- // std::cout << "det " << denseMat.determinant() << std::endl;
- // std::cout << "norm inv " << denseMat.inverse().norm() << std::endl;
- // Eigen::VectorXd denseVec = Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size());
- // lambdasCl_ = denseMat.partialPivLu().solve(denseVec);
-
- // std::cout << "nIter " << alinsol->getIterations() << std::endl;
- // std::cout << "error " << alinsol->getError() << std::endl;
-
- std::cout << "lambdasCl_ " << std::setprecision(20) << lambdaCl.norm() << std::endl;
-
- // Eigen::VectorXd lambdasClTest = Eigen::VectorXd::Zero(dRphi_phi.cols());
- // Eigen::Map<Eigen::VectorXd> lambdasClTest_(lambdasClTest.data(), lambdasClTest.size());
- // alinsol->compute(dRphi_phi, Eigen::Map<Eigen::VectorXd>(dCl_phi.data(), dCl_phi.size()), lambdasClTest_);
-
- // std::cout << "complete " << std::setprecision(12) << "lambdaClTest " << lambdasClTest.norm() << std::endl;
-
- // tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdasClTest;
-
- // std::cout << "After solver. NIter " << alinsol->getIterations() << ". Error " << alinsol->getError() << std::endl;
- // writeMatrixToFile(lambdasCl_, "lambdasCl4.txt");
- // std::cout << "lambdasCl " << std::setprecision(12) << lambdasCl.norm() << std::endl;
-
- Eigen::VectorXd lambdaClPhi = lambdaCl.block(0, 0, isol->pbl->msh->nodes.size(), 1);
- // writeMatrixToFile(lambdaClPhi, "lambdasClPhi4.txt");
-
- // std::cout << "lambdasClphi " << std::setprecision(12) << lambdaClPhi.norm() << std::endl;
-
- // // Compute total derivative dCl_dAoA
- // std::cout << "dCl_AoA (partial) " << std::setprecision(12) << dCl_AoA << std::endl;
- // std::cout << "dRphi_AoA " << std::setprecision(12) << dRphi_AoA.transpose().norm() << std::endl;
- // std::cout << "pdCl_AoA " << std::setprecision(12) << dCl_AoA << std::endl;
- // std::cout << "dRphi_AoA " << std::setprecision(12) << dRphi_AoA.transpose().norm() << std::endl;
- // std::cout << "lambdaClPhi " << std::setprecision(12) << lambdaClPhi.norm() << std::endl;
- tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaClPhi;
- std::cout << "tdCl_AoA " << std::setprecision(12) << tdCl_AoA << std::endl;
-
- /******** CD ********/
- // Build right hand side for Cd
- // std::vector<double> rhsCd(A_adjoint.cols(), 0.0);
- // std::copy(dCd_phi.data(), dCd_phi.data() + dCd_phi.size(), rhsCd.begin());
-
- // // Solution vector for lambdasCd
- // Eigen::VectorXd lambdasCd(A_adjoint.cols());
- // Eigen::Map<Eigen::VectorXd> lambdasCd_(lambdasCd.data(), lambdasCd.size());
-
- // // Solve coupled system
- // alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCd.data(), rhsCd.size()), lambdasCd_);
- // Eigen::VectorXd lambdaCdPhi = lambdasCd.block(0, 0, isol->pbl->msh->nodes.size(), 1);
-
- // // Compute total derivative dCl_dAoA
- // tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCdPhi;
+ 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());
+
+ 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);
+
+ tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCdPhi; // 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"
+ 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;
+ 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];
+ tdCd_x[n->row](m) = lambdaCl_x[isol->pbl->nDim * (n->row) + m];
+ }
+ }
}
void CoupledAdjoint::gradientswrtInviscidFlow()
@@ -346,6 +331,16 @@ void CoupledAdjoint::gradientwrtAoA(){
// All others are zero.
}
+void CoupledAdjoint::gradientwrtMesh(){
+ 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 = 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();
+}
+
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
@@ -385,38 +380,58 @@ void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_drho){
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
+ 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]));
+ }
+ }
+ // Set matrix
+ dRblw_x.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_x.prune(0.);
+}
void CoupledAdjoint::transposeMatrices(){
- // Transpose all matrices
+ // Transpose all matrices from the LHS of the adjoint system. Not the others.
dRphi_phi = dRphi_phi.transpose();
dRphi_M = dRphi_M.transpose();
dRphi_rho = dRphi_rho.transpose();
dRphi_v = dRphi_v.transpose();
dRphi_blw = dRphi_blw.transpose();
+ //dRphi_x = dRphi_x.transpose();
dRM_phi = dRM_phi.transpose();
dRM_M = dRM_M.transpose();
dRM_rho = dRM_rho.transpose();
dRM_v = dRM_v.transpose();
dRM_blw = dRM_blw.transpose();
+ //dRrho_x = dRrho_x.transpose();
dRrho_phi = dRrho_phi.transpose();
dRrho_M = dRrho_M.transpose();
dRrho_rho = dRrho_rho.transpose();
dRrho_v = dRrho_v.transpose();
dRrho_blw = dRrho_blw.transpose();
+ //dRrho_x = dRrho_x.transpose();
dRv_phi = dRv_phi.transpose();
dRv_M = dRv_M.transpose();
dRv_rho = dRv_rho.transpose();
dRv_v = dRv_v.transpose();
dRv_blw = dRv_blw.transpose();
+ //dRv_x = dRv_x.transpose();
dRblw_phi = dRblw_phi.transpose();
dRblw_M = dRblw_M.transpose();
dRblw_rho = dRblw_rho.transpose();
dRblw_v = dRblw_v.transpose();
dRblw_blw = dRblw_blw.transpose();
+ // dRblw_x = dRblw_x.transpose();
+
+ // dRx_x = dRx_x.transpose();
// writeMatrixToFile(dRphi_phi, "dRphi_phi.txt");
// writeMatrixToFile(dRphi_M, "dRphi_M.txt");
@@ -450,35 +465,35 @@ void CoupledAdjoint::transposeMatrices(){
// // (Cout all norms)
- std::cout << "norm(dRphi_phi) " << std::setprecision(12) << dRphi_phi.norm() << std::endl;
- std::cout << "norm(dRphi_M) " << std::setprecision(12) << dRphi_M.norm() << std::endl;
- std::cout << "norm(dRphi_rho) " << std::setprecision(12) << dRphi_rho.norm() << std::endl;
- std::cout << "norm(dRphi_v) " << std::setprecision(12) << dRphi_v.norm() << std::endl;
- std::cout << "norm(dRphi_blw) " << std::setprecision(12) << dRphi_blw.norm() << std::endl;
-
- std::cout << "norm(dRM_phi) " << std::setprecision(12) << dRM_phi.norm() << std::endl;
- std::cout << "norm(dRM_M) " << std::setprecision(12) << dRM_M.norm() << std::endl;
- std::cout << "norm(dRM_rho) " << std::setprecision(12) << dRM_rho.norm() << std::endl;
- std::cout << "norm(dRM_v) " << std::setprecision(12) << dRM_v.norm() << std::endl;
- std::cout << "norm(dRM_blw) " << std::setprecision(12) << dRM_blw.norm() << std::endl;
-
- std::cout << "norm(dRrho_phi) " << std::setprecision(12) << dRrho_phi.norm() << std::endl;
- std::cout << "norm(dRrho_M) " << std::setprecision(12) << dRrho_M.norm() << std::endl;
- std::cout << "norm(dRrho_rho) " << std::setprecision(12) << dRrho_rho.norm() << std::endl;
- std::cout << "norm(dRrho_v) " << std::setprecision(12) << dRrho_v.norm() << std::endl;
- std::cout << "norm(dRrho_blw) " << std::setprecision(12) << dRrho_blw.norm() << std::endl;
-
- std::cout << "norm(dRv_phi) " << std::setprecision(12) << dRv_phi.norm() << std::endl;
- std::cout << "norm(dRv_M) " << std::setprecision(12) << dRv_M.norm() << std::endl;
- std::cout << "norm(dRv_rho) " << std::setprecision(12) << dRv_rho.norm() << std::endl;
- std::cout << "norm(dRv_v) " << std::setprecision(12) << dRv_v.norm() << std::endl;
- std::cout << "norm(dRv_blw) " << std::setprecision(12) << dRv_blw.norm() << std::endl;
-
- std::cout << "norm(dRblw_phi) " << std::setprecision(12) << dRblw_phi.norm() << std::endl;
- std::cout << "norm(dRblw_M) " << std::setprecision(20) << dRblw_M.norm() << std::endl;
- std::cout << "norm(dRblw_rho) " << std::setprecision(20) << dRblw_rho.norm() << std::endl;
- std::cout << "norm(dRblw_v) " << std::setprecision(20) << dRblw_v.norm() << std::endl;
- std::cout << "norm(dRblw_blw) " << std::setprecision(12) << dRblw_blw.norm() << std::endl;
+ // std::cout << "norm(dRphi_phi) " << std::setprecision(12) << dRphi_phi.norm() << std::endl;
+ // std::cout << "norm(dRphi_M) " << std::setprecision(12) << dRphi_M.norm() << std::endl;
+ // std::cout << "norm(dRphi_rho) " << std::setprecision(12) << dRphi_rho.norm() << std::endl;
+ // std::cout << "norm(dRphi_v) " << std::setprecision(12) << dRphi_v.norm() << std::endl;
+ // std::cout << "norm(dRphi_blw) " << std::setprecision(12) << dRphi_blw.norm() << std::endl;
+
+ // std::cout << "norm(dRM_phi) " << std::setprecision(12) << dRM_phi.norm() << std::endl;
+ // std::cout << "norm(dRM_M) " << std::setprecision(12) << dRM_M.norm() << std::endl;
+ // std::cout << "norm(dRM_rho) " << std::setprecision(12) << dRM_rho.norm() << std::endl;
+ // std::cout << "norm(dRM_v) " << std::setprecision(12) << dRM_v.norm() << std::endl;
+ // std::cout << "norm(dRM_blw) " << std::setprecision(12) << dRM_blw.norm() << std::endl;
+
+ // std::cout << "norm(dRrho_phi) " << std::setprecision(12) << dRrho_phi.norm() << std::endl;
+ // std::cout << "norm(dRrho_M) " << std::setprecision(12) << dRrho_M.norm() << std::endl;
+ // std::cout << "norm(dRrho_rho) " << std::setprecision(12) << dRrho_rho.norm() << std::endl;
+ // std::cout << "norm(dRrho_v) " << std::setprecision(12) << dRrho_v.norm() << std::endl;
+ // std::cout << "norm(dRrho_blw) " << std::setprecision(12) << dRrho_blw.norm() << std::endl;
+
+ // std::cout << "norm(dRv_phi) " << std::setprecision(12) << dRv_phi.norm() << std::endl;
+ // std::cout << "norm(dRv_M) " << std::setprecision(12) << dRv_M.norm() << std::endl;
+ // std::cout << "norm(dRv_rho) " << std::setprecision(12) << dRv_rho.norm() << std::endl;
+ // std::cout << "norm(dRv_v) " << std::setprecision(12) << dRv_v.norm() << std::endl;
+ // std::cout << "norm(dRv_blw) " << std::setprecision(12) << dRv_blw.norm() << std::endl;
+
+ // std::cout << "norm(dRblw_phi) " << std::setprecision(12) << dRblw_phi.norm() << std::endl;
+ // std::cout << "norm(dRblw_M) " << std::setprecision(20) << dRblw_M.norm() << std::endl;
+ // std::cout << "norm(dRblw_rho) " << std::setprecision(20) << dRblw_rho.norm() << std::endl;
+ // std::cout << "norm(dRblw_v) " << std::setprecision(20) << dRblw_v.norm() << std::endl;
+ // std::cout << "norm(dRblw_blw) " << std::setprecision(12) << dRblw_blw.norm() << std::endl;
// // Print all matrices after transpose
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 5f69e3a5c2ec4f6a30d9642b2aa1a607031fc75b..3ae603da7eb523f8da19ef8243770e714620d478 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -41,47 +41,61 @@ public:
std::shared_ptr<blast::Driver> vsol; ///< Viscous solver
std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero)
- double dCl_AoA; ///< Partial derivative of the lift coefficient with respect to the angle of attack
- double dCd_AoA; ///< Partial erivative of the drag coefficient with respect to the angle of attack
+ double tdCl_AoA; ///< Total derivative of the lift coefficient with respect to the angle of attack
+ double tdCd_AoA; ///< Total erivative of the drag coefficient with respect to the angle of attack
- double tdCl_AoA; ///< Total derivative of the lift coefficient with respect to the angle of attack
- double tdCd_AoA; ///< Total erivative of the drag coefficient with respect to the angle of attack
-private:
- size_t sizeBlowing;
+ std::vector<Eigen::Vector3d> tdCl_x; ///< Total derivative of the lift coefficient with respect to the mesh coordinates
+ std::vector<Eigen::Vector3d> tdCd_x; ///< Total derivative of the drag coefficient with respect to the mesh coordinates
+private:
Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_phi; ///< Inviscid flow Jacobian
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_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
- Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi; ///< Derivative of the Mach number residual with respect to the inviscid flow
- 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_blw; ///< Derivative of the Mach number residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi; ///< Derivative of the Mach number residual with respect to the inviscid flow
+ 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_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
Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_phi; ///< Derivative of the density residual with respect to the inviscid flow
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_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
- Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_phi; ///< Derivative of the velocity residual with respect to the inviscid flow
- 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_blw; ///< Derivative of the velocity residual with respect to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_phi; ///< Derivative of the velocity residual with respect to the inviscid flow
+ 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_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> 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_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
+
+ // Objective functions
+ Eigen::VectorXd dCl_phi; ///< Derivative of the lift coefficient with respect to the inviscid flow
+ Eigen::VectorXd dCd_phi; ///< Derivative of the drag coefficient with respect to the inviscid flow
+
+ // Angle of attack
+ Eigen::VectorXd dRphi_AoA; ///< Derivative of the inviscid flow residual with respect to the angle of attack
+ double dCl_AoA; ///< Partial derivative of the lift coefficient with respect to the angle of attack
+ double dCd_AoA; ///< Partial erivative of the drag coefficient with respect to the angle of attack
- Eigen::VectorXd dRphi_AoA; ///< Derivative of the inviscid flow residual with respect to the angle of attack
- Eigen::VectorXd dCl_phi; ///< Derivative of the lift coefficient with respect to the inviscid flow
- Eigen::VectorXd dCd_phi; ///< Derivative of the drag coefficient with respect to the inviscid flow
+ // Mesh
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_x; ///< Derivative of the mesh residual with respect to the mesh coordinates
+ Eigen::VectorXd dCl_x; ///< Partial derivative of the lift coefficient with respect to the mesh coordinates
+ Eigen::VectorXd dCd_x; ///< Partial derivative of the drag coefficient with respect to the mesh coordinates
public:
CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
@@ -93,6 +107,7 @@ public:
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);
private:
void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
@@ -104,6 +119,7 @@ private:
void gradientswrtBlowingVelocity();
void gradientwrtAoA();
+ void gradientwrtMesh();
void transposeMatrices();
void writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename);
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index 0ccb893868c27b0f83dea26cacd55db066a6c03c..9ae18616fde6d7071861850b7cf83caea47b385b 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -2,6 +2,7 @@
#include "DBoundaryLayer.h"
#include "DSolver.h"
#include <iomanip>
+#include <memory>
#define ANSI_COLOR_RED "\x1b[1;31m"
#define ANSI_COLOR_GREEN "\x1b[1;32m"
@@ -282,7 +283,16 @@ void Driver::checkWaves(size_t iPoint, BoundaryLayer *bl)
*/
void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
{
- // Averages solution on transition marker.
+ // if (bl->name == "upper")
+ // std::cout << "xtrFt = " << bl->xtrF << " iPoint " << iPoint << std::endl;
+ // if (bl->name == "lower"){
+ // std::cout << bl->mesh->getnMarkers() << std::endl;
+ // std::cout << "xtrFl = " << bl->xtrF << " iPoint " << iPoint << std::endl;
+ // if (iPoint != 64 && bl->xtrF!=-1){
+ // std::cout << bl->mesh->getnMarkers() << std::endl;
+ // throw;}
+ // }
+ // Averages solution on transition marker
size_t nVar = bl->getnVar();
// Save laminar solution.
@@ -298,8 +308,15 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
bl->regime[k] = 1; // Set Turbulent regime.
// Compute transition location.
- bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1);
- bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);
+ if (bl->xtrF == -1){
+ bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1);
+ bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);
+ }
+ else{
+ bl->xtr = bl->xtrF;
+ bl->xoctr = bl->xtrF; // TODO: Make sure that xoc works.
+ //bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);
+ }
// Percentage of the subsection corresponding to a turbulent flow.
double avTurb = (bl->mesh->getx(iPoint) - bl->xtr) / (bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1));
@@ -323,11 +340,11 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
std::cout << "Warning: Transition point turbulent computation terminated with exit code " << exitCode << std::endl;
// Average both solutions (Now solution @ iPoint is the turbulent solution).
- // bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta.
- // bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H.
- // bl->u[iPoint * nVar + 2] = 9.; // N.
- // bl->u[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue.
- // bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct.
+ bl->u[iPoint * nVar + 0] = avLam * lamSol[0] + avTurb * bl->u[(iPoint)*nVar + 0]; // Theta.
+ bl->u[iPoint * nVar + 1] = avLam * lamSol[1] + avTurb * bl->u[(iPoint)*nVar + 1]; // H.
+ bl->u[iPoint * nVar + 2] = 9.; // N.
+ bl->u[iPoint * nVar + 3] = avLam * lamSol[3] + avTurb * bl->u[(iPoint)*nVar + 3]; // ue.
+ bl->u[iPoint * nVar + 4] = avTurb * bl->u[(iPoint)*nVar + 4]; // Ct.
bl->u[iPoint*nVar+2] = 9.;
// Recompute closures. (The turbulent BC @ iPoint - 1 does not influence laminar closure relations).
@@ -501,7 +518,7 @@ int Driver::outputStatus(double maxMach) const
std::cout << std::setw(10) << "xTr Upper"
<< std::setw(12) << "xTr Lower"
<< std::setw(12) << "Cd" << std::endl;
- std::cout << std::fixed << std::setprecision(2);
+ std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(10) << meanTransitions[0]
<< std::setw(12) << meanTransitions[1];
std::cout << std::fixed << std::setprecision(5);
diff --git a/blast/src/DDriver.h b/blast/src/DDriver.h
index 2937d3dc529e53ed3da7127fc151dd100732626b..3efe285fae39e52ff6502c8f5d9de88cdd068960 100644
--- a/blast/src/DDriver.h
+++ b/blast/src/DDriver.h
@@ -44,12 +44,13 @@ public:
// Getters.
double getxtr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xtr; };
- double getxoctr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xoctr; };
+ double getxoctr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xoctr; };
double getAvrgxoctr(size_t const iRegion) const;
double getRe() const { return Re; };
std::vector<std::vector<double>> getSolution(size_t iSec);
// Setters.
+ void setxtrF(size_t const iSec,size_t const iReg, double const xtrF) const {sections[iSec][iReg]->setxtrF(xtrF);};
void setGlobMesh(size_t iSec, size_t iRegion, std::vector<double> _x, std::vector<double> _y, std::vector<double> _z);
void setVelocities(size_t iSec, size_t iRegion, std::vector<double> _vx, std::vector<double> _vy, std::vector<double> _vz);
void setMe(size_t iSec, size_t iRegion, std::vector<double> _Me);
diff --git a/blast/tests/dart/adjointVII.py b/blast/tests/dart/adjointVII.py
index 5d10de6c1ea22dae8ebbd801490adf36edb1c751..065e06a797465019027b8a4d9614f60fb1a00a7f 100644
--- a/blast/tests/dart/adjointVII.py
+++ b/blast/tests/dart/adjointVII.py
@@ -79,7 +79,7 @@ def cfgInviscid(nthrds, verb):
'z_ref' : 0.0, # reference point for moment computation (z)
# Numerical
'LSolver' : 'PARDISO', # inner solver (Pardiso, MUMPS or GMRES)
- 'G_fill' : 1, # fill-in factor for GMRES preconditioner
+ 'G_fill' : 2, # fill-in factor for GMRES preconditioner
'G_tol' : 1e-5, # tolerance for GMRES
'G_restart' : 50, # restart for GMRES
'Rel_tol' : 1e-6, # relative tolerance on solver residual
@@ -98,7 +98,7 @@ def cfgBlast(verb):
'iterPrint': 5, # int, number of iterations between outputs
'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
'sections' : [0],
- 'xtrF' : [0.2, 0.2],# Forced transition location
+ 'xtrF' : [0.2, 0.2], # Forced transition location
'interpolator' : 'Matching',
'nDim' : 2
}
@@ -113,8 +113,6 @@ import numpy as np
def main():
- alpha = 2*np.pi/180
-
# Timer.
tms = fwk.Timers()
tms['total'].start()
@@ -123,6 +121,8 @@ def main():
icfg = cfgInviscid(args.k, args.verb)
vcfg = cfgBlast(args.verb)
+ alpha = icfg['AoA']*np.pi/180
+
tms['pre'].start()
coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
tms['pre'].stop()
@@ -138,11 +138,12 @@ def main():
tms['TotalAdjoint'].start()
tms['dartAdjoint'].start()
+ """for iReg in range(2):
+ vSolver.setxtrF(0, iReg, vSolver.getxtr(0, iReg))"""
iSolverAPI.adjointSolver.run()
tms['dartAdjoint'].stop()
dCl_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dClAoa
dCd_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dCdAoa
- print('Compute derivative')
tms['FDDerivatives'].start()
pyAdjoint.computeDerivatives()
tms['FDDerivatives'].stop()
@@ -335,12 +336,21 @@ def main():
tms['coupledFD'].start()
+ vcfg['xtrF'] = [vSolver.getxtr(0, i) for i in range(2)]
+
for aoa in aoaV:
coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
iSolverAPI.argOut['pbl'].update(aoa)
- coupler.run()
+ aeroCoeff = coupler.run()
clV.append(iSolverAPI.getCl())
cdV.append(iSolverAPI.getCd())
+ from matplotlib import pyplot as plt
+ #plt.plot(aeroCoeff[:,0], label='cl')
+ plt.plot(aeroCoeff[:,1], 'o', label='cdw')
+ plt.plot(aeroCoeff[:,2], 'x-', label='cdf+cdp')
+ plt.plot(aeroCoeff[:,3], '^-', label='cdf')
+ plt.legend()
+ plt.show()
dCl_dalpha_FD_coupled = (clV[-1] - clV[0]) / (2*dalpha)
dCd_dalpha_FD_coupled = (cdV[-1] - cdV[0]) / (2*dalpha)
tms['coupledFD'].stop()
diff --git a/blast/tests/dart/adjointVIImesh.py b/blast/tests/dart/adjointVIImesh.py
new file mode 100644
index 0000000000000000000000000000000000000000..ec40bb93f7a5d6f22325ba73e5d6196b75e9dba3
--- /dev/null
+++ b/blast/tests/dart/adjointVIImesh.py
@@ -0,0 +1,382 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2022 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+# @author Paul Dechamps
+# @date 2022
+# Test the blast implementation. The test case is a compressible attached transitional flow.
+# Tested functionalities;
+# - Time integration.
+# - Two time-marching techniques (agressive and safe).
+# - Transition routines.
+# - Forced transition.
+# - Compressible flow routines.
+# - Laminar and turbulent flow.
+#
+# Untested functionalities.
+# - Separation routines.
+# - Multiple failure case at one iteration.
+# - Response to unconverged inviscid solver.
+
+# Imports.
+
+import blast.utils as viscUtils
+import blast
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import math
+
+def cfgInviscid(nthrds, verb):
+ import os.path
+ # Parameters
+ return {
+ # Options
+ 'scenario' : 'aerodynamic',
+ 'task' : 'optimization',
+ 'Threads' : nthrds, # number of threads
+ '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
+ 'Dim' : 2, # problem dimension
+ 'Format' : 'gmsh', # save format (vtk or gmsh)
+ # Markers
+ 'Fluid' : 'field', # name of physical group containing the fluid
+ 'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+ 'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+ '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,
+ 'Upstream' : 'upstream',
+ # Freestream
+ 'M_inf' : 0.2, # freestream Mach number
+ 'AoA' : 2., # freestream angle of attack
+ # Geometry
+ 'S_ref' : 1., # reference surface length
+ 'c_ref' : 1., # reference chord length
+ 'x_ref' : .25, # reference point for moment computation (x)
+ 'y_ref' : 0.0, # reference point for moment computation (y)
+ 'z_ref' : 0.0, # reference point for moment computation (z)
+ # Numerical
+ 'LSolver' : 'PARDISO', # 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
+ 'Rel_tol' : 1e-6, # relative tolerance on solver residual
+ 'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+ 'Max_it' : 20, # solver maximum number of iterations
+ }
+
+def cfgBlast(verb):
+ return {
+ 'Re' : 1e7, # Freestream Reynolds number
+ '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
+ 'iterPrint': 5, # int, number of iterations between outputs
+ 'resetInv' : True, # bool, flag to reset the inviscid calculation at every iteration.
+ 'sections' : [0],
+ 'xtrF' : [0.2, 0.2], # Forced transition location
+ 'interpolator' : 'Matching',
+ 'nDim' : 2
+ }
+
+import dart.utils as floU
+import dart.default as floD
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+import numpy as np
+
+def main():
+
+ # Timer.
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ alpha = icfg['AoA']*np.pi/180
+
+ tms['pre'].start()
+ coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ tms['pre'].stop()
+
+ coupledAdjointSolver = blast.CoupledAdjoint(iSolverAPI.adjointSolver, vSolver)
+ from blast.interfaces.dart.DartAdjointInterface import DartAdjointInterface as DartAdjointinterface
+ pyAdjoint = DartAdjointinterface(coupledAdjointSolver, iSolverAPI, vSolver)
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ aeroCoeffs = coupler.run()
+ tms['solver'].stop()
+
+ tms['TotalAdjoint'].start()
+ tms['dartAdjoint'].start()
+ """for iReg in range(2):
+ vSolver.setxtrF(0, iReg, vSolver.getxtr(0, iReg))"""
+ iSolverAPI.adjointSolver.run()
+ tms['dartAdjoint'].stop()
+ dCl_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dClAoa
+ dCd_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dCdAoa
+ tms['FDDerivatives'].start()
+ pyAdjoint.computeDerivatives()
+ tms['FDDerivatives'].stop()
+ tms['CoupledAdjoint'].start()
+ coupledAdjointSolver.run()
+ tms['CoupledAdjoint'].stop()
+ tms['TotalAdjoint'].stop()
+ print(coupledAdjointSolver.tdCl_AoA)
+
+ # Get matrix
+ dCl_x_mat = np.zeros((iSolverAPI.argOut['msh'].nodes.size(), 3))
+ for n in iSolverAPI.argOut['msh'].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]
+
+ quit()
+
+ ################# ONLY DART #################
+
+ 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' : 1., 'msTe' : 0.01, 'msLe' : 0.01}
+ 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()
+
+ # 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()
+
+ # 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)
+
+ quit()
+
+ # ONLY DART
+
+ # define flow variables
+ 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' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
+ 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()
+
+ # set the problem and add fluid, initial/boundary conditions
+ tms['pre'].start()
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil',vsc=True, dbc=True)
+ tms['pre'].stop()
+
+ # 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()
+
+ dalpha = 1e-2
+ aoaV = [alpha-dalpha, alpha+dalpha]
+ clV = []
+ for aoa in aoaV:
+ pbl.update(aoa)
+ solver.reset()
+ solver.run()
+ clV.append(solver.Cl)
+ print(clV)
+ dCl_dalpha_FD = (clV[-1] - clV[0]) / (2*dalpha)
+ print('dCl_dalpha_FD', dCl_dalpha_FD)
+ print('dCl_dalpha_adjointDart', adjoint.dClAoa)
+
+ quit()
+
+ args = parseargs()
+ icfg = cfgInviscid(args.k, args.verb)
+ vcfg = cfgBlast(args.verb)
+
+ tms['pre'].start()
+ coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
+ tms['pre'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+ tms['solver'].start()
+ #aeroCoeffs = coupler.run()
+ iSolverAPI.run()
+ quit()
+ tms['solver'].stop()
+
+ # Display results.
+ print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+ print(' Re M alpha Cl Cd Cdp Cdf Cm')
+ print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(vcfg['Re']/1e6, iSolverAPI.getMinf(), iSolverAPI.getAoA()*180/math.pi, iSolverAPI.getCl(), vSolver.Cdt, vSolver.Cdp, vSolver.Cdf, iSolverAPI.getCm()))
+
+ # Write results to file.
+ vSolution = viscUtils.getSolution(vSolver)
+ vSolution['Cdt_int'] = vSolution['Cdf'] + iSolverAPI.getCd()
+ tms['total'].stop()
+
+ print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+ print('CPU statistics')
+ print(tms)
+ print('SOLVERS statistics')
+ print(coupler.tms)
+
+ # Test solution
+ print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+ tests = CTests()
+ tests.add(CTest('Cl', iSolverAPI.getCl(), 0.230, 5e-2)) # XFOIL 0.2325
+ tests.add(CTest('Cd wake', vSolution['Cdt_w'], 0.0056, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cd integral', vSolution['Cdt_int'], 0.0059, 1e-3, forceabs=True)) # XFOIL 0.00531
+ tests.add(CTest('Cdf', vSolution['Cdf'], 0.00465, 1e-3, forceabs=True)) # XFOIL 0.00084, Cdf = 0.00447
+ tests.add(CTest('xtr_top', vSolution['xtrT'], 0.20, 0.03, forceabs=True)) # XFOIL 0.1877
+ tests.add(CTest('xtr_bot', vSolution['xtrB'], 0.40, 0.01, forceabs=True)) # XFOIL 0.4998
+ tests.add(CTest('Iterations', len(aeroCoeffs), 17, 0, forceabs=True))
+ tests.run()
+
+ # Show results
+ if not args.nogui:
+ iCp = viscUtils.read('Cp_inviscid.dat')
+ vCp = viscUtils.read('Cp_viscous.dat')
+ xfoilCp = viscUtils.read('../../blast/models/references/cpXfoil_n0012.dat', delim=None, skip = 1)
+ xfoilCpInv = viscUtils.read('../../blast/models/references/cpXfoilInv_n0012.dat', delim=None, skip = 1)
+ plotcp = {'curves': [np.column_stack((vCp[:,0], vCp[:,3])),
+ np.column_stack((xfoilCp[:,0], xfoilCp[:,1])),
+ np.column_stack((iCp[:,0], iCp[:,3])),
+ np.column_stack((xfoilCpInv[:,0], xfoilCpInv[:,1]))],
+ 'labels': ['Blast (VII)', 'XFOIL VII', 'DART (inviscid)', 'XFOIL (inviscid)'],
+ 'lw': [3, 3, 2, 2],
+ 'color': ['firebrick', 'darkblue', 'firebrick', 'darkblue'],
+ 'ls': ['-', '-', '--', '--'],
+ 'reverse': True,
+ 'xlim':[0, 1],
+ 'yreverse': True,
+ 'legend': True,
+ 'xlabel': '$x/c$',
+ 'ylabel': '$c_p$'
+ }
+ viscUtils.plot(plotcp)
+
+ xfoilBl = viscUtils.read('../../blast/models/references/blXfoil_n0012.dat', delim=None, skip = 1)
+ plotcf = {'curves': [np.column_stack((vSolution['x'], vSolution['Cf'])),
+ np.column_stack((xfoilBl[:,1], xfoilBl[:,6]))],
+ 'labels': ['Blast', 'XFOIL'],
+ 'lw': [3, 3],
+ 'color': ['firebrick', 'darkblue'],
+ 'ls': ['-', '-'],
+ 'reverse': True,
+ 'xlim':[0, 1],
+ 'legend': True,
+ 'xlabel': '$x/c$',
+ 'ylabel': '$c_f$'
+ }
+ viscUtils.plot(plotcf)
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/blast/tests/dart/blwVelAF.dat b/blast/tests/dart/blwVelAF.dat
new file mode 100644
index 0000000000000000000000000000000000000000..facb64847e6af8b64545210a691620c731950600
--- /dev/null
+++ b/blast/tests/dart/blwVelAF.dat
@@ -0,0 +1,204 @@
+0 0.0015472017424161876
+1 0.0016773600429199993
+2 0.0018863569150160103
+3 0.001763439686266834
+4 0.001568199033448518
+5 0.0017402171824456903
+6 0.0014963297550599833
+7 0.0015198940629153157
+8 0.00126595490150758
+9 0.0011088256951547733
+10 0.001199767413947563
+11 0.001621608990058811
+12 0.0011689248903552197
+13 0.0002906323277981281
+14 0.0010067181753626955
+15 0.001222399531472063
+16 0.0010051533604559643
+17 0.0011054346379531501
+18 0.0009602169553236365
+19 0.0007691866722240861
+20 0.0009154655803897023
+21 0.0007143373179063372
+22 0.0011206069742154586
+23 -0.005880582515548099
+24 -0.0011891420749166079
+25 0.001083967866621263
+26 0.0021072781026718968
+27 0.002440540379804815
+28 0.002665570053442378
+29 0.0027130680212543823
+30 0.002782135614650298
+31 0.002761228852589572
+32 0.002835275907979123
+33 0.002730070774617385
+34 0.0027546041065057768
+35 0.002670695375120869
+36 0.002706689441739721
+37 0.0026474326498508935
+38 0.0027307747087445154
+39 0.0025710532146726753
+40 0.0025928574563471966
+41 0.002632481717826049
+42 0.0025873097987907418
+43 0.002591021968753447
+44 0.002513663301881237
+45 0.002634427401364173
+46 0.002522454918362463
+47 0.0026016772142720508
+48 0.0024733203316227885
+49 0.0025777882116510224
+50 0.0024707858062207357
+51 0.0026049426038305843
+52 0.0024385454193030184
+53 0.0024713191149579807
+54 0.002533716826142388
+55 0.0024874096608077842
+56 0.002551693243932139
+57 0.0024221379226876017
+58 0.0025982053387708674
+59 0.002411197163924921
+60 0.002594073120420692
+61 0.0025671102441198277
+62 0.00238186515236727
+63 0.0025724099185452996
+64 0.002423632102964626
+65 0.002664969073070211
+66 0.002357062567689345
+67 0.0026035573568808463
+68 0.002446702449225442
+69 0.0026208415301456943
+70 0.0024802283040171167
+71 0.0026904401881843124
+72 0.002496806005225793
+73 0.0027198909706144334
+74 0.0025210434743429856
+75 0.0027302434451465076
+76 0.0025591948405131928
+77 0.002863280152124771
+78 0.002665827075573164
+79 0.002903936731973565
+80 0.0027196107159141603
+81 0.0030067740112287933
+82 0.0028367298748751264
+83 0.003283778658718294
+84 0.0029890261051858957
+85 0.003095545741001796
+86 0.003381777291854361
+87 0.003418898062943104
+88 0.0034491335545482685
+89 0.003921041830961181
+90 0.003656615166711191
+91 0.004588382907007209
+92 0.004368001962680508
+93 0.004995362722026925
+94 0.005520030439969944
+95 0.005674627143943854
+96 0.00766628827926289
+97 0.007428439345187231
+98 0.011460250475623998
+99 0.01150693830574791
+100 0.021640744962243292
+101 0.01767973003707986
+102 0.015288996803945342
+103 0.018036265639936305
+104 0.00952936974546074
+105 0.009325491704701932
+106 0.0060851553131460925
+107 0.0062012430631393535
+108 0.004660576230784334
+109 0.004499987219842787
+110 0.004088762930665726
+111 0.0036191610824508704
+112 0.003764704186012373
+113 0.0030922033864582473
+114 0.0032542272029183147
+115 0.0029301006345121617
+116 0.0028818025326883176
+117 0.002872412599605049
+118 0.0026695578896779757
+119 0.0025762637420334925
+120 0.0027755937775263625
+121 0.0024917686373251815
+122 0.002588026732374773
+123 0.00241042943124517
+124 0.0025143677223275116
+125 0.0023783609210405034
+126 0.0024882005429060785
+127 0.002312829117084624
+128 0.0023953375008876924
+129 0.0022887744276228446
+130 0.002401338451244961
+131 0.0022843539607342476
+132 0.0023879330756457727
+133 0.002276310412001205
+134 0.0023413660657052183
+135 0.00226392525888348
+136 0.0023425875999032687
+137 0.0022260389981621606
+138 0.0023833542310969062
+139 0.0022769163395256462
+140 0.002339937766692311
+141 0.002268997755567468
+142 0.002347430319576633
+143 0.0024011401462303854
+144 0.0023118652435184447
+145 0.00240402293259424
+146 0.002339242285062878
+147 0.002390602499244685
+148 0.0023978736417541893
+149 0.0024207040925448894
+150 0.0023988074744847937
+151 0.0023861009441744235
+152 0.002446923332771523
+153 0.002383049215921073
+154 0.002348827676217905
+155 0.0022268703430344907
+156 0.002053806911459491
+157 0.0016666719963497245
+158 0.0009628716319828143
+159 -0.0005131147134452871
+160 -0.0033366046667752365
+161 -0.004966149804816334
+162 0.0008695919401765184
+163 0.00042551218042470677
+164 0.0005106963026773093
+165 0.0008197741640752258
+166 0.0006860586667983148
+167 0.0006735422110081416
+168 0.0006511586878900912
+169 0.0007538061520637727
+170 0.0007050565029041196
+171 0.0008367372560230479
+172 0.0007144597137843431
+173 0.0006687364284581794
+174 0.0006541066894374085
+175 0.0006847665758617066
+176 0.0006662131546093521
+177 0.000886984664124033
+178 0.000769413194522826
+179 0.0006625485852260795
+180 0.0007456469925613395
+181 0.000776275989234347
+182 0.000771150090255642
+183 0.000761864128033165
+184 0.0007594178293172634
+185 0.0007475080496973561
+186 0.0009374102544388237
+187 0.0008119031666674904
+188 0.0009712748062082091
+189 0.0009779312675692188
+190 0.0005625597878524078
+191 0.0007835974978463935
+192 0.0011577817508170807
+193 0.0010147076280330508
+194 0.0009555579196725242
+195 0.001029407335740798
+196 0.0012341821753167799
+197 0.0011783608496652087
+198 0.001219566363561317
+199 0.0013118334705232056
+200 0.0012328892530859045
+201 0.0012856275826508294
+202 0.0011786735893489458
+203 0.0015763669806176381