From fcf75af35840779e440c49caaa56394ad224f0c9 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Sun, 12 Sep 2021 22:54:14 +0200
Subject: [PATCH] Update adjoint (correction + refactoring)
---
dart/interfaces/mphys.py | 28 +-
dart/src/wAdjoint.cpp | 611 +++++++++++++++++++++------------------
dart/src/wAdjoint.h | 41 ++-
dart/tests/adjoint.py | 10 +-
dart/tests/adjoint3.py | 24 +-
5 files changed, 385 insertions(+), 329 deletions(-)
diff --git a/dart/interfaces/mphys.py b/dart/interfaces/mphys.py
index 64591d5..53081bc 100644
--- a/dart/interfaces/mphys.py
+++ b/dart/interfaces/mphys.py
@@ -136,7 +136,7 @@ class DartMorpher(om.ImplicitComponent):
def linearize(self, inputs, outputs, partials):
"""Compute the mesh jacobian matrix
"""
- self.adj.omLinearizeMesh()
+ self.adj.linearizeMesh()
def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
"""Compute the partials derivatives and perform the matrix-vector product
@@ -144,7 +144,7 @@ class DartMorpher(om.ImplicitComponent):
if mode == 'rev':
if 'xv' in d_residuals:
if 'xv' in d_outputs:
- d_outputs['xv'] += self.adj.omComputeJacVecRmeshMesh(d_residuals['xv']) # dX = dRx_dX^T * dRx
+ d_outputs['xv'] += self.adj.computeJacVecMesh(d_residuals['xv']) # dX = dRx_dX^T * dRx
if 'x_aero' in d_inputs:
for i in range(len(self.bnd.nodes)):
for j in range(self.dim):
@@ -156,7 +156,7 @@ class DartMorpher(om.ImplicitComponent):
"""Solve for the partials gradients of the mesh residuals
"""
if mode == 'rev':
- d_residuals['xv'] = self.adj.omSolveMesh(d_outputs['xv']) # dRx = dRx_dX^-T * dX
+ d_residuals['xv'] = self.adj.solveMesh(d_outputs['xv']) # dRx = dRx_dX^-T * dX
elif mode == 'fwd':
raise NotImplementedError('DartMorpher - forward mode not implemented!\n')
@@ -234,9 +234,9 @@ class DartSolver(om.ImplicitComponent):
raise NotImplementedError('DartSolver - cannot compute the residuals!\n')
def linearize(self, inputs, outputs, partials):
- """Compute the flow jacobian matrix
+ """Compute the flow jacobian matrix (and the other partials of the flow residuals)
"""
- self.adj.omLinearizeFlow()
+ self.adj.linearizeFlow()
def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
"""Compute the partials derivatives and perform the matrix-vector product
@@ -244,11 +244,11 @@ class DartSolver(om.ImplicitComponent):
if mode == 'rev':
if 'phi' in d_residuals:
if 'phi' in d_outputs:
- d_outputs['phi'] += self.adj.omComputeJacVecResidualFlow(d_residuals['phi']) # dPhi = dRphi_dPhi^T * dRphi
+ d_outputs['phi'] += self.adj.computeJacVecFlow(d_residuals['phi']) # dU = dRu_dU^T * dRu
if 'aoa' in d_inputs:
- d_inputs['aoa'] += self.adj.omComputeJacVecResidualAoa(d_residuals['phi']) # dAoa = dRphi_dAoa^T * dRphi
+ d_inputs['aoa'] += self.adj.computeJacVecFlowAoa(d_residuals['phi']) # dA = dRu_dA^T * dRu
if 'xv' in d_inputs:
- d_inputs['xv'] += self.adj.omComputeJacVecResidualMesh(d_residuals['phi']) # dX = dRphi_dX^T * dRphi
+ d_inputs['xv'] += self.adj.computeJacVecFlowMesh(d_residuals['phi']) # dX = dRu_dX^T * dRu
elif mode == 'fwd':
raise NotImplementedError('DartSolver - forward mode not implemented!\n')
@@ -256,7 +256,7 @@ class DartSolver(om.ImplicitComponent):
"""Solve for the partials gradients of the flow residuals
"""
if mode == 'rev':
- d_residuals['phi'] = self.adj.omSolveFlow(d_outputs['phi']) # dRphi = dRphi_dPhi^-T * dPhi
+ d_residuals['phi'] = self.adj.solveFlow(d_outputs['phi']) # dRu = dRu_dU^-T * dU
elif mode == 'fwd':
raise NotImplementedError('DartSolver - forward mode not implemented!\n')
@@ -304,9 +304,9 @@ class DartLoads(om.ExplicitComponent):
if mode == 'rev':
if 'f_aero' in d_outputs:
if 'xv' in d_inputs:
- d_inputs['xv'] += self.qinf * np.asarray(self.adj.omComputeJacVecLoadsMesh(d_outputs['f_aero'], self.bnd)) # dX = dL_dX^T * dL
+ d_inputs['xv'] += self.qinf * np.asarray(self.adj.computeJacVecLoadsMesh(d_outputs['f_aero'], self.bnd)) # dX = dL_dX^T * dL
if 'phi' in d_inputs:
- d_inputs['phi'] += self.qinf * np.asarray(self.adj.omComputeJacVecLoadsFlow(d_outputs['f_aero'], self.bnd)) # dPhi = dL_dPhi^T * dL
+ d_inputs['phi'] += self.qinf * np.asarray(self.adj.computeJacVecLoadsFlow(d_outputs['f_aero'], self.bnd)) # dU = dL_dU^T * dL
elif mode == 'fwd':
raise NotImplementedError('DartLoads - forward mode not implemented!\n')
@@ -368,15 +368,15 @@ class DartCoefficients(om.ExplicitComponent):
"""Compute the partials derivatives
"""
# Gradients wrt to angle of attack
- dCfAoa = self.adj.omComputeGradientCoefficientsAoa()
+ dCfAoa = self.adj.computeGradientCoefficientsAoa()
partials['cl', 'aoa'] = dCfAoa[0]
partials['cd', 'aoa'] = dCfAoa[1]
# Gradients wrt to mesh
- dCfMsh = self.adj.omComputeGradientCoefficientsMesh()
+ dCfMsh = self.adj.computeGradientCoefficientsMesh()
partials['cl', 'xv'] = dCfMsh[0]
partials['cd', 'xv'] = dCfMsh[1]
# Gradients wrt to flow variables
- dCfPhi = self.adj.omComputeGradientCoefficientsFlow()
+ dCfPhi = self.adj.computeGradientCoefficientsFlow()
partials['cl', 'phi'] = dCfPhi[0]
partials['cd', 'phi'] = dCfPhi[1]
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index 50e6bdf..103dce0 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -48,12 +48,8 @@
using namespace tbox;
using namespace dart;
-#define ANSI_COLOR_YELLOW "\x1b[1;33m"
-#define ANSI_COLOR_RESET "\x1b[0m"
-
/**
* @brief Initialize the adjoint solver
- * @authors Adrien Crovato
*/
Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph, std::shared_ptr<tbox::LinearSolver> _linsol) : sol(_sol), morph(_morph), linsol(_linsol)
{
@@ -71,9 +67,6 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
dCdMsh.resize(sol->pbl->msh->nodes.size(), Eigen::Vector3d::Zero());
dClAoa = 0;
dCdAoa = 0;
-
- // Setup rows (unknown index)
- this->mapRows();
}
/**
@@ -92,186 +85,329 @@ void Adjoint::run()
<< "Number of threads: " << nthreads << "\n"
<< std::endl;
- // Solve adjoint flow equation
+ // Compute partial gradients of flow residuals and solve flow adjoint
tms["0-AdjFlo"].start();
- this->solve();
+ this->linearizeFlow(); // dRu/dU, dRu/dA, dRu/dX
+ std::vector<std::vector<double>> dCf_U = this->computeGradientCoefficientsFlow(); // {dCl, dCd}/dU
+ lamClFlo = this->solveFlow(dCf_U[0]); // lambdaClU = dRu/dU^-T * dCl/dU
+ lamCdFlo = this->solveFlow(dCf_U[1]); // lambdaCdU = dRu/dU^-T * dCd/dU
tms["0-AdjFlo"].stop();
- // Compute total gradient of loads with respect to angle of attack
+ // Compute gradients with respect to AoA
tms["1-GradAoA"].start();
- this->computeGradientAoa();
+ // partials
+ std::vector<double> dCf_A = this->computeGradientCoefficientsAoa(); // {dCl, dCd}/dA
+ // totals
+ dClAoa = dCf_A[0] - this->computeJacVecFlowAoa(lamClFlo); // dCl/dA - dRu/dA^T * lambdaClU
+ dCdAoa = dCf_A[1] - this->computeJacVecFlowAoa(lamCdFlo); // dCd/dA - dRu/dA^T * lambdaCdU
tms["1-GradAoA"].stop();
- // Compute total gradient of loads with respect to mesh nodes (including deformation)
+ // Compute gradients with respect to mesh coordinates
tms["2-GradMsh"].start();
- this->computeGradientMesh();
+ // compute mesh Jacobian
+ this->linearizeMesh(); // dRx/dX
+ // partials
+ std::vector<std::vector<double>> dCf_X = this->computeGradientCoefficientsMesh(); // {dCl, dCd}/dX
+ std::vector<double> dCl_XU = this->computeJacVecFlowMesh(lamClFlo);
+ std::vector<double> dCd_XU = this->computeJacVecFlowMesh(lamCdFlo);
+ for (size_t i = 0; i < sol->pbl->nDim * sol->pbl->msh->nodes.size(); ++i)
+ {
+ dCf_X[0][i] -= dCl_XU[i]; // dCl/dX - dRu/dX^T * lambdaClU
+ dCf_X[1][i] -= dCd_XU[i]; // dCd/dX - dRu/dX^T * lambdaCdU
+ }
+ // compute mesh adjoint (totals)
+ std::vector<double> lamClMsh = this->solveMesh(dCf_X[0]); // lambdaClX = dRx/dX^T * (dCl/dX - dRu/dX^T * lambdaClU)
+ std::vector<double> lamCdMsh = this->solveMesh(dCf_X[1]); // lambdaCdX = dRx/dX^T * (dCd/dX - dRu/dX^T * lambdaCdU)
+ // store solution
+ for (auto n : sol->pbl->msh->nodes)
+ {
+ for (int m = 0; m < sol->pbl->nDim; m++)
+ {
+ dClMsh[n->row](m) = lamClMsh[sol->pbl->nDim * (n->row) + m];
+ dCdMsh[n->row](m) = lamCdMsh[sol->pbl->nDim * (n->row) + m];
+ }
+ }
tms["2-GradMsh"].stop();
+ // Display
+ // flow adjoint
+ std::cout << std::setw(15) << std::left << "Flow adjoint"
+ << std::setw(15) << std::right << "Res[lambdaCl]"
+ << std::setw(15) << std::right << "Res[lambdaCd]" << std::endl;
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(15) << std::left << ""
+ << std::setw(15) << std::right << log10((dRu_U * Eigen::Map<Eigen::VectorXd>(lamClFlo.data(), lamClFlo.size()) - Eigen::Map<Eigen::VectorXd>(dCf_U[0].data(), dCf_U[0].size())).norm())
+ << std::setw(15) << std::right << log10((dRu_U * Eigen::Map<Eigen::VectorXd>(lamCdFlo.data(), lamCdFlo.size()) - Eigen::Map<Eigen::VectorXd>(dCf_U[1].data(), dCf_U[1].size())).norm()) << std::endl;
+ // aoa gradients
+ std::cout << std::setw(15) << std::left << "AoA gradients"
+ << std::setw(15) << std::right << "dCl_dAoA"
+ << std::setw(15) << std::right << "dCd_dAoA" << std::endl;
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(15) << std::left << ""
+ << std::setw(15) << std::right << dClAoa
+ << std::setw(15) << std::right << dCdAoa << std::endl;
+ // mesh adjoint
+ std::cout << std::setw(15) << std::left << "Mesh adjoint"
+ << std::setw(15) << std::right << "Res[lambdaCl]"
+ << std::setw(15) << std::right << "Res[lambdaCd]" << std::endl;
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(15) << std::left << ""
+ << std::setw(15) << std::right << log10((dRx_X * Eigen::Map<Eigen::VectorXd>(lamClMsh.data(), lamClMsh.size()) - Eigen::Map<Eigen::VectorXd>(dCf_X[0].data(), dCf_X[0].size())).norm())
+ << std::setw(15) << std::right << log10((dRx_X * Eigen::Map<Eigen::VectorXd>(lamCdMsh.data(), lamCdMsh.size()) - Eigen::Map<Eigen::VectorXd>(dCf_X[1].data(), dCf_X[1].size())).norm()) << std::endl;
+ // mesh gradients
+ std::cout << std::setw(15) << std::left << "Mesh gradients"
+ << std::setw(15) << std::right << "Norm[dCl_dX]"
+ << std::setw(15) << std::right << "Norm[dCd_dX]" << std::endl;
+ std::cout << std::fixed << std::setprecision(5);
+ std::cout << std::setw(15) << std::left << ""
+ << std::setw(15) << std::right << Eigen::Map<Eigen::VectorXd>(lamClMsh.data(), lamClMsh.size()).norm()
+ << std::setw(15) << std::right << Eigen::Map<Eigen::VectorXd>(lamCdMsh.data(), lamCdMsh.size()).norm() << std::endl;
+
// Display timers
- if (verbose > 0)
+ if (verbose > 1)
std::cout << "Adjoint solver CPU" << std::endl
<< tms;
std::cout << std::endl;
}
/**
- * @brief Solve the adjoint steady transonic Full Potential Equation
+ * @brief Compute the mesh Jacobian
+ *
+ * J = dRx/dX
*/
-void Adjoint::solve()
+void Adjoint::linearizeMesh()
{
- // Build partial gradients of residual wrt to potential
- tms["00-AdjFloJac"].start();
- Eigen::SparseMatrix<double, Eigen::RowMajor> dR(sol->pbl->msh->nodes.size(), sol->pbl->msh->nodes.size());
- Eigen::SparseMatrix<double, Eigen::RowMajor> dR_(sol->pbl->msh->nodes.size(), sol->pbl->msh->nodes.size());
- sol->buildJac(dR_);
- dR = dR_.transpose();
- tms["00-AdjFloJac"].stop();
- tms["01-AdjFloRes"].start();
- // Build partial gradients of loads wrt potential
- Eigen::VectorXd dCl(sol->pbl->msh->nodes.size()), dCd(sol->pbl->msh->nodes.size());
- this->buildGradientLoadsFlow(dCl, dCd);
- tms["01-AdjFloRes"].stop();
-
- // Solve linear adjoint equations: dR {lambdaCl, lambdaCd} = {dCl, dCd}
- tms["02-AdjFloSol"].start();
- Eigen::Map<Eigen::VectorXd> lamCl_(lamClFlo.data(), lamClFlo.size()), lamCd_(lamCdFlo.data(), lamCdFlo.size());
- linsol->compute(dR, Eigen::Map<Eigen::VectorXd>(dCl.data(), dCl.size()), lamCl_);
- linsol->compute(dR, Eigen::Map<Eigen::VectorXd>(dCd.data(), dCd.size()), lamCd_);
- tms["02-AdjFloSol"].stop();
+ dRx_X = Eigen::SparseMatrix<double, Eigen::RowMajor>(sol->pbl->nDim * sol->pbl->msh->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
+ morph->build(dRx_X); // J = dRx/dX
+}
- // Display
- std::cout << std::setw(15) << std::left << "Flow adjoint"
- << std::setw(15) << std::right << "Res[lambdaCl]"
- << std::setw(15) << std::right << "Res[lambdaCd]" << std::endl;
- std::cout << std::fixed << std::setprecision(5);
- std::cout << std::setw(15) << std::left << ""
- << std::setw(15) << std::right << log10((dR * lamCl_ - dCl).norm())
- << std::setw(15) << std::right << log10((dR * lamCd_ - dCd).norm()) << std::endl;
+/**
+ * @brief Solve for the mesh residuals given a gradient (wrt to mesh coordinates)
+ *
+ * solve for dRx: dRx_dX^T * dRx = dX
+ */
+std::vector<double> Adjoint::solveMesh(std::vector<double> const &dX)
+{
+ // Create maps
+ Eigen::VectorXd dR(sol->pbl->nDim * sol->pbl->msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> dR_(dR.data(), dR.size());
+ Eigen::VectorXd dX_ = Eigen::Map<const Eigen::VectorXd>(dX.data(), dX.size());
+ // Solve
+ linsol->compute(dRx_X.transpose(), Eigen::Map<Eigen::VectorXd>(dX_.data(), dX_.size()), dR_);
+ return std::vector<double>(dR.data(), dR.data() + dR.size());
}
/**
- * @brief Compute total gradients of the loads with respect to angle of attack
- * @todo one gradient for each surface node
+ * @brief Compute the matrix-vector product of the mesh Jacobian with the mesh residual
+ *
+ * dX = dRx/dX^T * dRx
+ * @todo not tested
*/
-void Adjoint::computeGradientAoa()
+std::vector<double> Adjoint::computeJacVecMesh(std::vector<double> const &dR)
{
- // Compute partial gradients with respect to AoA
- tms["10-GradAoARes"].start();
- Eigen::VectorXd dR = Eigen::VectorXd::Zero(sol->pbl->msh->nodes.size());
- this->buildGradientResidualAoa(dR);
- tms["10-GradAoARes"].stop();
- tms["11-GradAoALoad"].start();
- double dL = 0, dD = 0;
- this->buildGradientLoadsAoa(dL, dD);
- tms["11-GradAoALoad"].stop();
- // Compute total gradients
- dClAoa = (dL - Eigen::Map<Eigen::VectorXd>(lamClFlo.data(), lamClFlo.size()).transpose() * dR);
- dCdAoa = (dD - Eigen::Map<Eigen::VectorXd>(lamCdFlo.data(), lamCdFlo.size()).transpose() * dR);
+ Eigen::VectorXd JdR = dRx_X.transpose() * Eigen::Map<const Eigen::VectorXd>(dR.data(), dR.size());
+ return std::vector<double>(JdR.data(), JdR.data() + JdR.size());
+}
- // Display
- std::cout << std::setw(15) << std::left << "AoA gradients"
- << std::setw(15) << std::right << "dCl_dAoA"
- << std::setw(15) << std::right << "dCd_dAoA" << std::endl;
- std::cout << std::fixed << std::setprecision(5);
- std::cout << std::setw(15) << std::left << ""
- << std::setw(15) << std::right << dClAoa
- << std::setw(15) << std::right << dCdAoa << std::endl;
+/**
+ * @brief Compute the partial gradients of the flow residuals
+ *
+ * dRu_U = dRu / dU
+ * dRu_A = dRu / dA
+ * dRu_X = dRu / dX
+ */
+void Adjoint::linearizeFlow()
+{
+ // Partials wrt flow (Jacobian)
+ dRu_U = Eigen::SparseMatrix<double, Eigen::RowMajor>(sol->pbl->msh->nodes.size(), sol->pbl->msh->nodes.size());
+ sol->buildJac(dRu_U);
+ // Partials wrt AoA
+ dRu_A = Eigen::VectorXd::Zero(sol->pbl->msh->nodes.size());
+ this->buildGradientResidualAoa(dRu_A);
+ // Partials wrt mesh
+ dRu_X = Eigen::SparseMatrix<double, Eigen::RowMajor>(sol->pbl->msh->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
+ this->buildGradientResidualMesh(dRu_X);
}
/**
- * @brief Compute total gradients of the loads with respect to mesh nodes
+ * @brief Solve for the flow residuals given a gradient (wrt to flow variables)
+ *
+ * solve for dRu: dRu/dU^T * dRu = dU
*/
-void Adjoint::computeGradientMesh()
+std::vector<double> Adjoint::solveFlow(std::vector<double> const &dU)
{
- // Build partial gradients of flow residual wrt mesh
- tms["20-GradMshRes"].start();
- Eigen::SparseMatrix<double, Eigen::RowMajor> dR(sol->pbl->msh->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
- this->buildGradientResidualMesh(dR);
- tms["20-GradMshRes"].stop();
- // Build partial gradients of loads wrt mesh
- tms["21-GradMshLoad"].start();
- Eigen::RowVectorXd dCl = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size()), dCd = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size());
- this->buildGradientLoadsMesh(dCl, dCd);
- tms["21-GradMshLoad"].stop();
-
- // Build partial gradients of mesh deformation residual wrt mesh
- tms["22-GradMshAdj"].start();
- Eigen::SparseMatrix<double, Eigen::RowMajor> Kmesh(sol->pbl->nDim * sol->pbl->msh->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
- Eigen::SparseMatrix<double, Eigen::RowMajor> Kmesh_(sol->pbl->nDim * sol->pbl->msh->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
- morph->build(Kmesh_);
- Kmesh = Kmesh_.transpose();
- // Build RHS of mesh deformation adjoint equations: {dCl, dCd} = {dCl, dCd}_x - {lambdaCl, lambdaCd}_flow * dR_x
- Eigen::VectorXd dClXyz = (dCl - Eigen::Map<Eigen::VectorXd>(lamClFlo.data(), lamClFlo.size()).transpose() * dR).transpose();
- Eigen::VectorXd dCdXyz = (dCd - Eigen::Map<Eigen::VectorXd>(lamCdFlo.data(), lamCdFlo.size()).transpose() * dR).transpose();
- tms["22-GradMshAdj"].stop();
- // Solve mesh deformation adjoint equations: K {lambdaCl, lambdaCd}_mesh = {dCl, dCd}
- tms["23-GradMshSol"].start();
- Eigen::VectorXd lamClMsh = Eigen::VectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size());
- Eigen::VectorXd lamCdMsh = Eigen::VectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size());
- Eigen::Map<Eigen::VectorXd> lamClMsh_(lamClMsh.data(), lamClMsh.size()), lamCdMsh_(lamCdMsh.data(), lamCdMsh.size());
- linsol->compute(Kmesh, Eigen::Map<Eigen::VectorXd>(dClXyz.data(), dClXyz.size()), lamClMsh_);
- linsol->compute(Kmesh, Eigen::Map<Eigen::VectorXd>(dCdXyz.data(), dCdXyz.size()), lamCdMsh_);
- tms["23-GradMshSol"].stop();
- // Store solution
- for (auto n : sol->pbl->msh->nodes)
+ // Create maps
+ Eigen::VectorXd dR(sol->pbl->msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> dR_(dR.data(), dR.size());
+ Eigen::VectorXd dU_ = Eigen::Map<const Eigen::VectorXd>(dU.data(), dU.size());
+ // Solve
+ linsol->compute(dRu_U.transpose(), Eigen::Map<Eigen::VectorXd>(dU_.data(), dU_.size()), dR_);
+ return std::vector<double>(dR.data(), dR.data() + dR.size());
+}
+
+/**
+ * @brief Compute the matrix-vector product of the flow Jacobian with the flow residuals
+ *
+ * dU = dRu/dU^T * dRu
+ * @todo not tested
+ */
+std::vector<double> Adjoint::computeJacVecFlow(std::vector<double> const &dR)
+{
+ Eigen::VectorXd JdR = dRu_U.transpose() * Eigen::Map<const Eigen::VectorXd>(dR.data(), dR.size());
+ return std::vector<double>(JdR.data(), JdR.data() + JdR.size());
+}
+
+/**
+ * @brief Compute the matrix-vector product of the gradients of the flow residuals (wrt angle of attack) with the flow residuals
+ *
+ * dA = dRu/dA^T * dRu
+ */
+double Adjoint::computeJacVecFlowAoa(std::vector<double> const &dR)
+{
+ return dRu_A.transpose() * Eigen::Map<const Eigen::VectorXd>(dR.data(), dR.size());
+}
+
+/**
+ * @brief Compute the matrix-vector product of the gradients of the flow residuals (wrt mesh coordinates) with the flow residuals
+ *
+ * dX = dRu/dX^T * dRu
+ */
+std::vector<double> Adjoint::computeJacVecFlowMesh(std::vector<double> const &dR)
+{
+ Eigen::VectorXd dX = dRu_X.transpose() * Eigen::Map<const Eigen::VectorXd>(dR.data(), dR.size());
+ return std::vector<double>(dX.data(), dX.data() + dX.size());
+}
+
+/**
+ * @brief Compute the matrix-vector product of the gradients of the loads (wrt flow variables) with the loads, on a given boundary
+ *
+ * dU = dL/dU^T * dU
+ * @todo only for FSI, not tested by run()
+ * @todo linearize outside (check om)
+ */
+std::vector<double> Adjoint::computeJacVecLoadsFlow(std::vector<double> const &dL, Boundary const &bnd)
+{
+ // Build gradients
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dL_U = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd.nodes.size(), sol->pbl->msh->nodes.size());
+ this->buildGradientLoadsFlow(dL_U, bnd);
+ // Compute product
+ Eigen::VectorXd dU = dL_U.transpose() * Eigen::Map<const Eigen::VectorXd>(dL.data(), dL.size());
+ return std::vector<double>(dU.data(), dU.data() + dU.size());
+}
+
+/**
+ * @brief Compute the matrix-vector product of the gradients of the loads (wrt the mesh coordinates) with the loads, on a given boundary
+ *
+ * dX = dL/dX^T * dL
+ * @todo only for FSI, not tested by run()
+ * @todo linearize outside (check om)
+ */
+std::vector<double> Adjoint::computeJacVecLoadsMesh(std::vector<double> const &dL, Boundary const &bnd)
+{
+ // Build gradients
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dL_X = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd.nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
+ this->buildGradientLoadsMesh(dL_X, bnd);
+ // Compute product
+ Eigen::VectorXd dX = dL_X.transpose() * Eigen::Map<const Eigen::VectorXd>(dL.data(), dL.size()); // dMsh = dLoads_dMsh^T * dLoads
+ return std::vector<double>(dX.data(), dX.data() + dX.size());
+}
+
+/**
+ * @brief Compute the gradients of the load coefficients wrt flow variables
+ *
+ * dCl/dU, dCd/dU
+ */
+std::vector<std::vector<double>> Adjoint::computeGradientCoefficientsFlow()
+{
+ // Compute gradients of lift and drag
+ Eigen::RowVectorXd dCl = Eigen::RowVectorXd::Zero(sol->pbl->msh->nodes.size()), dCd = Eigen::RowVectorXd::Zero(sol->pbl->msh->nodes.size());
+ for (auto bnd : sol->pbl->bnds)
{
- for (int m = 0; m < sol->pbl->nDim; m++)
+ // gradients of loads
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dL = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd->nodes.size(), sol->pbl->msh->nodes.size());
+ this->buildGradientLoadsFlow(dL, *bnd);
+ // project to lift and drag directions
+ Eigen::RowVector3d dirL = sol->pbl->dirL.eval().transpose(), dirD = sol->pbl->dirD.eval().transpose();
+ for (size_t i = 0; i < bnd->nodes.size(); ++i)
{
- dClMsh[n->row](m) = lamClMsh(sol->pbl->nDim * (n->row) + m);
- dCdMsh[n->row](m) = lamCdMsh(sol->pbl->nDim * (n->row) + m);
+ dCl += dirL * dL.middleRows(i * 3, 3);
+ dCd += dirD * dL.middleRows(i * 3, 3);
}
}
+ return std::vector<std::vector<double>>{std::vector<double>(dCl.data(), dCl.data() + dCl.size()), std::vector<double>(dCd.data(), dCd.data() + dCd.size())};
+}
- // Display
- std::cout << std::setw(15) << std::left << "Mesh adjoint"
- << std::setw(15) << std::right << "Res[lambdaCl]"
- << std::setw(15) << std::right << "Res[lambdaCd]" << std::endl;
- std::cout << std::fixed << std::setprecision(5);
- std::cout << std::setw(15) << std::left << ""
- << std::setw(15) << std::right << log10((Kmesh * lamClMsh - dClXyz).norm())
- << std::setw(15) << std::right << log10((Kmesh * lamCdMsh - dCdXyz).norm()) << std::endl;
- std::cout << std::setw(15) << std::left << "Mesh gradients"
- << std::setw(15) << std::right << "Norm[dCl_dX]"
- << std::setw(15) << std::right << "Norm[dCd_dX]" << std::endl;
- std::cout << std::setw(15) << std::left << ""
- << std::setw(15) << std::right << lamClMsh.norm()
- << std::setw(15) << std::right << lamCdMsh.norm() << std::endl;
+/**
+ * @brief Compute the gradients of the load coefficients wrt angle of attack
+ *
+ * dCl/dA, dCd/dA
+ */
+std::vector<double> Adjoint::computeGradientCoefficientsAoa()
+{
+ double dCl, dCd;
+ this->buildGradientLoadsAoa(dCl, dCd);
+ return std::vector<double>{dCl, dCd};
+}
+
+/**
+ * @brief Compute the gradients of the load coefficients wrt mesh coordinates
+ *
+ * dCl/dX, dCd/dX
+ */
+std::vector<std::vector<double>> Adjoint::computeGradientCoefficientsMesh()
+{
+ // Compute the gradients of the lift and drag
+ Eigen::RowVectorXd dCl = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size()), dCd = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size());
+ for (auto bnd : sol->pbl->bnds)
+ {
+ // gradients of loads
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dL = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
+ this->buildGradientLoadsMesh(dL, *bnd);
+ // project to lift and drag directions
+ Eigen::RowVector3d dirL = sol->pbl->dirL.eval().transpose(), dirD = sol->pbl->dirD.eval().transpose();
+ for (size_t i = 0; i < bnd->nodes.size(); ++i)
+ {
+ dCl += dirL * dL.middleRows(i * 3, 3);
+ dCd += dirD * dL.middleRows(i * 3, 3);
+ }
+ }
+ return std::vector<std::vector<double>>{std::vector<double>(dCl.data(), dCl.data() + dCl.size()), std::vector<double>(dCd.data(), dCd.data() + dCd.size())};
}
/**
- * @brief Build the gradients of the lift and drag with respect to the flow variables
- * @todo one gradient for each surface node
+ * @brief Build the gradients of the loads with respect to the flow variables one a given boundary
*/
-void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD)
+void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd)
{
// Multithread
+ tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
tbb::spin_mutex mutex;
- // Reset and compute gradient of lift and drag
- dL.setZero();
- dD.setZero();
- for (auto sur : sol->pbl->bnds)
- {
- tbb::parallel_for_each(sur->groups[0]->tag->elems.begin(), sur->groups[0]->tag->elems.end(), [&](Element *e) {
- // Associated volume element
- Element *eV = sur->svMap.at(e);
- // Build
- Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, *sol->pbl->medium->cP);
- // Assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t i = 0; i < e->nodes.size(); ++i)
+ // Build gradients
+ std::deque<Eigen::Triplet<double>> T;
+ tbb::parallel_for_each(bnd.groups[0]->tag->elems.begin(), bnd.groups[0]->tag->elems.end(), [&](Element *e) {
+ // Associated volume element
+ Element *eV = bnd.svMap.at(e);
+ // Build
+ Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, *sol->pbl->medium->cP);
+ // Assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t i = 0; i < e->nodes.size(); ++i)
+ {
+ int rowi = bnd.nMap.at(e->nodes[i]);
+ for (int m = 0; m < 3; ++m)
{
- Eigen::RowVectorXd AeL = sol->pbl->dirL.eval().transpose() * Ae.block(i * 3, 0, 3, eV->nodes.size());
- Eigen::RowVectorXd AeD = sol->pbl->dirD.eval().transpose() * Ae.block(i * 3, 0, 3, eV->nodes.size());
for (size_t j = 0; j < eV->nodes.size(); ++j)
{
- Node *nodj = eV->nodes[j];
- dL(nodj->row) += AeL(j);
- dD(nodj->row) += AeD(j);
+ int rowj = eV->nodes[j]->row;
+ T.push_back(Eigen::Triplet<double>(3 * rowi + m, rowj, Ae(3 * i + m, j)));
}
}
- });
- }
+ }
+ });
+ dL.setFromTriplets(T.begin(), T.end());
+ dL.prune(0.);
+ dL.makeCompressed();
}
/**
@@ -280,6 +416,7 @@ void Adjoint::buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD)
void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
{
// Multithread
+ tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
tbb::spin_mutex mutex;
// Reset gradient
@@ -305,9 +442,8 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
/**
* @brief Compute gradients of the lift and drag with respect to angle of attack
- * @todo one gradient for each surface node
*/
-void Adjoint::buildGradientLoadsAoa(double &dL, double &dD)
+void Adjoint::buildGradientLoadsAoa(double &dCl, double &dCd)
{
// Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
Eigen::Vector3d Cf(0, 0, 0);
@@ -315,16 +451,17 @@ void Adjoint::buildGradientLoadsAoa(double &dL, double &dD)
for (size_t i = 0; i < bnd->nodes.size(); ++i)
Cf += bnd->nLoads[i];
// Rotate to gradient of flow direction
- dD = Cf.dot(sol->pbl->dirD.evalGrad());
- dL = Cf.dot(sol->pbl->dirL.evalGrad());
+ dCl = Cf.dot(sol->pbl->dirL.evalGrad());
+ dCd = Cf.dot(sol->pbl->dirD.evalGrad());
}
/**
- * @brief Compute gradients of the residual with respect to mesh nodes
+ * @brief Compute gradients of the residual with respect to mesh coordinates
*/
void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
{
// Multithread
+ tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
tbb::spin_mutex mutex;
// List of triplets
@@ -348,14 +485,9 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
size_t rowi = e->nodes[i]->row;
for (size_t j = 0; j < e->nodes.size(); ++j)
{
+ int rowj = e->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[e->nodes[j]->row].size(); ++jj)
- {
- int rowj = sol->pbl->nDim * (arows[e->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(sol->rows[rowi], rowj, Ae1(i, sol->pbl->nDim * j + n)));
- }
- }
+ T.push_back(Eigen::Triplet<double>(sol->rows[rowi], sol->pbl->nDim * rowj + n, Ae1(i, sol->pbl->nDim * j + n)));
}
}
// Assembly (supersonic)
@@ -367,14 +499,9 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
size_t rowi = e->nodes[i]->row;
for (size_t j = 0; j < eU->nodes.size(); ++j)
{
+ int rowj = eU->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[eU->nodes[j]->row].size(); ++jj)
- {
- int rowj = sol->pbl->nDim * (arows[eU->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(sol->rows[rowi], rowj, Ae2(i, sol->pbl->nDim * j + n)));
- }
- }
+ T.push_back(Eigen::Triplet<double>(sol->rows[rowi], sol->pbl->nDim * rowj + n, Ae2(i, sol->pbl->nDim * j + n)));
}
}
}
@@ -395,36 +522,21 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
int rowi = we->wkN[i].second->row;
for (size_t j = 0; j < we->nColUp; ++j)
{
+ int rowj = we->volUpE->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[we->volUpE->nodes[j]->row].size(); ++jj)
- {
- int rowj = sol->pbl->nDim * (arows[we->volUpE->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(rowi, rowj, Aeu(i, sol->pbl->nDim * j + n)));
- }
- }
+ T.push_back(Eigen::Triplet<double>(rowi, sol->pbl->nDim * rowj + n, Aeu(i, sol->pbl->nDim * j + n)));
}
for (size_t j = 0; j < we->nColLw; ++j)
{
+ int rowj = we->volLwE->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[we->volLwE->nodes[j]->row].size(); ++jj)
- {
- int rowj = sol->pbl->nDim * (arows[we->volLwE->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(rowi, rowj, -Ael(i, sol->pbl->nDim * j + n)));
- }
- }
+ T.push_back(Eigen::Triplet<double>(rowi, sol->pbl->nDim * rowj + n, -Ael(i, sol->pbl->nDim * j + n)));
}
for (size_t j = 0; j < we->surUpE->nodes.size(); ++j)
{
+ int rowj = we->surUpE->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[we->surUpE->nodes[j]->row].size(); ++jj)
- {
- int rowj = sol->pbl->nDim * (arows[we->surUpE->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(rowi, rowj, Aes(i, sol->pbl->nDim * j + n)));
- }
- }
+ T.push_back(Eigen::Triplet<double>(rowi, sol->pbl->nDim * rowj + n, Aes(i, sol->pbl->nDim * j + n)));
}
}
});
@@ -443,25 +555,15 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
int rowi = ke->teN[i].second->row;
for (size_t j = 0; j < ke->nCol; ++j)
{
+ int rowj = ke->volE->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[ke->volE->nodes[j]->row].size(); ++jj)
- {
- int rowj = sol->pbl->nDim * (arows[ke->volE->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(rowi, rowj, Ae(i, sol->pbl->nDim * j + n)));
- }
- }
+ T.push_back(Eigen::Triplet<double>(rowi, sol->pbl->nDim * rowj + n, Ae(i, sol->pbl->nDim * j + n)));
}
for (size_t j = 0; j < ke->surE->nodes.size(); ++j)
{
+ int rowj = ke->surE->nodes[j]->row;
for (int n = 0; n < sol->pbl->nDim; ++n)
- {
- for (size_t jj = 0; jj < arows[ke->surE->nodes[j]->row].size(); ++jj)
- {
- int rowj = sol->pbl->nDim * (arows[ke->surE->nodes[j]->row][jj]) + n;
- T.push_back(Eigen::Triplet<double>(rowi, rowj, Aes(i, sol->pbl->nDim * j + n)));
- }
- }
+ T.push_back(Eigen::Triplet<double>(rowi, sol->pbl->nDim * rowj + n, Aes(i, sol->pbl->nDim * j + n)));
}
}
});
@@ -473,110 +575,49 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
}
/**
- * @brief Compute gradients of the lift and drag with respect to mesh nodes
- * @todo one gradient for each surface node
+ * @brief Compute gradients of the loads with respect to mesh coordinates, on a given boundary
*/
-void Adjoint::buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd &dD)
+void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd)
{
// Multithread
+ tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
tbb::spin_mutex mutex;
- // Reset and compute gradient of lift and drag
- dL.setZero();
- dD.setZero();
- for (auto bnd : sol->pbl->bnds)
- {
- tbb::parallel_for_each(bnd->groups[0]->tag->elems.begin(), bnd->groups[0]->tag->elems.end(), [&](Element *e) {
- // Volume element
- Element *eV = bnd->svMap.at(e);
- // Build
- Eigen::MatrixXd Aes, Aev;
- std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, *sol->pbl->medium->cP, sol->pbl->nDim);
- // Assembly
- tbb::spin_mutex::scoped_lock lock(mutex);
- for (size_t i = 0; i < e->nodes.size(); ++i)
+ // Build gradients
+ std::deque<Eigen::Triplet<double>> T;
+ tbb::parallel_for_each(bnd.groups[0]->tag->elems.begin(), bnd.groups[0]->tag->elems.end(), [&](Element *e) {
+ // Volume element
+ Element *eV = bnd.svMap.at(e);
+ // Build
+ Eigen::MatrixXd Aes, Aev;
+ std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, *sol->pbl->medium->cP, sol->pbl->nDim);
+ // Assembly
+ tbb::spin_mutex::scoped_lock lock(mutex);
+ for (size_t i = 0; i < e->nodes.size(); ++i)
+ {
+ int rowi = bnd.nMap.at(e->nodes[i]);
+ for (int m = 0; m < 3; ++m)
{
- // rotate to flow direction
- Eigen::RowVectorXd AesL = sol->pbl->dirL.eval().transpose() * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
- Eigen::RowVectorXd AesD = sol->pbl->dirD.eval().transpose() * Aes.block(i * 3, 0, 3, sol->pbl->nDim * e->nodes.size());
- Eigen::RowVectorXd AevL = sol->pbl->dirL.eval().transpose() * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
- Eigen::RowVectorXd AevD = sol->pbl->dirD.eval().transpose() * Aev.block(i * 3, 0, 3, sol->pbl->nDim * eV->nodes.size());
- for (int m = 0; m < sol->pbl->nDim; ++m)
+ // surface
+ for (size_t j = 0; j < e->nodes.size(); ++j)
{
- // surface
- for (size_t j = 0; j < e->nodes.size(); ++j)
- {
- size_t rowj = e->nodes[j]->row;
- for (size_t jj = 0; jj < arows[rowj].size(); ++jj)
- {
- dL.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AesL.col(j * sol->pbl->nDim + m);
- dD.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AesD.col(j * sol->pbl->nDim + m);
- }
- }
- // volume
- for (size_t j = 0; j < eV->nodes.size(); ++j)
- {
- size_t rowj = eV->nodes[j]->row;
- for (size_t jj = 0; jj < arows[rowj].size(); ++jj)
- {
- dL.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AevL.col(j * sol->pbl->nDim + m);
- dD.col((arows[rowj][jj]) * sol->pbl->nDim + m) += AevD.col(j * sol->pbl->nDim + m);
- }
- }
+ int rowj = e->nodes[j]->row;
+ for (int n = 0; n < sol->pbl->nDim; ++n)
+ T.push_back(Eigen::Triplet<double>(3 * rowi + m, sol->pbl->nDim * rowj + n, Aes(3 * i + m, sol->pbl->nDim * j + n)));
}
- }
- });
- }
-}
-
-/**
- * @brief Create the rows vector for the adjoint mesh
- * @todo duplicated from wake
- */
-void Adjoint::mapRows()
-{
- // Collect each row associated to each node
- arows.resize(sol->pbl->msh->nodes.size());
- for (size_t i = 0; i < sol->pbl->msh->nodes.size(); ++i)
- arows[i].push_back(sol->pbl->msh->nodes[i]->row);
-
- // Add duplicated row on the wake
- for (auto wake : sol->pbl->wBCs)
- {
- // Collect nodes on the wake
- std::vector<Node *> nodesUp, nodesLw;
- for (auto e : wake->groups[0]->tag->elems)
- for (auto n : e->nodes)
- nodesUp.push_back(n);
- for (auto e : wake->groups[1]->tag->elems)
- for (auto n : e->nodes)
- nodesLw.push_back(n);
- std::sort(nodesUp.begin(), nodesUp.end());
- auto it = std::unique(nodesUp.begin(), nodesUp.end());
- nodesUp.resize(std::distance(nodesUp.begin(), it));
- std::sort(nodesLw.begin(), nodesLw.end());
- it = std::unique(nodesLw.begin(), nodesLw.end());
- nodesLw.resize(std::distance(nodesLw.begin(), it));
- // Create map between upper and lower sides, but only retain different nodes
- std::map<Node *, Node *> nodMap;
- for (auto nUp : nodesUp)
- {
- for (auto nLw : nodesLw)
- {
- if ((nUp->pos - nLw->pos).norm() <= 1e-14 && nUp != nLw)
+ // volume
+ for (size_t j = 0; j < eV->nodes.size(); ++j)
{
- nodMap[nUp] = nLw;
- break;
+ size_t rowj = eV->nodes[j]->row;
+ for (int n = 0; n < sol->pbl->nDim; ++n)
+ T.push_back(Eigen::Triplet<double>(3 * rowi + m, sol->pbl->nDim * rowj + n, Aev(3 * i + m, sol->pbl->nDim * j + n)));
}
}
}
- // Create rows
- for (auto p : nodMap)
- {
- arows[p.first->row].push_back(p.second->row);
- arows[p.second->row].push_back(p.first->row);
- }
- }
+ });
+ dL.setFromTriplets(T.begin(), T.end());
+ dL.prune(0.);
+ dL.makeCompressed();
}
/**
diff --git a/dart/src/wAdjoint.h b/dart/src/wAdjoint.h
index ecac9d7..a378133 100644
--- a/dart/src/wAdjoint.h
+++ b/dart/src/wAdjoint.h
@@ -32,6 +32,7 @@ namespace dart
/**
* @brief Adjoint solver class
* @authors Adrien Crovato
+ * @todo add interface for Eigen to avoid using Eigen::Map
*/
class DART_API Adjoint : public fwk::wSharedObject
{
@@ -51,8 +52,11 @@ public:
double dCdAoa; ///< drag gradient with respect to angle of attack
private:
- fwk::Timers tms; ///< internal timers
- std::vector<std::vector<int>> arows; ///< rows (unknown index) with duplicated dof on wake
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_X; ///< mesh Jacobian
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRu_U; ///< flow Jacobian
+ Eigen::VectorXd dRu_A; ///< partial gradients of flow residual wrt aoa
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRu_X; ///< partial gradients of flow residual wrt mesh
+ fwk::Timers tms; ///< internal timers
public:
Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph, std::shared_ptr<tbox::LinearSolver> _linsol);
@@ -61,26 +65,37 @@ public:
void run();
void save(tbox::MshExport *mshWriter, int n = 0);
+ // Partial gradients of the mesh residuals
+ void linearizeMesh();
+ std::vector<double> solveMesh(std::vector<double> const &dMesh);
+ std::vector<double> computeJacVecMesh(std::vector<double> const &dR);
+ // Partial gradients of the flow residuals
+ void linearizeFlow();
+ std::vector<double> solveFlow(std::vector<double> const &dFlow);
+ std::vector<double> computeJacVecFlow(std::vector<double> const &dR);
+ double computeJacVecFlowAoa(std::vector<double> const &dR);
+ std::vector<double> computeJacVecFlowMesh(std::vector<double> const &dR);
+ // Partial gradients of the loads
+ std::vector<double> computeJacVecLoadsFlow(std::vector<double> const &dLoads, Boundary const &bnd);
+ std::vector<double> computeJacVecLoadsMesh(std::vector<double> const &dLoads, Boundary const &bnd);
+ // Partial gradients of the coefficients
+ std::vector<std::vector<double>> computeGradientCoefficientsFlow();
+ std::vector<double> computeGradientCoefficientsAoa();
+ std::vector<std::vector<double>> computeGradientCoefficientsMesh();
+
#ifndef SWIG
virtual void write(std::ostream &out) const override;
#endif
private:
- // Adjoint solution and total gradients
- void solve();
- void computeGradientAoa();
- void computeGradientMesh();
// Partial gradients wrt with respect to flow variables
- void buildGradientLoadsFlow(Eigen::VectorXd &dL, Eigen::VectorXd &dD);
+ void buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd);
// Partial gradients with respect to angle of attack
void buildGradientResidualAoa(Eigen::VectorXd &dR);
- void buildGradientLoadsAoa(double &dL, double &dD);
- // Partial gradients with respect to mesh
+ void buildGradientLoadsAoa(double &dCl, double &dCd);
+ // Partial gradients with respect to mesh coordinates
void buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR);
- void buildGradientLoadsMesh(Eigen::RowVectorXd &dL, Eigen::RowVectorXd &dD);
-
- // @todo
- void mapRows();
+ void buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd);
};
} // namespace dart
diff --git a/dart/tests/adjoint.py b/dart/tests/adjoint.py
index 254fcf7..7151ad4 100644
--- a/dart/tests/adjoint.py
+++ b/dart/tests/adjoint.py
@@ -88,7 +88,7 @@ def main():
dClAoA = 0
dCdAoA = 0
for n in bnd.nodes:
- dx = drot.dot(np.array([n.pos[0] - 1, n.pos[1]]))
+ dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1]]))
dClAoA += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1]]).dot(dx)
dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1]]).dot(dx)
@@ -114,15 +114,15 @@ def main():
if M_inf == 0. and alpha == 0*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 6.9, 1e-2))
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.0, 1e-3))
- tests.add(CTest('dCl_dMsh', dClX, 61.302, 1e-3))
- tests.add(CTest('dCd_dMsh', dCdX, 0.0749, 1e-3))
+ tests.add(CTest('dCl_dMsh', dClX, 55.439, 1e-3))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.482, 1e-3))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, 6.9, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.0, 1e-3))
elif M_inf == 0.7 and alpha == 2*np.pi/180:
tests.add(CTest('dCl_dAoA', adjoint.dClAoa, 11.8, 1e-2)) # FD 11.835 (step = 1e-5)
tests.add(CTest('dCd_dAoA', adjoint.dCdAoa, 0.127, 1e-2)) # 0.12692 (step = 1e-5)
- tests.add(CTest('dCl_dMsh', dClX, 99.89, 1e-3))
- tests.add(CTest('dCd_dMsh', dCdX, 0.747, 1e-3))
+ tests.add(CTest('dCl_dMsh', dClX, 82.67, 1e-3))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.908, 1e-3))
tests.add(CTest('dCl_dAoA (msh)', dClAoA, 11.8, 1e-2))
tests.add(CTest('dCd_dAoA (msh)', dCdAoA, 0.127, 1e-2))
else:
diff --git a/dart/tests/adjoint3.py b/dart/tests/adjoint3.py
index 2aa430e..2ccb806 100644
--- a/dart/tests/adjoint3.py
+++ b/dart/tests/adjoint3.py
@@ -89,7 +89,7 @@ def main():
dClAoA = 0
dCdAoA = 0
for n in bnd.nodes:
- dx = drot.dot(np.array([n.pos[0] - 1, n.pos[1], n.pos[2]]))
+ dx = drot.dot(np.array([n.pos[0] - 0.25, n.pos[1], n.pos[2]]))
dClAoA += np.array([adjoint.dClMsh[n.row][0], adjoint.dClMsh[n.row][1], adjoint.dClMsh[n.row][2]]).dot(dx)
dCdAoA += np.array([adjoint.dCdMsh[n.row][0], adjoint.dCdMsh[n.row][1], adjoint.dCdMsh[n.row][2]]).dot(dx)
@@ -110,17 +110,17 @@ def main():
if M_inf == 0.3 and alpha == 3*np.pi/180 and span == 1:
tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 2.5, 5e-2)) # FD 2.526 (step 1e-5)
tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.146, 5e-2)) # FD 0.1460 (step 1e-5)
- tests.add(CTest('dCl_dMsh', dClX, 3.778, 5e-2))
- tests.add(CTest('dCd_dMsh', dCdX, 0.192, 5e-2))
- tests.add(CTest('dCl/dAoA (msh)', dClAoA, 2.7, 5e-2))
- tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.151, 5e-2))
- elif M_inf == 0.7 and alpha == 5*np.pi/180 and span == 2: # mesh size must be 0.01
- tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 5.0, 1e-2))
- tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.427, 1e-2))
- tests.add(CTest('dCl_dMsh', dClX, 4.773, 1e-3))
- tests.add(CTest('dCd_dMsh', dCdX, 0.327, 1e-3))
- tests.add(CTest('dCl/dAoA (msh)', dClAoA, 5.2, 1e-2))
- tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.441, 1e-2))
+ tests.add(CTest('dCl_dMsh', dClX, 2.016, 5e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.133, 5e-2))
+ tests.add(CTest('dCl/dAoA (msh)', dClAoA, 2.5, 5e-2))
+ tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.146, 5e-2))
+ elif M_inf == 0.7 and alpha == 5*np.pi/180 and span == 2:
+ tests.add(CTest('dCl/dAoA', adjoint.dClAoa, 4.9, 5e-2))
+ tests.add(CTest('dCd/dAoA', adjoint.dCdAoa, 0.454, 5e-2))
+ tests.add(CTest('dCl_dMsh', dClX, 2.571, 5e-2))
+ tests.add(CTest('dCd_dMsh', dCdX, 0.204, 5e-2))
+ tests.add(CTest('dCl/dAoA (msh)', dClAoA, 4.9, 5e-2))
+ tests.add(CTest('dCd/dAoA (msh)', dCdAoA, 0.454, 5e-2))
else:
raise Exception('Test not defined for this flow')
tests.run()
--
GitLab