From fd2279f39229900cec6a745d967bb500bb624534 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Sat, 24 Feb 2024 13:51:08 +0100
Subject: [PATCH] (feat) Build adjoint system and solve it
Tested only for phi influence (same results as dart). The rest is implemented but untested until the bug for dirichlet bc is fixed.
---
blast/_src/blastw.i | 2 +
blast/interfaces/dart/DartAdjointInterface.py | 83 ++++
blast/src/CMakeLists.txt | 8 +-
blast/src/DCoupledAdjoint.cpp | 438 +++++++++++-------
blast/src/DCoupledAdjoint.h | 23 +-
blast/src/DDriver.cpp | 2 +-
blast/tests/dart/adjointVII.py | 122 ++++-
7 files changed, 489 insertions(+), 189 deletions(-)
create mode 100644 blast/interfaces/dart/DartAdjointInterface.py
diff --git a/blast/_src/blastw.i b/blast/_src/blastw.i
index ac1b2e6..732a113 100644
--- a/blast/_src/blastw.i
+++ b/blast/_src/blastw.i
@@ -58,3 +58,5 @@ threads="1"
%shared_ptr(blast::CoupledAdjoint);
%include "DCoupledAdjoint.h"
+%immutable blast::CoupledAdjoint::tdCl_AoA;
+%immutable blast::CoupledAdjoint::tdCd_AoA;
diff --git a/blast/interfaces/dart/DartAdjointInterface.py b/blast/interfaces/dart/DartAdjointInterface.py
new file mode 100644
index 0000000..cb3d632
--- /dev/null
+++ b/blast/interfaces/dart/DartAdjointInterface.py
@@ -0,0 +1,83 @@
+import numpy as np
+
+class DartAdjointInterface():
+ def __init__(self, _coupledAdjointSolver, _iSolverAPI, _vSolver):
+ self.coupledAdjointSolver = _coupledAdjointSolver
+ self.iSolverAPI = _iSolverAPI
+ self.vSolver = _vSolver
+
+ def computeDerivatives(self):
+
+ 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-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 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
+
+ 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/CMakeLists.txt b/blast/src/CMakeLists.txt
index 394f4cb..f6964d9 100644
--- a/blast/src/CMakeLists.txt
+++ b/blast/src/CMakeLists.txt
@@ -18,9 +18,13 @@ FILE(GLOB SRCS *.h *.cpp *.inl *.hpp)
ADD_LIBRARY(blast SHARED ${SRCS})
MACRO_DebugPostfix(blast)
-TARGET_INCLUDE_DIRECTORIES(blast PUBLIC ${PROJECT_SOURCE_DIR}/blast/src)
+TARGET_INCLUDE_DIRECTORIES(blast PUBLIC ${PROJECT_SOURCE_DIR}/blast/src
+ ${PROJECT_SOURCE_DIR}/modules/dartflo/ext/amfe/tbox/src
+ ${PROJECT_SOURCE_DIR}/modules/dartflo/ext/amfe/fwk/src
+ ${PROJECT_SOURCE_DIR}/modules/dartflo/dart/src
+ ${PYTHON_INCLUDE_PATH})
-TARGET_LINK_LIBRARIES(blast tbox)
+TARGET_LINK_LIBRARIES(blast dart tbox fwk)
INSTALL(TARGETS blast DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 34db14b..2b54e0b 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/blast/src/DCoupledAdjoint.cpp
@@ -15,12 +15,12 @@
*/
#include "DCoupledAdjoint.h"
-#include "../../modules/dartflo/dart/src/wAdjoint.h"
-#include "../../modules/dartflo/dart/src/wNewton.h"
-#include "../../modules/dartflo/dart/src/wSolver.h"
-#include "../../modules/dartflo/dart/src/wProblem.h"
-#include "../../modules/dartflo/dart/src/wBlowing.h"
-#include "../../modules/dartflo/dart/src/wFluid.h"
+#include "wAdjoint.h"
+#include "wNewton.h"
+#include "wSolver.h"
+#include "wProblem.h"
+#include "wBlowing.h"
+#include "wFluid.h"
#include "wMshData.h"
#include "wNode.h"
@@ -45,136 +45,132 @@ CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint, std::sh
}
else
alinsol = isol->linsol;
-
+
sizeBlowing = iadjoint->getSizeBlowing();
-
- dCldAoA = 0.0;
- dCddAoA = 0.0;
+
+ 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
+
+ dRphi_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNd);
+ dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
+ dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNd);
+ dRv_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNd);
+ dRblw_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNd);
+
+ dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
+ dRM_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs);
+
+ dRphi_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
+ dRM_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRrho_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
+ dRv_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nNs);
+ dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nNs);
+
+ dRphi_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim*nNs);
+ dRM_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
+ dRrho_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim*nNs);
+ dRv_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nDim*nNs);
+ dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim*nNs);
+
+ dRphi_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nEs);
+ dRM_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
+ dRrho_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
+ dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim*nNs, nEs);
+ dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nEs);
+
+ dCl_phi = Eigen::VectorXd::Zero(nNd);
+ dCd_phi = Eigen::VectorXd::Zero(nNd);
+
+ dCl_AoA = 0.0;
+ dCd_AoA = 0.0;
}
void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint)
{
- // Determine the sizes of the smaller matrices
- int rows = matrixAdjoint.rows();
- int cols = matrixAdjoint.cols();
-
// Create a list of triplets for the larger matrix
std::vector<Eigen::Triplet<double>> triplets;
// Iterate over the rows and columns of the vector
- int rowOffset = 0;
- for (const auto& row : matrices) {
- int colOffset = 0;
- for (const auto& matrix : row) {
- // Add the triplets of the matrix to the list of triplets for the larger matrix
- for (int k = 0; k < matrix->outerSize(); ++k) {
- for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(*matrix, k); it; ++it) {
- triplets.push_back(Eigen::Triplet<double>(it.row() + rowOffset, it.col() + colOffset, it.value()));
- }
+int rowOffset = 0;
+for (const auto& row : matrices) {
+ int colOffset = 0;
+ for (const auto& matrix : row) {
+ // Add the triplets of the matrix to the list of triplets for the larger matrix
+ for (int k = 0; k < matrix->outerSize(); ++k) {
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(*matrix, k); it; ++it) {
+ triplets.push_back(Eigen::Triplet<double>(it.row() + rowOffset, it.col() + colOffset, it.value()));
}
- colOffset += matrix->cols();
}
- rowOffset += row[0]->rows();
+ colOffset += matrix->cols();
}
-
+ rowOffset += row[0]->rows();
+}
// Create a new matrix from the deque of triplets
matrixAdjoint.setFromTriplets(triplets.begin(), triplets.end());
+ matrixAdjoint.prune(0.);
+ matrixAdjoint.makeCompressed();
}
void CoupledAdjoint::run()
-{
- Eigen::SparseMatrix<double, Eigen::RowMajor> A = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> B = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> C = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> D = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> E = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> F = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> G = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> H = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
- Eigen::SparseMatrix<double, Eigen::RowMajor> I = Eigen::SparseMatrix<double, Eigen::RowMajor>(3,3);
-
- std::deque<Eigen::Triplet<double>> T;
- T.push_back(Eigen::Triplet<double>(0, 0, 1.0));
- T.push_back(Eigen::Triplet<double>(1, 1, 2.0));
- A.setFromTriplets(T.begin(), T.end());
- //std::cout << "A\n" <<A << std::endl;
- T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 3.0));
- T.push_back(Eigen::Triplet<double>(0, 1, 2.0));
- B.setFromTriplets(T.begin(), T.end());
- //std::cout << "B\n" <<B << std::endl;
- /*T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 0.0));
- T.push_back(Eigen::Triplet<double>(1, 1, 0.0));
- C.setFromTriplets(T.begin(), T.end());*/
- std::cout << "C\n" << C << std::endl;
- T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
- T.push_back(Eigen::Triplet<double>(2, 1, 1.0));
- D.setFromTriplets(T.begin(), T.end());
- //std::cout << "D\n" << D << std::endl;
-
- T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
- T.push_back(Eigen::Triplet<double>(1, 1, 1.0));
- E.setFromTriplets(T.begin(), T.end());
- //std::cout << "E\n" << E << std::endl;
-
- T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
- T.push_back(Eigen::Triplet<double>(1, 1, 1.0));
- F.setFromTriplets(T.begin(), T.end());
- //std::cout << "F\n" << F << std::endl;
-
- T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
- T.push_back(Eigen::Triplet<double>(1, 1, 4.0));
- G.setFromTriplets(T.begin(), T.end());
- //std::cout << "G\n" << G << std::endl;
-
- T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 2.0));
- T.push_back(Eigen::Triplet<double>(1, 1, 1.0));
- H.setFromTriplets(T.begin(), T.end());
- //std::cout << "H\n" << H << std::endl;
-
- T.clear();
- T.push_back(Eigen::Triplet<double>(0, 0, 5.0));
- T.push_back(Eigen::Triplet<double>(1, 2, 1.0));
- I.setFromTriplets(T.begin(), T.end());
- //std::cout << "I\n" << I << std::endl;
-
- std::cout << "gradientswrtInviscidFlow\n";
- dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->nodes.size(), isol->pbl->msh->nodes.size());
- dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->nodes.size(), isol->pbl->msh->nodes.size());
+{
gradientswrtInviscidFlow();
- std::cout << " coucou " << std::endl;
- std::cout << "dRM_phi\n" << dRM_phi.norm() << std::endl;
- std::cout << "dRrho_phi\n" << dRrho_phi.norm() << std::endl;
+ gradientswrtBlowingVelocity();
+ transposeMatrices();
+ gradientwrtAoA();
+
+ // Adjoint matrix
// Store pointers to the matrices in a 2D vector
std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> matrices = {
- {&A, &B, &C, &D, &E},
- {&D, &E, &F, &I, &G},
- {&G, &H, &I, &A, &B},
- {&B, &E, &F, &I, &G},
- {&B, &A, &I, &C, &B}
+ {&dRphi_phi, &dRM_phi, &dRv_phi, &dRrho_phi, &dRblw_phi},
+ {&dRphi_M, &dRM_M, &dRv_M, &dRrho_M, &dRblw_M},
+ {&dRphi_rho, &dRM_rho, &dRv_rho, &dRrho_rho, &dRblw_rho},
+ {&dRphi_v, &dRM_v, &dRv_v, &dRrho_v, &dRblw_v},
+ {&dRphi_blw, &dRM_blw, &dRv_blw, &dRrho_blw, &dRblw_blw}
};
- // Create a larger sparse matrix and set it from the list of triplets
int rows = 0; int cols = 0;
- for (const auto& row : matrices) {
- rows += row[0]->rows(); // All matrices in a row have the same number of rows
- cols = std::max(cols, static_cast<int>(row.size() * row[0]->cols())); // All matrices in a column have the same number of columns
- }
+ for (const auto& row : matrices) rows += row[0]->rows(); // All matrices in a row have the same number of rows
+ for (const auto& mat1 : matrices[0]) cols += mat1->cols(); // All matrices in a column have the same number of columns
Eigen::SparseMatrix<double, Eigen::RowMajor> A_adjoint(rows, cols);
-
buildAdjointMatrix(matrices, A_adjoint);
- std::cout << "largeMatrix\n" << A_adjoint << std::endl;
+ /******** CL ********/
+ // Build right hand side for Cl
+ std::vector<double> rhsCl(A_adjoint.cols(), 0.0);
+ std::copy(dCl_phi.data(), dCl_phi.data() + dCl_phi.size(), rhsCl.begin());
+
+ // Solution vector for lambdasCl
+ Eigen::VectorXd lambdasCl(A_adjoint.cols());
+ Eigen::Map<Eigen::VectorXd> lambdasCl_(lambdasCl.data(), lambdasCl.size());
// Solve coupled system
- // alinsol->compute(AdjointMatrix.transpose(), Eigen::Map<Eigen::VectorXd>(dU_.data(), dU_.size()), lambdas);
+ alinsol->compute(A_adjoint, Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size()), lambdasCl_);
+ Eigen::VectorXd lambdaClPhi = lambdasCl.block(0, 0, isol->pbl->msh->nodes.size(), 1);
+
+ // Compute total derivative dCl_dAoA
+ 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());
+
+ // 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;
}
void CoupledAdjoint::gradientswrtInviscidFlow()
@@ -182,96 +178,180 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
//**** dRphi_phi ****//
dRphi_phi = iadjoint->getRu_U();
- //**** 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> drho = std::vector<double>(pbl->bBCs[0]->nodes.size(), 0.0); // Differential of the density
-
- std::vector<Eigen::Triplet<double>> tripletsdM;
- std::vector<Eigen::Triplet<double>> tripletsdrho;
-
- 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.
- size_t jj = 0;
- for (auto srfNode : pbl->bBCs[0]->nodes){
- // Average of the Mach 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[jj] /= pbl->fluid->neMap[srfNode].size();
- drho[jj] /= 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));
- 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());
- // Remove zeros
- dRM_phi.prune(0.);
- dRrho_phi.prune(0.);
-
- // // Objective functionss
- // dCl_dPhi =
- // dCd_dPhi =
+ // //**** 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
+
+ // std::vector<Eigen::Triplet<double>> tripletsdM;
+ // std::vector<Eigen::Triplet<double>> tripletsdrho;
+ // std::vector<Eigen::Triplet<double>> tripletsdv;
+
+ // 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.
+ // size_t jj = 0;
+ // for (auto srfNode : pbl->bBCs[0]->nodes){
+ // // Average of the Mach 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);
+ // 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));
+ // }
+ // // 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));
+ // 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.);
+
+ // Objective functionss
+ 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());
}
void CoupledAdjoint::gradientswrtMachNumber(){
- // dRphi_dM = 0
- // dRM_dM =
- // dRrho_dM = 0
- // dRv_dM = 0
- // dRblw_dM =
-
- // // Objective functions
- // dCl_dM =
- // dCd_dM =
+ // dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->msh->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->nDim*isol->pbl->bBCs[0]->nodes.size(), isol->pbl->bBCs[0]->nodes.size());
+ // dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(isol->pbl->bBCs[0]->uE.size(), isol->pbl->bBCs[0]->nodes.size());
}
void CoupledAdjoint::gradientswrtDensity(){
// dRphi_dRho = 0
// dRM_dRho = 0
- // dRrho_dRho =
+ // dRrho_dRho =
// dRv_dRho = 0
// dRblw_dRho =
- // dCl_dRho =
- // dCd_dRho =
+ // dCl_dRho =
+ // dCd_dRho =
}
void CoupledAdjoint::gradientswrtVelocity(){
// dRphi_dV = 0
// dRM_dV = 0
// dRrho_dV = 0
- // dRv_dV =
- // dRblw_dV =
+ // dRv_dV =
+ // dRblw_dV =
- // dCl_dV =
- // dCd_dV =
+ // dCl_dV =
+ // dCd_dV =
}
void CoupledAdjoint::gradientswrtBlowingVelocity(){
- // dRphi_dBlw =
- // dRM_dBlw = 0
- // dRrho_dBlw = 0
- // dRv_dBlw = 0
- // dRblw_dBlw =
-
- // dCl_dBlw =
- // dCd_dBlw =
-}
\ No newline at end of file
+ //dRphi_blw = iadjoint->getRu_Blw();
+
+ // All others are zero.
+}
+
+void CoupledAdjoint::gradientwrtAoA(){
+ std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa();
+ dCl_AoA = dCx_AoA[0]; dCd_AoA = dCx_AoA[1];
+
+ dRphi_AoA = iadjoint->getRu_A();
+
+ // All others are zero.
+}
+
+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
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_dM.size(); i++){
+ for (size_t j = 0; j < _dRblw_dM[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dM[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_M.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_M.prune(0.);
+}
+void CoupledAdjoint::setdRblw_v(std::vector<std::vector<double>> _dRblw_dv){
+ dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_dv.size(), _dRblw_dv[0].size());
+ // Convert std::vector<double> to Eigen::Triplets
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_dv.size(); i++){
+ for (size_t j = 0; j < _dRblw_dv[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_dv[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_v.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_v.prune(0.);
+}
+void CoupledAdjoint::setdRblw_rho(std::vector<std::vector<double>> _dRblw_drho){
+ dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(_dRblw_drho.size(), _dRblw_drho[0].size());
+ // Convert std::vector<double> to Eigen::Triplets
+ std::vector<Eigen::Triplet<double>> triplets;
+ for (size_t i = 0; i < _dRblw_drho.size(); i++){
+ for (size_t j = 0; j < _dRblw_drho[i].size(); j++){
+ triplets.push_back(Eigen::Triplet<double>(i, j, _dRblw_drho[i][j]));
+ }
+ }
+ // Set matrix
+ dRblw_rho.setFromTriplets(triplets.begin(), triplets.end());
+ dRblw_rho.prune(0.);
+}
+
+void CoupledAdjoint::transposeMatrices(){
+ // Transpose all matrices
+ 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();
+
+ 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_phi = dRrho_phi.transpose();
+ dRrho_M = dRrho_M.transpose();
+ dRrho_rho = dRrho_rho.transpose();
+ dRrho_v = dRrho_v.transpose();
+ dRrho_blw = dRrho_blw.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();
+
+ 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();
+}
diff --git a/blast/src/DCoupledAdjoint.h b/blast/src/DCoupledAdjoint.h
index 73aa0b5..93973c7 100644
--- a/blast/src/DCoupledAdjoint.h
+++ b/blast/src/DCoupledAdjoint.h
@@ -19,7 +19,7 @@
#include "blast.h"
#include "wObject.h"
-#include "../../modules/dartflo/dart/src/wAdjoint.h"
+#include "dart.h"
#include <Eigen/Sparse>
#include <iostream>
@@ -39,9 +39,12 @@ public:
std::shared_ptr<dart::Newton> isol; ///< Inviscid Newton solver
std::shared_ptr<blast::Driver> vsol; ///< Viscous solver
std::shared_ptr<tbox::LinearSolver> alinsol; ///< linear solver (adjoint aero)
- double dCldAoA;
- double dCddAoA;
+ 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
private:
size_t sizeBlowing;
@@ -75,12 +78,20 @@ private:
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::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
+
public:
CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint, std::shared_ptr<blast::Driver> vSolver);
- virtual ~CoupledAdjoint () {std::cout << "~CoupledAdjoint()\n";};
+ virtual ~CoupledAdjoint () {std::cout << "~blast::CoupledAdjoint()\n";};
void run();
+ 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);
+
private:
void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor>*>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
@@ -89,6 +100,10 @@ private:
void gradientswrtDensity();
void gradientswrtVelocity();
void gradientswrtBlowingVelocity();
+
+ void gradientwrtAoA();
+
+ void transposeMatrices();
};
} // namespace blast
#endif // DDARTADJOINT
diff --git a/blast/src/DDriver.cpp b/blast/src/DDriver.cpp
index a1fcba0..b937b4f 100644
--- a/blast/src/DDriver.cpp
+++ b/blast/src/DDriver.cpp
@@ -162,7 +162,7 @@ int Driver::run(unsigned int const couplIter)
{
tSolver->setCFL0(CFL0);
tSolver->setitFrozenJac(5);
- tSolver->setinitSln(false);
+ tSolver->setinitSln(true);
// Keyword annoncing iteration safety level.
if (verbose > 1)
diff --git a/blast/tests/dart/adjointVII.py b/blast/tests/dart/adjointVII.py
index 6b2f8bb..a689678 100644
--- a/blast/tests/dart/adjointVII.py
+++ b/blast/tests/dart/adjointVII.py
@@ -105,7 +105,18 @@ def cfgBlast(verb):
'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():
+
+ alpha = 2*np.pi/180
+
# Timer.
tms = fwk.Timers()
tms['total'].start()
@@ -118,15 +129,120 @@ def main():
coupler, iSolverAPI, vSolver = viscUtils.initBlast(icfg, vcfg)
tms['pre'].stop()
- adjointSolver = blast.CoupledAdjoint(iSolverAPI.adjointSolver, vSolver)
+ 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()
#iSolverAPI.run()
- adjointSolver.run()
- print("sucess", adjointSolver)
+ 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
+ coupledAdjointSolver.run()
+ print("sucess", coupledAdjointSolver)
+
+ 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)
+ coupler.run()
+ clV.append(iSolverAPI.getCl())
+ cdV.append(iSolverAPI.getCd())
+ dCl_dalpha_FD_coupled = (clV[-1] - clV[0]) / (2*dalpha)
+ dCd_dalpha_FD_coupled = (cdV[-1] - cdV[0]) / (2*dalpha)
+
+ 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)
+ iSolverAPI.run()
+ clV.append(iSolverAPI.getCl())
+ cdV.append(iSolverAPI.getCd())
+ dCl_dalpha_FD_dart = (clV[-1] - clV[0]) / (2*dalpha)
+ 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(tms)
+
+ 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()
--
GitLab