(feat) Adjoint: Validating inviscid part.

Bug with linear solver?

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

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();
+    }
 }

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

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)