Skip to content
Snippets Groups Projects
Commit fcf75af3 authored by Adrien Crovato's avatar Adrien Crovato
Browse files

Update adjoint (correction + refactoring)

parent 2f3ee55a
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
dart/interfaces/mphys.py
+ 14
− 14
View file @ fcf75af3
......@@ -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]
......
dart/src/wAdjoint.cpp
+ 326
− 285
View file @ fcf75af3
This diff is collapsed.
dart/src/wAdjoint.h
+ 28
− 13
View file @ fcf75af3
......@@ -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
......
dart/tests/adjoint.py
+ 5
− 5
View file @ fcf75af3
......@@ -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:
......
dart/tests/adjoint3.py
+ 12
− 12
View file @ fcf75af3
......@@ -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()
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment