From 71e1331b04962e61000a9e4a0a0cfc32a6324e10 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Tue, 27 Feb 2024 18:20:30 +0100
Subject: [PATCH] (feat) Adjoint: Validating inviscid part.
Bug with linear solver?
---
blast/interfaces/dart/DartAdjointInterface.py | 234 +++++++++++++-----
blast/src/DCoupledAdjoint.cpp | 212 ++++++++++++----
blast/src/DCoupledAdjoint.h | 2 +
blast/tests/dart/adjointVII.py | 187 ++++++++++++--
4 files changed, 510 insertions(+), 125 deletions(-)
diff --git a/blast/interfaces/dart/DartAdjointInterface.py b/blast/interfaces/dart/DartAdjointInterface.py
index cb3d632..f21e495 100644
--- a/blast/interfaces/dart/DartAdjointInterface.py
+++ b/blast/interfaces/dart/DartAdjointInterface.py
@@ -7,77 +7,179 @@ class DartAdjointInterface():
self.vSolver = _vSolver
def computeDerivatives(self):
+ nDim = self.iSolverAPI.solver.pbl.nDim
+ nNd = self.iSolverAPI.solver.pbl.msh.nodes.size()
+ nEs = len(self.iSolverAPI.iBnd[0].elemsCoord[:,0])
+ nNs = len(self.iSolverAPI.iBnd[0].nodesCoord[:,0])
+ print(f"{nDim=}")
+ print(f"{nNd=}")
+ print(f"{nEs=}")
+ print(f"{nNs=}")
+ dRblw_M = np.zeros((nEs, nNs))
+ dRblw_rho = np.zeros((nEs, nNs))
+ dRblw_v = np.zeros((nEs, nDim*nNs))
+
+ delta = 1e-4
+
+ saveBlw = np.zeros(nEs)
+ 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)
+ for i in range(nNs):
+ saveMach[i] = self.iSolverAPI.solver.M[i]
+ saverho[i] = self.iSolverAPI.solver.rho[i]
+ for j in range(nDim):
+ savev[i*nDim + j] = self.iSolverAPI.solver.U[i][j]
+
+ # Mach
+ for iNo in range(nNs):
+ self.iSolverAPI.solver.M[i] += delta
+ blowing = self.__runViscous()
+ dRblw_M[:,iNo] = (blowing-saveBlw)/delta
+ self.iSolverAPI.solver.M[i] = saveMach[i]
+
+ # Rho
+ for iNo in range(nNs):
+ self.iSolverAPI.solver.rho[i] += delta
+ blowing = self.__runViscous()
+ dRblw_rho[:,iNo] = (blowing-saveBlw)/delta
+ self.iSolverAPI.solver.rho[i] = saverho[i]
+
+ # V
+ for iNo in range(nNs):
+ for iDim in range(nDim):
+ self.iSolverAPI.solver.U[iNo][iDim] += delta
+ blowing = self.__runViscous()
+ dRblw_v[:,iNo*nDim + iDim] = (blowing-saveBlw)/delta
+ self.iSolverAPI.solver.U[iNo][iDim] = savev[iNo*nDim + iDim]
+
+
+
+
+ def __runViscous(self):
+ self.iSolverAPI.imposeInvBC(0)
+ # Run viscous solver.
+ exitCode = self.vSolver.run(0)
+ # Impose blowing boundary condition in the inviscid solver.
+ self.iSolverAPI.imposeBlwVel()
+
+ blowing = np.zeros(len(self.iSolverAPI.iBnd[0].elemsCoord[:,0]))
+ for i in range(len(self.iSolverAPI.iBnd[0].elemsCoord[:,0])):
+ blowing[i] = self.iSolverAPI.blw[0].getU(i)
+ return blowing
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
- nRun = 0
- # Store blowing velocity
- saveBlw = self.__getBlowingVelocity()
- saveMach, saveRho, saveV = self.__getInviscidState()
- dRblw_dM = np.zeros((len(saveBlw), 0))
- deltaMach = 1e-4
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ # nRun = 0
+
+ # # Store blowing velocity
+ # saveBlw = self.__getBlowingVelocity()
+
+ # saveMach, saveRho, saveV = self.__getInviscidState()
+
+ # dRblw_dM = np.zeros((len(saveBlw), 0))
+ # deltaMach = 1e-4
+
+ # # Finite difference
+ # for i in range(len(saveMach)):
+ # self.iSolverAPI.iBnd[0].M[i] = saveMach[i] + deltaMach
+ # self.iSolverAPI.setViscousSolver(1)
+ # self.vSolver.run(1)
+ # nRun+=1
+ # self.iSolverAPI.imposeBlwVel()
+ # newBlowing = self.__getBlowingVelocity()
+ # dRblw_dM = np.column_stack((dRblw_dM, (newBlowing - saveBlw) / deltaMach))
+ # self.iSolverAPI.iBnd[0].M[i] = saveMach[i]
+
+ # dRblw_dRho = np.zeros((len(saveBlw), 0))
+ # deltaRho = 1e-4
- # Finite difference
- for i in range(len(saveMach)):
- self.iSolverAPI.iBnd[0].M[i] = saveMach[i] + deltaMach
- self.iSolverAPI.setViscousSolver(1)
- self.vSolver.run(1)
- nRun+=1
- self.iSolverAPI.imposeBlwVel()
- newBlowing = self.__getBlowingVelocity()
- dRblw_dM = np.column_stack((dRblw_dM, (newBlowing - saveBlw) / deltaMach))
- self.iSolverAPI.iBnd[0].M[i] = saveMach[i]
-
- dRblw_dRho = np.zeros((len(saveBlw), 0))
- deltaRho = 1e-8
+ # # Finite difference
+ # for i in range(len(saveRho)):
+ # self.iSolverAPI.iBnd[0].Rho[i] = saveRho[i] + deltaRho
+ # self.iSolverAPI.setViscousSolver(1)
+ # self.vSolver.run(1)
+ # nRun+=1
+ # self.iSolverAPI.imposeBlwVel()
+ # newBlowing = self.__getBlowingVelocity()
+ # dRblw_dRho = np.column_stack((dRblw_dRho, (newBlowing - saveBlw) / deltaRho))
+ # self.iSolverAPI.iBnd[0].Rho[i] = saveRho[i]
+
+ # dRblw_dv = np.zeros((len(saveBlw), 0))
+ # deltaV = 1e-4
- # Finite difference
- for i in range(len(saveRho)):
- self.iSolverAPI.iBnd[0].Rho[i] = saveRho[i] + deltaRho
- self.iSolverAPI.setViscousSolver(1)
- self.vSolver.run(1)
- nRun+=1
- self.iSolverAPI.imposeBlwVel()
- newBlowing = self.__getBlowingVelocity()
- dRblw_dRho = np.column_stack((dRblw_dRho, (newBlowing - saveBlw) / deltaRho))
- self.iSolverAPI.iBnd[0].Rho[i] = saveRho[i]
-
- dRblw_dv = np.zeros((len(saveBlw), 0))
- deltaV = 1e-4
+ # # Finite difference
+ # for j in range(self.iSolverAPI.solver.pbl.nDim):
+ # for i in range(len(saveV[:,0])):
+ # self.iSolverAPI.iBnd[0].V[i,j] = saveV[i,j] + deltaV
+ # self.iSolverAPI.setViscousSolver(1)
+ # self.vSolver.run(1)
+ # nRun+=1
+ # self.iSolverAPI.imposeBlwVel()
+ # newBlowing = self.__getBlowingVelocity()
+ # dRblw_dv = np.column_stack((dRblw_dv, (newBlowing - saveBlw) / deltaV))
+ # self.iSolverAPI.iBnd[0].V[i,j] = saveV[i,j]
- # Finite difference
- for j in range(self.iSolverAPI.solver.pbl.nDim):
- for i in range(len(saveV[:,0])):
- self.iSolverAPI.iBnd[0].V[i,j] = saveV[i,j] + deltaV
- self.iSolverAPI.setViscousSolver(1)
- self.vSolver.run(1)
- nRun+=1
- self.iSolverAPI.imposeBlwVel()
- newBlowing = self.__getBlowingVelocity()
- dRblw_dv = np.column_stack((dRblw_dv, (newBlowing - saveBlw) / deltaV))
- self.iSolverAPI.iBnd[0].V[i,j] = saveV[i,j]
-
- print('Total runs for finite difference', nRun)
- self.coupledAdjointSolver.setdRblw_M(dRblw_dM)
- self.coupledAdjointSolver.setdRblw_v(dRblw_dv)
- self.coupledAdjointSolver.setdRblw_rho(dRblw_dRho)
-
- def __getBlowingVelocity(self):
- blowingVelocity = np.zeros(0)
- for k in range(len(self.iSolverAPI.blw[0].tag.elems)):
- blowingVelocity = np.append(blowingVelocity, self.iSolverAPI.blw[0].getU(k))
- return blowingVelocity
+ # dRblw_dM[abs(dRblw_dM) < 1e-10] = 0.
+ # dRblw_dRho[abs(dRblw_dRho) < 1e-10] = 0.
+ # dRblw_dv[abs(dRblw_dv) < 1e-10] = 0.
+
+ # print('Total runs for finite difference', nRun)
+ # self.coupledAdjointSolver.setdRblw_M(dRblw_dM)
+ # self.coupledAdjointSolver.setdRblw_v(dRblw_dv)
+ # self.coupledAdjointSolver.setdRblw_rho(dRblw_dRho)
+
+ # def __getBlowingVelocity(self):
+ # blowingVelocity = np.zeros(0)
+ # for k in range(len(self.iSolverAPI.blw[0].tag.elems)):
+ # blowingVelocity = np.append(blowingVelocity, self.iSolverAPI.blw[0].getU(k))
+ # return blowingVelocity
- def __getInviscidState(self):
- V = np.zeros((len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]), 3))
- M = np.zeros(len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]))
- Rho = np.zeros(len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]))
- for n in range(2):
- for iRow, row in enumerate(self.iSolverAPI.iBnd[n].nodesCoord[:,3]):
- row=int(row)
- for iDim in range(3):
- V[iRow,iDim] = self.iSolverAPI.solver.U[row][iDim]
- M[iRow] = self.iSolverAPI.solver.M[row]
- Rho[iRow] = self.iSolverAPI.solver.rho[row]
- return M, Rho, V
+ # def __getInviscidState(self):
+ # V = np.zeros((len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]), 3))
+ # M = np.zeros(len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]))
+ # Rho = np.zeros(len(self.iSolverAPI.iBnd[0].nodesCoord[:,3]))
+ # for n in range(2):
+ # for iRow, row in enumerate(self.iSolverAPI.iBnd[n].nodesCoord[:,3]):
+ # row=int(row)
+ # for iDim in range(3):
+ # V[iRow,iDim] = self.iSolverAPI.solver.U[row][iDim]
+ # M[iRow] = self.iSolverAPI.solver.M[row]
+ # Rho[iRow] = self.iSolverAPI.solver.rho[row]
+ # return M, Rho, V
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 2b54e0b..b9bac11 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -28,6 +28,8 @@
#include "wGmres.h"
#include <Eigen/Sparse>
+#include<Eigen/SparseLU>
+#include <fstream>
#include <deque>
#include <algorithm>
@@ -48,10 +50,16 @@ CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::sh
sizeBlowing = iadjoint->getSizeBlowing();
- auto nDim = isol->pbl->nDim; // Number of dimensions of the problem
- auto nNs = isol->pbl->bBCs[0]->nodes.size(); // Number of surface nodes
- auto nNd = isol->pbl->msh->nodes.size(); // Number of domain nodes
- auto nEs = isol->pbl->bBCs[0]->uE.size(); // Number of elements on the surface
+ // Initalize all derivatives
+ this->reset();
+}
+
+void CoupledAdjoint::reset()
+{
+ size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
+ size_t nNs = isol->pbl->bBCs[0]->nodes.size(); // Number of surface nodes
+ size_t nNd = isol->pbl->msh->nodes.size(); // Number of domain nodes
+ size_t nEs = isol->pbl->bBCs[0]->uE.size(); // Number of elements on the surface
dRphi_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNd);
dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
@@ -146,31 +154,43 @@ void CoupledAdjoint::run()
std::copy(dCl_phi.data(), dCl_phi.data() + dCl_phi.size(), rhsCl.begin());
// Solution vector for lambdasCl
- Eigen::VectorXd lambdasCl(A_adjoint.cols());
+ Eigen::VectorXd lambdasCl = Eigen::VectorXd::Zero(A_adjoint.cols());
Eigen::Map<Eigen::VectorXd> lambdasCl_(lambdasCl.data(), lambdasCl.size());
// Solve coupled system
+ std::cout << "A_adjoint " << A_adjoint.rows() << " " << A_adjoint.cols() << " " << std::setprecision(10) << A_adjoint.norm() << std::endl;
+ std::cout << "rhsCl " << std::setprecision(10) << Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()).norm() << std::endl;
+ std::cout << "lambdasCl_ " << std::setprecision(10) << lambdasCl_.norm() << std::endl;
+
+ std::cout << "Before solver. NIter " << alinsol->getIterations() << ". Error " << alinsol->getError() << std::endl;
+ alinsol->mit = 500;
alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()), lambdasCl_);
+ std::cout << "After solver. NIter " << alinsol->getIterations() << ". Error " << alinsol->getError() << std::endl;
+ std::cout << "lambdasCl " << std::setprecision(10) << lambdasCl.norm() << std::endl;
Eigen::VectorXd lambdaClPhi = lambdasCl.block(0, 0, isol->pbl->msh->nodes.size(), 1);
+ std::cout << "lambdasClphi " << std::setprecision(10) << lambdaClPhi.norm() << std::endl;
+
// Compute total derivative dCl_dAoA
+ std::cout << "dCl_AoA (partial) " << std::setprecision(10) << dCl_AoA << std::endl;
+ std::cout << "dRphi_AoA " << std::setprecision(10) << dRphi_AoA.transpose().norm() << std::endl;
tdCl_AoA = dCl_AoA - dRphi_AoA.transpose()*lambdaClPhi;
/******** 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());
+ // 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());
+ // // 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);
+ // // 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;
+ // // Compute total derivative dCl_dAoA
+ // tdCd_AoA = dCd_AoA - dRphi_AoA.transpose()*lambdaCdPhi;
}
void CoupledAdjoint::gradientswrtInviscidFlow()
@@ -178,15 +198,11 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
//**** dRphi_phi ****//
dRphi_phi = iadjoint->getRu_U();
- // //**** dRM_phi, dRrho_phi, dRv_phi ****//
+ //**** dRM_phi, dRrho_phi, dRv_phi ****//
// auto pbl = isol->pbl;
- // //dRM_dPhi = iadjoint->getRM_U();
// // Finite difference
- // double deltaPhi = 1e-4; // perturbation
// std::vector<double> Phi_save = isol->phi; // Differential of the independant variable (phi)
- // std::vector<double> dM = std::vector<double>(pbl->bBCs[0]->nodes.size(), 0.0); // Differential of the Mach number
- // std::vector<double> dv = std::vector<double>(pbl->nDim*pbl->bBCs[0]->nodes.size(), 0.0); // Differential of the Mach number
- // std::vector<double> drho = std::vector<double>(pbl->bBCs[0]->nodes.size(), 0.0); // Differential of the density
+ // double deltaPhi = 1e-6; // perturbation
// std::vector<Eigen::Triplet<double>> tripletsdM;
// std::vector<Eigen::Triplet<double>> tripletsdrho;
@@ -195,45 +211,42 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
// for (auto n : pbl->msh->nodes){
// Phi_save[n->row] = isol->phi[n->row] + deltaPhi;
- // std::fill(dM.begin(), dM.end(), 0.0);
- // std::fill(drho.begin(), drho.end(), 0.0);
-
- // // Recompute Mach number on surface nodes.
+ // // Recompute the quantities on surface nodes.
// size_t jj = 0;
// for (auto srfNode : pbl->bBCs[0]->nodes){
- // // Average of the Mach on the elements adjacent to the considered node.
+
+ // double dM = 0.; double drho = 0.;
+ // Eigen::Vector3d dv = Eigen::Vector3d::Zero();
+ // // Average of the quantity on the elements adjacent to the considered node.
// for (auto e : pbl->fluid->neMap[srfNode]){
- // dM[jj] += pbl->fluid->mach->eval(*e, Phi_save, 0);
- // drho[jj] += pbl->fluid->rho->eval(*e, Phi_save, 0);
+ // dM += pbl->fluid->mach->eval(*e, Phi_save, 0);
+ // drho += pbl->fluid->rho->eval(*e, Phi_save, 0);
// Eigen::VectorXd grad = e->computeGradient(Phi_save, 0);
- // for (size_t i = 0; i < grad.size(); ++i)
- // dv[jj*pbl->nDim + i] += grad(i);
- // }
- // dM[jj] /= pbl->fluid->neMap[srfNode].size();
- // drho[jj] /= pbl->fluid->neMap[srfNode].size();
- // for (size_t i = 0; i < pbl->nDim; ++i){
- // dv[jj*pbl->nDim + i] /= pbl->fluid->neMap[srfNode].size();
- // tripletsdv.push_back(Eigen::Triplet<double>(jj*pbl->nDim + i, n->row, dv[jj*pbl->nDim + i] - isol->U[srfNode->row][i]/deltaPhi));
+ // for (int i = 0; i < grad.size(); ++i)
+ // dv(i) += grad(i);
// }
+ // dM /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ // drho /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ // dv /= static_cast<double>(pbl->fluid->neMap[srfNode].size());
+
// // Form triplets
- // tripletsdM.push_back(Eigen::Triplet<double>(jj, n->row, (dM[jj] - isol->M[srfNode->row])/deltaPhi));
- // tripletsdrho.push_back(Eigen::Triplet<double>(jj, n->row, (drho[jj] - isol->rho[srfNode->row])/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 (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));
// jj++;
// }
// // Reset phi
// Phi_save[n->row] = isol->phi[n->row];
// }
-
// // Fill matrices with gradients
// dRM_phi.setFromTriplets(tripletsdM.begin(), tripletsdM.end());
// dRrho_phi.setFromTriplets(tripletsdrho.begin(), tripletsdrho.end());
// dRv_phi.setFromTriplets(tripletsdv.begin(), tripletsdv.end());
// // Remove zeros
- // dRM_phi.prune(0.);
- // dRrho_phi.prune(0.);
- // dRv_phi.prune(0.);
+ // dRM_phi.prune(0.); dRrho_phi.prune(0.); dRv_phi.prune(0.);
- // Objective functionss
+ //**** dCl_phi, dCd_phi ****//
std::vector<std::vector<double>> dCx_phi = iadjoint->computeGradientCoefficientsFlow();
dCl_phi = Eigen::VectorXd::Map(dCx_phi[0].data(), dCx_phi[0].size());
dCd_phi = Eigen::VectorXd::Map(dCx_phi[1].data(), dCx_phi[1].size());
@@ -269,7 +282,7 @@ void CoupledAdjoint::gradientswrtVelocity(){
}
void CoupledAdjoint::gradientswrtBlowingVelocity(){
- //dRphi_blw = iadjoint->getRu_Blw();
+ // dRphi_blw = iadjoint->getRu_Blw();
// All others are zero.
}
@@ -354,4 +367,115 @@ void CoupledAdjoint::transposeMatrices(){
dRblw_rho = dRblw_rho.transpose();
dRblw_v = dRblw_v.transpose();
dRblw_blw = dRblw_blw.transpose();
+
+ // writeMatrixToFile(dRphi_phi, "dRphi_phi.txt");
+ // writeMatrixToFile(dRphi_M, "dRphi_M.txt");
+ // writeMatrixToFile(dRphi_rho, "dRphi_rho.txt");
+ // writeMatrixToFile(dRphi_v, "dRphi_v.txt");
+ // writeMatrixToFile(dRphi_blw, "dRphi_blw.txt");
+
+ // writeMatrixToFile(dRM_phi, "dRM_phi.txt");
+ // writeMatrixToFile(dRM_M, "dRM_M.txt");
+ // writeMatrixToFile(dRM_rho, "dRM_rho.txt");
+ // writeMatrixToFile(dRM_v, "dRM_v.txt");
+ // writeMatrixToFile(dRM_blw, "dRM_blw.txt");
+
+ // writeMatrixToFile(dRrho_phi, "dRrho_phi.txt");
+ // writeMatrixToFile(dRrho_M, "dRrho_M.txt");
+ // writeMatrixToFile(dRrho_rho, "dRrho_rho.txt");
+ // writeMatrixToFile(dRrho_v, "dRrho_v.txt");
+ // writeMatrixToFile(dRrho_blw, "dRrho_blw.txt");
+
+ // writeMatrixToFile(dRv_phi, "dRv_phi.txt");
+ // writeMatrixToFile(dRv_M, "dRv_M.txt");
+ // writeMatrixToFile(dRv_rho, "dRv_rho.txt");
+ // writeMatrixToFile(dRv_v, "dRv_v.txt");
+ // writeMatrixToFile(dRv_blw, "dRv_blw.txt");
+
+ // writeMatrixToFile(dRblw_phi, "dRblw_phi.txt");
+ // writeMatrixToFile(dRblw_M, "dRblw_M.txt");
+ // writeMatrixToFile(dRblw_rho, "dRblw_rho.txt");
+ // writeMatrixToFile(dRblw_v, "dRblw_v.txt");
+ // writeMatrixToFile(dRblw_blw, "dRblw_blw.txt");
+
+ // // (Cout all norms)
+
+ std::cout << "norm(dRphi_phi) " << dRphi_phi.norm() << std::endl;
+ std::cout << "norm(dRphi_M) " << dRphi_M.norm() << std::endl;
+ std::cout << "norm(dRphi_rho) " << dRphi_rho.norm() << std::endl;
+ std::cout << "norm(dRphi_v) " << dRphi_v.norm() << std::endl;
+ std::cout << "norm(dRphi_blw) " << dRphi_blw.norm() << std::endl;
+
+ std::cout << "norm(dRM_phi) " << dRM_phi.norm() << std::endl;
+ std::cout << "norm(dRM_M) " << dRM_M.norm() << std::endl;
+ std::cout << "norm(dRM_rho) " << dRM_rho.norm() << std::endl;
+ std::cout << "norm(dRM_v) " << dRM_v.norm() << std::endl;
+ std::cout << "norm(dRM_blw) " << dRM_blw.norm() << std::endl;
+
+ std::cout << "norm(dRrho_phi) " << dRrho_phi.norm() << std::endl;
+ std::cout << "norm(dRrho_M) " << dRrho_M.norm() << std::endl;
+ std::cout << "norm(dRrho_rho) " << dRrho_rho.norm() << std::endl;
+ std::cout << "norm(dRrho_v) " << dRrho_v.norm() << std::endl;
+ std::cout << "norm(dRrho_blw) " << dRrho_blw.norm() << std::endl;
+
+ std::cout << "norm(dRv_phi) " << dRv_phi.norm() << std::endl;
+ std::cout << "norm(dRv_M) " << dRv_M.norm() << std::endl;
+ std::cout << "norm(dRv_rho) " << dRv_rho.norm() << std::endl;
+ std::cout << "norm(dRv_v) " << dRv_v.norm() << std::endl;
+ std::cout << "norm(dRv_blw) " << dRv_blw.norm() << std::endl;
+
+ std::cout << "norm(dRblw_phi) " << dRblw_phi.norm() << std::endl;
+ std::cout << "norm(dRblw_M) " << dRblw_M.norm() << std::endl;
+ std::cout << "norm(dRblw_rho) " << dRblw_rho.norm() << std::endl;
+ std::cout << "norm(dRblw_v) " << dRblw_v.norm() << std::endl;
+ std::cout << "norm(dRblw_blw) " << dRblw_blw.norm() << std::endl;
+
+ // // Print all matrices after transpose
+
+ // std::cout << "dRphi_phi " << dRphi_phi.rows() << " " << dRphi_phi.cols() << "\n" << dRphi_phi << std::endl;
+ // std::cout << "dRphi_M " << dRphi_M.rows() << " " << dRphi_M.cols() << "\n" << dRphi_M << std::endl;
+ // std::cout << "dRphi_rho " << dRphi_rho.rows() << " " << dRphi_rho.cols() << "\n" << dRphi_rho << std::endl;
+ // std::cout << "dRphi_v " << dRphi_v.rows() << " " << dRphi_v.cols() << "\n" << dRphi_v << std::endl;
+ // std::cout << "dRphi_blw " << dRphi_blw.rows() << " " << dRphi_blw.cols() << "\n" << dRphi_blw << std::endl;
+
+ // std::cout << "dRM_phi " << dRM_phi.rows() << " " << dRM_phi.cols() << "\n" << dRM_phi << std::endl;
+ // std::cout << "dRM_M " << dRM_M.rows() << " " << dRM_M.cols() << "\n" << dRM_M << std::endl;
+ // std::cout << "dRM_rho " << dRM_rho.rows() << " " << dRM_rho.cols() << "\n" << dRM_rho << std::endl;
+ // std::cout << "dRM_v " << dRM_v.rows() << " " << dRM_v.cols() << "\n" << dRM_v << std::endl;
+ // std::cout << "dRM_blw " << dRM_blw.rows() << " " << dRM_blw.cols() << "\n" << dRM_blw << std::endl;
+
+ // std::cout << "dRrho_phi " << dRrho_phi.rows() << " " << dRrho_phi.cols() << "\n" << dRrho_phi << std::endl;
+ // std::cout << "dRrho_M " << dRrho_M.rows() << " " << dRrho_M.cols() << "\n" << dRrho_M << std::endl;
+ // std::cout << "dRrho_rho " << dRrho_rho.rows() << " " << dRrho_rho.cols() << "\n" << dRrho_rho << std::endl;
+ // std::cout << "dRrho_v " << dRrho_v.rows() << " " << dRrho_v.cols() << "\n" << dRrho_v << std::endl;
+ // std::cout << "dRrho_blw " << dRrho_blw.rows() << " " << dRrho_blw.cols() << "\n" << dRrho_blw << std::endl;
+
+ // std::cout << "dRv_phi " << dRv_phi.rows() << " " << dRv_phi.cols() << "\n" << dRv_phi << std::endl;
+ // std::cout << "dRv_M " << dRv_M.rows() << " " << dRv_M.cols() << "\n" << dRv_M << std::endl;
+ // std::cout << "dRv_rho " << dRv_rho.rows() << " " << dRv_rho.cols() << "\n" << dRv_rho << std::endl;
+ // std::cout << "dRv_v " << dRv_v.rows() << " " << dRv_v.cols() << "\n" << dRv_v << std::endl;
+ // std::cout << "dRv_blw " << dRv_blw.rows() << " " << dRv_blw.cols() << "\n" << dRv_blw << std::endl;
+
+ // std::cout << "dRblw_phi " << dRblw_phi.rows() << " " << dRblw_phi.cols() << "\n" << dRblw_phi << std::endl;
+ // std::cout << "dRblw_M " << dRblw_M.rows() << " " << dRblw_M.cols() << "\n" << dRblw_M << std::endl;
+ // std::cout << "dRblw_rho " << dRblw_rho.rows() << " " << dRblw_rho.cols() << "\n" << dRblw_rho << std::endl;
+ // std::cout << "dRblw_v " << dRblw_v.rows() << " " << dRblw_v.cols() << "\n" << dRblw_v << std::endl;
+ // std::cout << "dRblw_blw " << dRblw_blw.rows() << " " << dRblw_blw.cols() << "\n" << dRblw_blw << std::endl;
+
+}
+
+void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename) {
+ std::ofstream file(filename);
+ if (file.is_open()) {
+ for (int i = 0; i < matrix.rows(); ++i) {
+ for (int j = 0; j < matrix.cols(); ++j) {
+ file << matrix(i, j);
+ if (j != matrix.cols() - 1) {
+ file << " ";
+ }
+ }
+ file << "\n";
+ }
+ file.close();
+ }
}
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 93973c7..8b2a61e 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -87,6 +87,7 @@ public:
virtual ~CoupledAdjoint () {std::cout << "~blast::CoupledAdjoint()\n";};
void run();
+ void reset();
void setdRblw_M(std::vector<std::vector<double>> dRblw_dM);
void setdRblw_v(std::vector<std::vector<double>> dRblw_dv);
@@ -104,6 +105,7 @@ private:
void gradientwrtAoA();
void transposeMatrices();
+ void writeMatrixToFile(const Eigen::MatrixXd& matrix, const std::string& filename);
};
} // namespace blast
#endif // DDARTADJOINT
diff --git a/blast/tests/dart/adjointVII.py b/blast/tests/dart/adjointVII.py
index a689678..5a9cf09 100644
--- a/blast/tests/dart/adjointVII.py
+++ b/blast/tests/dart/adjointVII.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' : 2, 'msTe' : 0.05, 'msLe' : 0.05}, # parameters for input file model
+ 'Pars' : {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}, # parameters for input file model
'Dim' : 2, # problem dimension
'Format' : 'gmsh', # save format (vtk or gmsh)
# Markers
@@ -135,25 +135,180 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
- aeroCoeffs = coupler.run()
+ #aeroCoeffs = coupler.run()
tms['solver'].stop()
- #iSolverAPI.run()
+ iSolverAPI.run()
tms['FDDerivatives'].start()
- #pyAdjoint.computeDerivatives()
tms['FDDerivatives'].stop()
tms['CoupledAdjoint'].start()
iSolverAPI.adjointSolver.run()
tms['CoupledAdjoint'].stop()
dCl_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dClAoa
dCd_dAlpha_DartAdjoint = iSolverAPI.adjointSolver.dCdAoa
+ #pyAdjoint.computeDerivatives()
coupledAdjointSolver.run()
+ print(coupledAdjointSolver.tdCl_AoA)
+
+ quit()
+
+ print('loop')
+ dRphi_phi = []
+ dRphi_M = []
+ dRphi_rho = []
+ dRphi_v = []
+ dRphi_blw = []
+ dRM_phi = []
+ dRM_M = []
+ dRM_rho = []
+ dRM_v = []
+ dRM_blw = []
+ dRrho_phi = []
+ dRrho_M = []
+ dRrho_rho = []
+ dRrho_v = []
+ dRrho_blw = []
+ dRv_phi = []
+ dRv_M = []
+ dRv_rho = []
+ dRv_v = []
+ dRv_blw = []
+ dRblw_phi = []
+ dRblw_M = []
+ dRblw_rho = []
+ dRblw_v = []
+ dRblw_blw = []
+ A_adjoint = []
+ rhsCl = []
+ for i in range(10):
+ coupledAdjointSolver.reset()
+ #pyAdjoint.computeDerivatives()
+ coupledAdjointSolver.run()
+ dRphi_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_phi.txt'))
+ dRphi_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_M.txt'))
+ dRphi_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_rho.txt'))
+ dRphi_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_v.txt'))
+ dRphi_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRphi_blw.txt'))
+ dRM_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_phi.txt'))
+ dRM_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_M.txt'))
+ dRM_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_rho.txt'))
+ dRM_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_v.txt'))
+ dRM_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRM_blw.txt'))
+ dRrho_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_phi.txt'))
+ dRrho_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_M.txt'))
+ dRrho_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_rho.txt'))
+ dRrho_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_v.txt'))
+ dRrho_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRrho_blw.txt'))
+ dRv_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_phi.txt'))
+ dRv_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_M.txt'))
+ dRv_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_rho.txt'))
+ dRv_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_v.txt'))
+ dRv_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRv_blw.txt'))
+ dRblw_phi.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_phi.txt'))
+ dRblw_M.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_M.txt'))
+ dRblw_rho.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_rho.txt'))
+ dRblw_v.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_v.txt'))
+ dRblw_blw.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/dRblw_blw.txt'))
+ A_adjoint.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/A_adjoint.txt'))
+ rhsCl.append(np.loadtxt('/Users/pauldechamps/lab/softwares/blaster/workspace/blast_tests_dart_adjointVII/rhsCl.txt'))
+ print(i, coupledAdjointSolver.tdCl_AoA)
+ np.set_printoptions(threshold=np.inf)
+
+ #print(dRblw_M[0] - dRblw_M[1])
+
+ for i in range(len(dRphi_phi)-1):
+ if np.any(dRphi_phi[i+1] - dRphi_phi[i] != 0):
+ print('dRphi_phi')
+ break
+ if np.any(dRphi_M[i+1] - dRphi_M[i] != 0):
+ print('dRphi_M')
+ break
+ if np.any(dRphi_rho[i+1] - dRphi_rho[i] != 0):
+ print('dRphi_rho')
+ break
+ if np.any(dRphi_v[i+1] - dRphi_v[i] != 0):
+ print('dRphi_v')
+ break
+ if np.any(dRphi_blw[i+1] - dRphi_blw[i] != 0):
+ print('dRphi_blw')
+ break
+ if np.any(dRM_phi[i+1] - dRM_phi[i] != 0):
+ print('dRM_phi')
+ break
+ if np.any(dRM_M[i+1] - dRM_M[i] != 0):
+ print('dRM_M')
+ break
+ if np.any(dRM_rho[i+1] - dRM_rho[i] != 0):
+ print('dRM_rho')
+ break
+ if np.any(dRM_v[i+1] - dRM_v[i] != 0):
+ print('dRM_v')
+ break
+ if np.any(dRM_blw[i+1] - dRM_blw[i] != 0):
+ print('dRM_blw')
+ break
+ if np.any(dRrho_phi[i+1] - dRrho_phi[i] != 0):
+ print('dRrho_phi')
+ break
+ if np.any(dRrho_M[i+1] - dRrho_M[i] != 0):
+ print('dRrho_M')
+ break
+ if np.any(dRrho_rho[i+1] - dRrho_rho[i] != 0):
+ print('dRrho_rho')
+ break
+ if np.any(dRrho_v[i+1] - dRrho_v[i] != 0):
+ print('dRrho_v')
+ break
+ if np.any(dRrho_blw[i+1] - dRrho_blw[i] != 0):
+ print('dRrho_blw')
+ break
+ if np.any(dRv_phi[i+1] - dRv_phi[i] != 0):
+ print('dRv_phi')
+ break
+ if np.any(dRv_M[i+1] - dRv_M[i] != 0):
+ print('dRv_M')
+ break
+ if np.any(dRv_rho[i+1] - dRv_rho[i] != 0):
+ print('dRv_rho')
+ break
+ if np.any(dRv_v[i+1] - dRv_v[i] != 0):
+ print('dRv_v')
+ break
+ if np.any(dRv_blw[i+1] - dRv_blw[i] != 0):
+ print('dRv_blw')
+ break
+ if np.any(dRblw_phi[i+1] - dRblw_phi[i] != 0):
+ print('dRblw_phi')
+ break
+ if np.any(dRblw_M[i+1] - dRblw_M[i] >= 1e-8):
+ print('dRblw_M')
+ print(np.any(dRblw_M[i+1] - dRblw_M[i] >= 1e-8))
+ #break
+ if np.any(dRblw_rho[i+1] - dRblw_rho[i] >= 1e-8):
+ print('dRblw_rho')
+ #break
+ if np.any(dRblw_v[i+1] - dRblw_v[i] >= 1e-8):
+ print('dRblw_v')
+ #break
+ if np.any(dRblw_blw[i+1] - dRblw_blw[i] >= 1e-8):
+ print('dRblw_blw')
+ break
+ if np.any(A_adjoint[i+1] - A_adjoint[i] != 0):
+ print('A_adjoint')
+ break
+ if np.any(rhsCl[i+1] - rhsCl[i] != 0):
+ print('rhsCl')
+ break
+ for i in rhsCl:
+ print(i)
print("sucess", coupledAdjointSolver)
+ quit()
dalpha = 1e-4
aoaV = [alpha-dalpha, alpha+dalpha]
clV = []
cdV = []
+
for aoa in aoaV:
coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
iSolverAPI.argOut['pbl'].update(aoa)
@@ -177,17 +332,19 @@ def main():
dCd_dalpha_FD_dart = (cdV[-1] - cdV[0]) / (2*dalpha)
print(' ------------------- RESULTS -------------------')
- print('CL')
- print('dCl_dalpha_adjointDart', dCl_dAlpha_DartAdjoint)
- print('dCl_dalpha_FD', dCl_dalpha_FD_dart)
- print('dCl_dalpha_FD', dCl_dalpha_FD_coupled)
- print('dCl_dalpha_coupledAdjoint', coupledAdjointSolver.tdCl_AoA)
-
- print('CD')
- print('dCd_dalpha_adjointDart', dCd_dAlpha_DartAdjoint)
- print('dCd_dalpha_FD', dCd_dalpha_FD_dart)
- print('dCd_dalpha_FD', dCd_dalpha_FD_coupled)
- print('dCd_dalpha_coupledAdjoint', coupledAdjointSolver.tdCd_AoA)
+ print('********************** CL **********************')
+ print('DART FD', dCl_dalpha_FD_dart)
+ print('DART ADJOINT', dCl_dAlpha_DartAdjoint)
+ print('------------------------------------------------')
+ print('COUPLED FD', dCl_dalpha_FD_coupled)
+ print('COUPLED ADJOINT', coupledAdjointSolver.tdCl_AoA)
+
+ print('********************** CD **********************')
+ print('DART FD', dCd_dalpha_FD_dart)
+ print('DART ADJOINT', dCd_dAlpha_DartAdjoint)
+ print('------------------------------------------------')
+ print('COUPLED FD ', dCd_dalpha_FD_coupled)
+ print('COUPLED ADJOINT', coupledAdjointSolver.tdCd_AoA)
print(tms)
--
GitLab