From 9116fa17743f0af9e385bfa6a7fe3dbea13ccd7f Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Tue, 25 Mar 2025 02:50:28 +0100
Subject: [PATCH] (feat) Added auto-diff adjoint solver for boundary layer
---
blast/_src/blastw.h | 1 +
blast/_src/blastw.i | 5 +
blast/interfaces/blSolversInterface.py | 29 +-
blast/src/blAdjoint.cpp | 496 ++++++++++
blast/src/blAdjoint.h | 92 ++
blast/src/blBody.h | 2 +-
blast/src/blCoupledAdjoint.cpp | 1168 ++++++++----------------
blast/src/blCoupledAdjoint.h | 289 +++---
blast/src/blDriver.h | 14 +-
blast/src/blNode.cpp | 2 +-
blast/src/blNode.h | 7 +-
blast/src/blast.h | 5 +
blast/tests/dart/t_adjoint_2D.py | 9 +-
13 files changed, 1142 insertions(+), 977 deletions(-)
create mode 100644 blast/src/blAdjoint.cpp
create mode 100644 blast/src/blAdjoint.h
diff --git a/blast/_src/blastw.h b/blast/_src/blastw.h
index 75d5d85..27cefda 100644
--- a/blast/_src/blastw.h
+++ b/blast/_src/blastw.h
@@ -21,5 +21,6 @@
#include "blNode.h"
// Solvers
+#include "blAdjoint.h"
#include "blCoupledAdjoint.h"
#include "blDriver.h"
diff --git a/blast/_src/blastw.i b/blast/_src/blastw.i
index 58432be..98c3c28 100644
--- a/blast/_src/blastw.i
+++ b/blast/_src/blastw.i
@@ -101,6 +101,11 @@ typedef double bldouble; // needed so that SWIG does not convert bldouble to dou
%immutable blast::Driver::verbose;
%include "blDriver.h"
+#ifdef USE_CODI
+%shared_ptr(blast::blAdjoint)
+%include "blAdjoint.h"
+#endif // USE_CODI
+
%shared_ptr(blast::CoupledAdjoint);
%immutable blast::CoupledAdjoint::tdCl_AoA;
%immutable blast::CoupledAdjoint::tdCd_AoA;
diff --git a/blast/interfaces/blSolversInterface.py b/blast/interfaces/blSolversInterface.py
index 66cec41..559a132 100644
--- a/blast/interfaces/blSolversInterface.py
+++ b/blast/interfaces/blSolversInterface.py
@@ -363,7 +363,8 @@ class SolversInterface:
for inod, nod in enumerate(pystruct.getNodes(reg.getName())):
pos = blast.blVector3d(nod[xIdx], nod[yIdx], nod[zIdx])
row = int(pystruct.getNodeConnectList(reg.getName())[inod]) if self.cfg['interpolator'] == 'Matching' else inod
- nodes.push_back(blast.blNode(inod, pos, row))
+ adjrow = int(pystruct.getNodeRows(reg.getName())[inod]) if self.cfg['interpolator'] == 'Matching' else inod
+ nodes.push_back(blast.blNode(inod, pos, row, adjrow))
# Set nodes
reg.setMesh(nodes)
# Set elements connectivity
@@ -698,7 +699,11 @@ class SolversInterface:
lower_side = lower_side[unique_lower]
chord = pts[ite, 0] - pts[ile, 0]
- x_interp = (chord/2) * (np.cos(zeta)+1) + pts[ile,0]
+ if self.cfg['cosine_spacing']:
+ xx = (chord/2) * (np.cos(zeta) + 1)
+ else:
+ xx = chord * np.linspace(1, 0, len(zeta))
+ x_interp = xx + pts[ile,0]
upper_interp = interp1d(upper_side[:,0], upper_side[:,2], kind='linear')(x_interp)
lower_interp = interp1d(lower_side[:,0], lower_side[:,2], kind='linear')(x_interp)
upper_final = np.column_stack((x_interp, np.full_like(x_interp, ys[isec]), upper_interp))
@@ -725,7 +730,11 @@ class SolversInterface:
ptsw = ptsw[ptsw[:,0].argsort()]
chordw = ptsw[-1, 0] - ptsw[0, 0]
- x_interp = np.flip((chordw/2) * (np.cos(zeta)+1) + ptsw[0,0])
+ if self.cfg['cosine_spacing']:
+ xx = (chordw/2) * (np.cos(zeta) + 1)
+ else:
+ xx = chordw * np.linspace(1, 0, len(zeta))
+ x_interp = np.flip( xx + ptsw[0,0])
interp = interp1d(ptsw[:,0], ptsw[:,2], kind='linear')(x_interp)
wake_final = np.column_stack((x_interp, np.full_like(x_interp, ys[isec]), interp))
@@ -836,6 +845,8 @@ class SolversInterface:
cfg['sections'] = {'airfoil': [0]}
else:
raise RuntimeError('Expected "sections" in cfg')
+ if 'cosine_spacing' not in cfg:
+ cfg['cosine_spacing'] = True
if len(cfg['sections']) != self.getnBodies():
raise RuntimeError(f'Expected the length of "sections" to be {self.getnBodies()}, got {len(cfg["sections"])}')
return cfg
@@ -1043,13 +1054,13 @@ class SolversInterface:
elems = region["elems"][:]
pyreg.setMesh(nodes, elems)
pyreg.setBlowingVelocity(region["blowing"])
+
# Coupling parameters are defined on upper/lower sides and wake
-
- for iReg in range(len(self.deltaStarExt[ibody][isec])):
- blregion = section[f"blInteraction{iReg}"]
- self.deltaStarExt[ibody][isec][iReg] = blregion["deltaStarExt"][:]
- self.vtExt[ibody][isec][iReg] = blregion["vtExt"][:]
- self.xxExt[ibody][isec][iReg] = blregion["xxExt"][:]
+ for iReg in range(len(self.deltaStarExt[ibody][isec])):
+ blregion = section[f"blInteraction{iReg}"]
+ self.deltaStarExt[ibody][isec][iReg] = blregion["deltaStarExt"][:]
+ self.vtExt[ibody][isec][iReg] = blregion["vtExt"][:]
+ self.xxExt[ibody][isec][iReg] = blregion["xxExt"][:]
# Abstract methods for the interface
def run()->int:
diff --git a/blast/src/blAdjoint.cpp b/blast/src/blAdjoint.cpp
new file mode 100644
index 0000000..64a592c
--- /dev/null
+++ b/blast/src/blAdjoint.cpp
@@ -0,0 +1,496 @@
+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "blast.h"
+#include "blAdjoint.h"
+#include "blDriver.h"
+#include "blBoundaryLayer.h"
+
+#include <iostream>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+using namespace blast;
+
+blAdjoint::blAdjoint(std::shared_ptr<blast::Driver> &vsolver, int _ndim) : vsol(vsolver), ndim(_ndim)
+{
+ // Fill adjoint map
+ adjointMatrixMap["dRdStar_M"] = &dRdStar_M;
+ adjointMatrixMap["dRdStar_vt"] = &dRdStar_vt;
+ adjointMatrixMap["dRdStar_rho"] = &dRdStar_rho;
+ adjointMatrixMap["dRdStar_dStar"] = &dRdStar_dStar;
+ adjointMatrixMap["dRdStar_blw"] = &dRdStar_blw;
+ adjointMatrixMap["dRdStar_x"] = &dRdStar_x;
+ adjointMatrixMap["dRblw_M"] = &dRblw_M;
+ adjointMatrixMap["dRblw_vt"] = &dRblw_vt;
+ adjointMatrixMap["dRblw_rho"] = &dRblw_rho;
+ adjointMatrixMap["dRblw_dStar"] = &dRblw_dStar;
+ adjointMatrixMap["dRblw_blw"] = &dRblw_blw;
+ adjointMatrixMap["dRblw_x"] = &dRblw_x;
+ adjointVectorMap["dCdf_M"] = &dCdf_M;
+ adjointVectorMap["dCdf_rho"] = &dCdf_rho;
+ adjointVectorMap["dCdf_vt"] = &dCdf_vt;
+ adjointVectorMap["dCdf_dStar"] = &dCdf_dStar;
+ adjointVectorMap["dCdf_blw"] = &dCdf_dStar;
+ adjointVectorMap["dCdf_x"] = &dCdf_x;
+}
+
+void blAdjoint::reset()
+{
+ // Compute total number of nodes and elements
+ size_t nNodes = 0;
+ size_t nElems = 0;
+ for (auto &body : vsol->bodies)
+ for (auto §ion : body->sections)
+ for (auto ®ion : section->regions)
+ {
+ nNodes += region->getnNodes();
+ nElems += region->getnElms();
+ }
+
+ dRdStar_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRdStar_vt = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRdStar_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRdStar_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRdStar_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nElems);
+ dRdStar_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes * ndim);
+
+ dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nElems, nNodes);
+ dRblw_vt = Eigen::SparseMatrix<double, Eigen::RowMajor>(nElems, nNodes);
+ dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nElems, nNodes);
+ dRblw_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nElems, nNodes);
+ dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nElems, nElems);
+ dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nElems, nNodes * ndim);
+
+ dRcf_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRcf_vt = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRcf_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRcf_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes);
+ dRcf_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNodes, nNodes * ndim);
+
+ dCdf_M = Eigen::VectorXd::Zero(nNodes);
+ dCdf_vt = Eigen::VectorXd::Zero(nNodes);
+ dCdf_rho = Eigen::VectorXd::Zero(nNodes);
+ dCdf_dStar = Eigen::VectorXd::Zero(nNodes);
+ dCdf_x = Eigen::VectorXd::Zero(nNodes * ndim);
+}
+
+void blAdjoint::run()
+{
+ this->reset();
+ this->solve();
+ this->evaluate(adjointMatrixMap, adjointVectorMap);
+}
+
+/**
+ * @brief Solve the forward boundary layer problem
+ * and record the operations
+ */
+int blAdjoint::solve()
+{
+ std::cout << "Running boundary layer solver..." << std::endl;
+ bltape &adjtape = bldouble::getTape();
+ adjtape.reset();
+ adjtape.setActive();
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ®ion : section->regions)
+ {
+ for (size_t inode = 0; inode < region->getnNodes(); inode++)
+ {
+ for (size_t idim = 0; idim < ndim; ++idim)
+ adjtape.registerInput(region->nodes[inode]->pos[idim]);
+ adjtape.registerInput(region->Me[inode]);
+ adjtape.registerInput(region->vt[inode]);
+ adjtape.registerInput(region->deltaStarExt[inode]);
+ }
+ }
+ section->computeLocalCoordinates();
+ }
+ }
+ int excode = vsol->run();
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ®ion : section->regions)
+ {
+ for (size_t inode = 0; inode < region->getnNodes(); inode++)
+ {
+ adjtape.registerOutput(region->deltaStar[inode]);
+ adjtape.registerOutput(region->cf[inode]);
+ }
+ }
+ section->computeLocalCoordinates();
+ }
+ }
+ adjtape.registerOutput(vsol->Cdf);
+ adjtape.setPassive();
+ return excode;
+}
+
+void blAdjoint::evaluate(std::map<std::string, Eigen::SparseMatrix<double, Eigen::RowMajor> *> &matrixMap,
+ std::map<std::string, Eigen::VectorXd *> &vectorMap)
+{
+
+ // Check if required fields are in the map
+ for (auto &item : adjointMatrixMap)
+ {
+ if (matrixMap.find(item.first) == matrixMap.end())
+ {
+ throw std::runtime_error(item.first + " not found");
+ }
+ }
+ for (auto &item : adjointVectorMap)
+ {
+ if (vectorMap.find(item.first) == vectorMap.end())
+ {
+ throw std::runtime_error(item.first + " not found");
+ }
+ }
+
+ // Gradient wrt blowing velocity
+ gradientBlowing(matrixMap);
+
+ bltape &adjtape = bldouble::getTape();
+ adjtape.clearAdjoints();
+
+ std::deque<Eigen::Triplet<double>> Trip_dstar_M;
+ std::deque<Eigen::Triplet<double>> Trip_dstar_rho;
+ std::deque<Eigen::Triplet<double>> Trip_dstar_vt;
+ std::deque<Eigen::Triplet<double>> Trip_dstar_dstar;
+ std::deque<Eigen::Triplet<double>> Trip_dstar_x;
+
+ // std::deque<Eigen::Triplet<double>> Trip_cf_M;
+ // std::deque<Eigen::Triplet<double>> Trip_cf_rho;
+ // std::deque<Eigen::Triplet<double>> Trip_cf_vt;
+ // std::deque<Eigen::Triplet<double>> Trip_cf_dstar;
+ // std::deque<Eigen::Triplet<double>> Trip_cf_x;
+
+ size_t jj = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ®ion : section->regions)
+ {
+ for (size_t inode = 0; inode < region->getnNodes(); inode++)
+ {
+ region->deltaStar[inode].gradient() = 1.0;
+ adjtape.evaluate();
+ size_t kk = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ®ion : section->regions)
+ {
+ for (size_t inode = 0; inode < region->getnNodes(); inode++)
+ {
+ int rowj = region->nodes[inode]->adjrow;
+ double gradM = region->Me[inode].getGradient();
+ double gradrho = region->rhoe[inode].getGradient();
+ double gradvt = region->vt[inode].getGradient();
+ double graddstar = region->deltaStarExt[inode].getGradient();
+ if (kk == jj)
+ graddstar -= 1.0; // because Rdstar = deltaStar - f(deltaStarExt)
+ Trip_dstar_M.push_back(Eigen::Triplet<double>(jj, kk, -gradM));
+ Trip_dstar_rho.push_back(Eigen::Triplet<double>(jj, kk, -gradrho));
+ Trip_dstar_vt.push_back(Eigen::Triplet<double>(jj, kk, -gradvt));
+ Trip_dstar_dstar.push_back(Eigen::Triplet<double>(jj, kk, -graddstar));
+ for (int idim = 0; idim < ndim; idim++)
+ {
+ double gradx = region->nodes[inode]->pos[idim].getGradient();
+ Trip_dstar_x.push_back(Eigen::Triplet<double>(jj, rowj * ndim + idim, -gradx));
+ }
+ kk++;
+ }
+ }
+ }
+ }
+ adjtape.clearAdjoints();
+ jj++;
+ }
+ }
+ }
+ }
+
+ // Friction coefficient
+ // jj = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ®ion : section->regions)
+ // {
+ // for (size_t inode = 0; inode < region->getnNodes(); inode++)
+ // {
+ // region->cf[inode].gradient() = 1.0;
+ // adjtape.evaluate();
+ // size_t kk = 0;
+ // for (auto &body : vsol->bodies)
+ // {
+ // for (auto §ion : body->sections)
+ // {
+ // for (auto ®ion : section->regions)
+ // {
+ // for (size_t inode = 0; inode < region->getnNodes(); inode++)
+ // {
+ // int rowj = region->nodes[inode]->adjrow;
+ // double gradM = region->Me[inode].getGradient();
+ // double gradrho = region->rhoe[inode].getGradient();
+ // double gradvt = region->vt[inode].getGradient();
+ // double graddstar = region->deltaStarExt[inode].getGradient();
+ // Trip_cf_M.push_back(Eigen::Triplet<double>(jj, kk, -gradM));
+ // Trip_cf_rho.push_back(Eigen::Triplet<double>(jj, kk, -gradrho));
+ // Trip_cf_vt.push_back(Eigen::Triplet<double>(jj, kk, -gradvt));
+ // Trip_cf_dstar.push_back(Eigen::Triplet<double>(jj, kk, -graddstar));
+ // for (int idim = 0; idim < ndim; idim++)
+ // {
+ // double gradx = region->nodes[inode]->pos[idim].getGradient();
+ // Trip_cf_x.push_back(Eigen::Triplet<double>(jj, rowj * ndim + idim, -gradx));
+ // }
+ // kk++;
+ // }
+ // }
+ // }
+ // }
+ // adjtape.clearAdjoints();
+ // jj++;
+ // }
+ // }
+ // }
+ // }
+
+ // Grad Cdf
+ vsol->Cdf.gradient() = 1.0;
+ adjtape.evaluate();
+ size_t kk = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ®ion : section->regions)
+ {
+ for (size_t inode = 0; inode < region->getnNodes(); inode++)
+ {
+ int rowj = region->nodes[inode]->adjrow;
+ double gradM = region->Me[inode].getGradient();
+ double gradrho = region->rhoe[inode].getGradient();
+ double gradvt = region->vt[inode].getGradient();
+ double graddstar = region->deltaStarExt[inode].getGradient();
+ vectorMap["dCdf_M"]->coeffRef(kk) = gradM;
+ vectorMap["dCdf_rho"]->coeffRef(kk) = gradrho;
+ vectorMap["dCdf_vt"]->coeffRef(kk) = gradvt;
+ vectorMap["dCdf_dStar"]->coeffRef(kk) = graddstar;
+ for (int idim = 0; idim < ndim; idim++)
+ {
+ double gradx = region->nodes[inode]->pos[idim].getGradient();
+ vectorMap["dCdf_x"]->coeffRef(rowj * ndim + idim) = gradx;
+ }
+ kk++;
+ }
+ }
+ }
+ }
+
+ matrixMap["dRdStar_M"]->setFromTriplets(Trip_dstar_M.begin(), Trip_dstar_M.end());
+ matrixMap["dRdStar_rho"]->setFromTriplets(Trip_dstar_rho.begin(), Trip_dstar_rho.end());
+ matrixMap["dRdStar_vt"]->setFromTriplets(Trip_dstar_vt.begin(), Trip_dstar_vt.end());
+ matrixMap["dRdStar_dStar"]->setFromTriplets(Trip_dstar_dstar.begin(), Trip_dstar_dstar.end());
+ matrixMap["dRdStar_x"]->setFromTriplets(Trip_dstar_x.begin(), Trip_dstar_x.end());
+ adjtape.printStatistics();
+
+ // dRcf_M.setFromTriplets(Trip_cf_M.begin(), Trip_cf_M.end());
+ // dRcf_rho.setFromTriplets(Trip_cf_rho.begin(), Trip_cf_rho.end());
+ // dRcf_vt.setFromTriplets(Trip_cf_vt.begin(), Trip_cf_vt.end());
+ // dRcf_dStar.setFromTriplets(Trip_cf_dstar.begin(), Trip_cf_dstar.end());
+ // dRcf_x.setFromTriplets(Trip_cf_x.begin(), Trip_cf_x.end());
+}
+
+void blAdjoint::gradientBlowing(std::map<std::string, Eigen::SparseMatrix<double, Eigen::RowMajor> *> &map)
+{
+ // Compute total number of nodes and elements
+ size_t nNodes = 0;
+ size_t nElems = 0;
+ for (auto &body : vsol->bodies)
+ for (auto §ion : body->sections)
+ for (auto ®ion : section->regions)
+ {
+ nNodes += region->getnNodes();
+ nElems += region->getnElms();
+ }
+
+ //**************************************************************************//
+ // dRblw_rho //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ std::deque<Eigen::Triplet<double>> Trip_blw_rho;
+ int offSetElms = 0;
+ int offSetNodes = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ std::vector<std::vector<bldouble>> grad = reg->evalGradwrtRho();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ Trip_blw_rho.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
+ -grad[i][j].getValue()));
+ offSetElms += reg->getnElms();
+ offSetNodes += reg->getnNodes();
+ }
+ }
+ }
+ map["dRblw_rho"]->setFromTriplets(Trip_blw_rho.begin(), Trip_blw_rho.end());
+
+ //**************************************************************************//
+ // dRblw_vt //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ std::deque<Eigen::Triplet<double>> Trip_blw_vt;
+ offSetElms = 0;
+ offSetNodes = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ std::vector<std::vector<bldouble>> grad = reg->evalGradwrtVt();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ Trip_blw_vt.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
+ -grad[i][j].getValue()));
+ offSetElms += reg->getnElms();
+ offSetNodes += reg->getnNodes();
+ }
+ }
+ }
+ blSparseMatrix dRblw_vt(nElems, nNodes);
+ map["dRblw_vt"]->setFromTriplets(Trip_blw_vt.begin(), Trip_blw_vt.end());
+
+ //**************************************************************************//
+ // dRblw_dStar //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ std::deque<Eigen::Triplet<double>> Trip_blw_dstar;
+ offSetElms = 0;
+ offSetNodes = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ std::vector<std::vector<bldouble>> grad = reg->evalGradwrtDeltaStar();
+ for (size_t i = 0; i < grad.size(); ++i)
+ for (size_t j = 0; j < grad[i].size(); ++j)
+ Trip_blw_dstar.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
+ -grad[i][j].getValue()));
+ offSetElms += reg->getnElms();
+ offSetNodes += reg->getnNodes();
+ }
+ }
+ }
+ map["dRblw_dStar"]->setFromTriplets(Trip_blw_dstar.begin(), Trip_blw_dstar.end());
+
+ //**************************************************************************//
+ // dRblw_blw //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ map["dRblw_blw"]->setIdentity();
+
+ std::deque<Eigen::Triplet<double>> T_blw_x;
+ offSetElms = 0;
+ int iReg = -1;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ std::vector<std::vector<bldouble>> gradX = reg->evalGradwrtX();
+ std::vector<std::vector<bldouble>> gradY = reg->evalGradwrtY();
+ for (size_t i = 0; i < gradX.size(); ++i)
+ for (size_t j = 0; j < gradX[i].size(); ++j)
+ {
+ int rowj = reg->nodes[j]->adjrow;
+ T_blw_x.push_back(Eigen::Triplet<double>(
+ i + offSetElms, rowj * ndim + 0, -gradX[i][j].getValue()));
+ T_blw_x.push_back(Eigen::Triplet<double>(
+ i + offSetElms, rowj * ndim + 1, -gradY[i][j].getValue()));
+ }
+ offSetElms += reg->getnElms();
+ }
+ }
+ }
+ map["dRblw_x"]->setFromTriplets(T_blw_x.begin(), T_blw_x.end());
+}
+
+/**
+ * @brief Write the matrices to files
+ * @param matrix Eigen matrix to write
+ * @param filename Name of the file to write
+ * @return void
+ *
+ * The matrices are written in a file with the following format:
+ * - The first line contains the column headers
+ * - The following lines contain the matrix data
+ * - The matrix data is written with a precision of 12 digits
+ */
+void blAdjoint::writeMatrixToFile(const blMatrixXd &matrix,
+ const std::string &filename)
+{
+ std::ofstream file(filename);
+ if (file.is_open())
+ {
+ // Write the column headers
+ for (int j = 0; j < matrix.cols(); ++j)
+ {
+ file << std::setw(15) << j;
+ if (j != matrix.cols() - 1)
+ {
+ file << " ";
+ }
+ }
+ file << "\n";
+
+ // Write the matrix data
+ for (int i = 0; i < matrix.rows(); ++i)
+ {
+ for (int j = 0; j < matrix.cols(); ++j)
+ {
+ file << std::fixed << std::setprecision(12) << std::setw(15)
+ << matrix(i, j);
+ if (j != matrix.cols() - 1)
+ {
+ file << " ";
+ }
+ }
+ file << "\n";
+ }
+ file.close();
+ }
+}
diff --git a/blast/src/blAdjoint.h b/blast/src/blAdjoint.h
new file mode 100644
index 0000000..8958e98
--- /dev/null
+++ b/blast/src/blAdjoint.h
@@ -0,0 +1,92 @@
+/*
+ * Copyright 2025 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef BLADJOINT_H
+#define BLADJOINT_H
+
+#include "blast.h"
+#include "wObject.h"
+#include "wTimers.h"
+#include <memory>
+#include <Eigen/Sparse>
+
+#include <iostream>
+#include <memory>
+
+namespace blast
+{
+
+/**
+ * @brief Adjoint of the boundary layer solver
+ * @author Paul Dechamps
+ */
+class BLAST_API blAdjoint : public fwk::wSharedObject
+{
+ friend class CoupledAdjoint;
+
+private:
+ std::shared_ptr<blast::Driver> &vsol; ///< Viscous solver
+
+ std::map<std::string, Eigen::SparseMatrix<double, Eigen::RowMajor> *> adjointMatrixMap;
+ std::map<std::string, Eigen::VectorXd *> adjointVectorMap;
+ size_t ndim; ///< Number of dimensions
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_M; ///< Partial delta* wrt Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt; ///< Partial delta* wrt velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_rho; ///< Partial delta* wrt density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar; ///< Partial delta* with respect to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_blw; ///< Partial delta* wrt blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_x; ///< Partial delta* wrt mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRcf_M; ///< Partial skin friction coefficient wrt Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRcf_rho; ///< Partial skin friction coefficient wrt density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRcf_vt; ///< Partial skin friction coefficient wrt velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRcf_dStar; ///< Partial skin friction coefficient with respect to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRcf_blw; ///< Partial skin friction coefficient with respect to blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRcf_x; ///< Partial skin friction coefficient wrt mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_M; ///< Partial blowing velocity wrt Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt; ///< Partial blowing velocity wrt velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_rho; ///< Partial blowing velocity wrt density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_dStar; //< Partial blowing velocity with respect to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_blw; ///< Partial blowing velocity wrt blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_x; ///< Partial blowing velocity wrt mesh coordinates
+
+ // Objective functions
+ Eigen::VectorXd dCdf_M; ///< Partial drag coefficient wrt Mach number
+ Eigen::VectorXd dCdf_rho; ///< Partial drag coefficient wrt density
+ Eigen::VectorXd dCdf_vt; ///< Partial drag coefficient wrt velocity
+ Eigen::VectorXd dCdf_x; ///< Partial drag coefficient wrt mesh coordinates
+ Eigen::VectorXd dCdf_dStar; ///< Partial drag coefficient with respect to delta*
+ Eigen::VectorXd dCdf_blw; ///< Partial drag coefficient with respect to blowing velocity
+
+public:
+ blAdjoint(std::shared_ptr<blast::Driver> &vsolver, int _ndim);
+ virtual ~blAdjoint() { std::cout << "~blast::Adjoint()\n"; }
+
+ int solve();
+ void evaluate(std::map<std::string, Eigen::SparseMatrix<double, Eigen::RowMajor> *> &matrixMap,
+ std::map<std::string, Eigen::VectorXd *> &vectorMap);
+ void run();
+ void reset();
+ void writeMatrixToFile(const blMatrixXd &matrix,
+ const std::string &filename);
+
+private:
+ void gradientBlowing(std::map<std::string, Eigen::SparseMatrix<double, Eigen::RowMajor> *> &matrixMap);
+};
+} // namespace blast
+#endif // BLADJOINT_H
diff --git a/blast/src/blBody.h b/blast/src/blBody.h
index c029eac..07bbe45 100644
--- a/blast/src/blBody.h
+++ b/blast/src/blBody.h
@@ -19,7 +19,7 @@ class BLAST_API Body : public fwk::wSharedObject
{
private:
std::string name; ///< Name of the body.
- bldouble span; ///< Body span. Used for drag computation.
+ bldouble span; ///< Body span. Used for drag computation.
public:
bldouble Cdt; ///< Total drag coefficient of the body
diff --git a/blast/src/blCoupledAdjoint.cpp b/blast/src/blCoupledAdjoint.cpp
index 1e8867f..82a39b6 100644
--- a/blast/src/blCoupledAdjoint.cpp
+++ b/blast/src/blCoupledAdjoint.cpp
@@ -17,6 +17,7 @@
#include "blCoupledAdjoint.h"
#include "blBody.h"
#include "blBoundaryLayer.h"
+#include "blAdjoint.h"
#include "blDriver.h"
#include "wAdjoint.h"
#include "wBlowing.h"
@@ -46,10 +47,11 @@
using namespace blast;
CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint,
- std::shared_ptr<blast::Driver> vSolver)
- : iadjoint(_iadjoint), vsol(vSolver)
+ std::shared_ptr<blast::blAdjoint> _vadjoint)
+ : iadjoint(_iadjoint), vadjoint(_vadjoint)
{
isol = iadjoint->sol;
+ vsol = vadjoint->vsol;
// Linear solvers (if GMRES, use the same, but with a thighter tolerance)
if (isol->linsol->type() == tbox::LSolType::GMRES_ILUT)
{
@@ -67,64 +69,64 @@ CoupledAdjoint::CoupledAdjoint(std::shared_ptr<dart::Adjoint> _iadjoint,
void CoupledAdjoint::reset()
{
- size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
+ size_t ndim = isol->pbl->nDim; // Number of dimensions of the problem
size_t nNd = isol->pbl->msh->nodes.size(); // Number of domain nodes
- size_t nNs = 0;
+ size_t nnodes = 0;
for (auto bBC : isol->pbl->bBCs)
- nNs += bBC->nodes.size(); // Number of surface nodes
- nNs += 1; // Duplicate stagnation point
- size_t nEs = 0;
+ nnodes += bBC->nodes.size(); // Number of surface nodes
+ nnodes += 1; // Duplicate stagnation point
+ size_t nelms = 0;
for (auto bBC : isol->pbl->bBCs)
- nEs += bBC->uE.size(); // Number of blowing elements
+ nelms += bBC->uE.size(); // Number of blowing elements
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);
- dRdStar_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(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);
- dRdStar_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(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);
- dRdStar_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(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);
- dRdStar_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNs);
- dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim * nNs);
-
- dRphi_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nNs);
- dRM_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
- dRrho_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
- dRv_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nNs);
- dRdStar_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nNs);
- dRblw_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, 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);
- dRdStar_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nEs);
- dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nEs);
-
- dRphi_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nDim * nNd);
- dRM_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd);
- dRrho_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd);
- dRv_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNs, nDim * nNd);
- dRdStar_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNs, nDim * nNd);
- dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nEs, nDim * nNd);
+ dRM_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nNd);
+ dRrho_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nNd);
+ dRv_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nnodes, nNd);
+ dRdStar_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nNd);
+ dRblw_phi = Eigen::SparseMatrix<double, Eigen::RowMajor>(nelms, nNd);
+
+ dRphi_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nnodes);
+ dRM_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRrho_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRv_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nnodes, nnodes);
+ dRdStar_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRblw_M = Eigen::SparseMatrix<double, Eigen::RowMajor>(nelms, nnodes);
+
+ dRphi_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nnodes);
+ dRM_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRrho_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRv_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nnodes, nnodes);
+ dRdStar_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRblw_rho = Eigen::SparseMatrix<double, Eigen::RowMajor>(nelms, nnodes);
+
+ dRphi_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, ndim * nnodes);
+ dRM_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, ndim * nnodes);
+ dRrho_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, ndim * nnodes);
+ dRv_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nnodes, ndim * nnodes);
+ dRdStar_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, ndim * nnodes);
+ dRblw_v = Eigen::SparseMatrix<double, Eigen::RowMajor>(nelms, ndim * nnodes);
+
+ dRphi_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nnodes);
+ dRM_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRrho_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRv_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nnodes, nnodes);
+ dRdStar_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nnodes);
+ dRblw_dStar = Eigen::SparseMatrix<double, Eigen::RowMajor>(nelms, nnodes);
+
+ dRphi_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, nelms);
+ dRM_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nelms);
+ dRrho_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nelms);
+ dRv_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nnodes, nelms);
+ dRdStar_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, nelms);
+ dRblw_blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(nelms, nelms);
+
+ dRphi_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nNd, ndim * nNd);
+ dRM_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, ndim * nNd);
+ dRrho_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, ndim * nNd);
+ dRv_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nnodes, ndim * nNd);
+ dRdStar_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nnodes, ndim * nNd);
+ dRblw_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nelms, ndim * nNd);
// Gradient wrt objective function
dRphi_AoA = Eigen::VectorXd::Zero(nNd);
@@ -133,9 +135,10 @@ void CoupledAdjoint::reset()
// Cl = f(phi) & Cd = Cdp + Cdf = f(phi, M, v)
dCl_phi = Eigen::VectorXd::Zero(nNd);
dCd_phi = Eigen::VectorXd::Zero(nNd);
- dCd_M = Eigen::VectorXd::Zero(nNs);
- dCd_v = Eigen::VectorXd::Zero(nDim * nNs);
- dCd_dStar = Eigen::VectorXd::Zero(nNs);
+ dCd_M = Eigen::VectorXd::Zero(nnodes);
+ dCd_rho = Eigen::VectorXd::Zero(nnodes);
+ dCd_v = Eigen::VectorXd::Zero(ndim * nnodes);
+ dCd_dStar = Eigen::VectorXd::Zero(nnodes);
// Angle of attack
dCl_AoA = 0.0;
@@ -144,19 +147,16 @@ void CoupledAdjoint::reset()
tdCd_AoA = 0.0;
// Mesh
- dRx_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(nDim * nNd, nDim * nNd);
- dCl_x = Eigen::VectorXd::Zero(nDim * nNd);
- dCdp_x = Eigen::VectorXd::Zero(nDim * nNd);
- dCdf_x = Eigen::VectorXd::Zero(nDim * nNd);
+ dRx_x = Eigen::SparseMatrix<double, Eigen::RowMajor>(ndim * nNd, ndim * nNd);
+ dCl_x = Eigen::VectorXd::Zero(ndim * nNd);
+ dCdp_x = Eigen::VectorXd::Zero(ndim * nNd);
+ dCdf_x = Eigen::VectorXd::Zero(ndim * nNd);
tdCl_x.resize(nNd, Eigen::Vector3d::Zero());
tdCd_x.resize(nNd, Eigen::Vector3d::Zero());
}
-void CoupledAdjoint::buildAdjointMatrix(
- std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>>
- &matrices,
- Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint)
+void CoupledAdjoint::buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint)
{
// Create a list of triplets for the larger matrix
std::deque<Eigen::Triplet<double>> triplets;
@@ -215,31 +215,69 @@ void CoupledAdjoint::run()
tms["0-Total"].start();
tms["1-Derivatives"].start();
tms2["1-adj_inviscid"].start();
- gradientswrtInviscidFlow();
+ gradientsInviscid();
tms2["1-adj_inviscid"].stop();
- tms2["2-adj_mach"].start();
- gradientswrtMachNumber();
- tms2["2-adj_mach"].stop();
- tms2["3-adj_density"].start();
- gradientswrtDensity();
- tms2["3-adj_density"].stop();
- tms2["4-adj_velocity"].start();
- gradientswrtVelocity();
- tms2["4-adj_velocity"].stop();
- tms2["5-adj_dStar"].start();
- gradientswrtDeltaStar();
- tms2["5-adj_dStar"].stop();
- tms2["6-adj_blw"].start();
- gradientswrtBlowingVelocity();
- tms2["6-adj_blw"].stop();
- tms2["7-adj_aoa"].start();
- gradientwrtAoA();
- tms2["7-adj_aoa"].stop();
- tms2["8-adj_mesh"].start();
- gradientwrtMesh();
- tms2["8-adj_mesh"].stop();
- tms["1-Derivatives"].stop();
-
+ tms2["2-adj_viscous"].start();
+ gradientsViscous();
+ tms2["2-adj_viscous"].stop();
+
+ // writeMatrixToFile(dRphi_phi, "dRphi_phi");
+ // writeMatrixToFile(dRphi_rho, "dRphi_rho");
+ // writeMatrixToFile(dRphi_v, "dRphi_v");
+ // writeMatrixToFile(dRphi_dStar, "dRphi_dStar");
+ // writeMatrixToFile(dRphi_blw, "dRphi_blw");
+ // writeMatrixToFile(dRphi_x, "dRphi_x");
+
+ // writeMatrixToFile(dRM_phi, "dRM_phi");
+ // writeMatrixToFile(dRM_M, "dRM_M");
+ // writeMatrixToFile(dRM_rho, "dRM_rho");
+ // writeMatrixToFile(dRM_v, "dRM_v");
+ // writeMatrixToFile(dRM_dStar, "dRM_dStar");
+ // writeMatrixToFile(dRM_blw, "dRM_blw");
+ // writeMatrixToFile(dRM_x, "dRM_x");
+
+ // writeMatrixToFile(dRrho_phi, "dRrho_phi");
+ // writeMatrixToFile(dRrho_M, "dRrho_M");
+ // writeMatrixToFile(dRrho_rho, "dRrho_rho");
+ // writeMatrixToFile(dRrho_v, "dRrho_v");
+ // writeMatrixToFile(dRrho_dStar, "dRrho_dStar");
+ // writeMatrixToFile(dRrho_blw, "dRrho_blw");
+ // writeMatrixToFile(dRrho_x, "dRrho_x");
+
+ // writeMatrixToFile(dRv_phi, "dRv_phi");
+ // writeMatrixToFile(dRv_M, "dRv_M");
+ // writeMatrixToFile(dRv_rho, "dRv_rho");
+ // writeMatrixToFile(dRv_v, "dRv_v");
+ // writeMatrixToFile(dRv_dStar, "dRv_dStar");
+ // writeMatrixToFile(dRv_blw, "dRv_blw");
+ // writeMatrixToFile(dRv_x, "dRv_x");
+
+ // writeMatrixToFile(dRdStar_phi, "dRdStar_phi");
+ // writeMatrixToFile(dRdStar_M, "dRdStar_M");
+ // writeMatrixToFile(dRdStar_rho, "dRdStar_rho");
+ // writeMatrixToFile(dRdStar_v, "dRdStar_v");
+ // writeMatrixToFile(dRdStar_dStar, "dRdStar_dStar");
+ // writeMatrixToFile(dRdStar_blw, "dRdStar_blw");
+ // writeMatrixToFile(dRdStar_x, "dRdStar_x");
+
+ // writeMatrixToFile(dRblw_phi, "dRb_phi");
+ // writeMatrixToFile(dRblw_M, "dRb_M");
+ // writeMatrixToFile(dRblw_rho, "dRb_rho");
+ // writeMatrixToFile(dRblw_v, "dRb_v");
+ // writeMatrixToFile(dRblw_dStar, "dRb_dStar");
+ // writeMatrixToFile(dRblw_blw, "dRb_blw");
+ // writeMatrixToFile(dRblw_x, "dRb_x");
+
+ // writeMatrixToFile(dCl_phi, "dCl_phi");
+
+ // writeMatrixToFile(dRphi_AoA, "dRphi_AoA");
+
+ // writeMatrixToFile(dCd_phi, "dCd_phi");
+ // writeMatrixToFile(dCd_M, "dCd_M");
+ // writeMatrixToFile(dCd_rho, "dCd_rho");
+ // writeMatrixToFile(dCd_v, "dCd_v");
+ // writeMatrixToFile(dCd_dStar, "dCd_dStar");
+ // writeMatrixToFile(dCdf_x, "dCdf_x");
transposeMatrices();
tms2["9-build"].start();
@@ -267,6 +305,8 @@ void CoupledAdjoint::run()
buildAdjointMatrix(matrices, A_adjoint);
tms2["9-build"].stop();
+ writeMatrixToFile(A_adjoint, "A_adjointdev.txt");
+
size_t phiIdx = 0;
size_t machIdx = dRphi_phi.cols();
size_t rhoIdx = machIdx + dRM_M.cols();
@@ -477,8 +517,10 @@ void CoupledAdjoint::run()
* Non-zero are: dRphi_phi, dRM_phi, dRrho_phi, dRv_phi
* Zeros are: dRdStar_phi, dRblw_phi, dRx_phi
*/
-void CoupledAdjoint::gradientswrtInviscidFlow()
+void CoupledAdjoint::gradientsInviscid()
{
+
+ int ndim = isol->pbl->nDim;
//**************************************************************************//
// dRphi_phi //
//--------------------------------------------------------------------------//
@@ -511,7 +553,7 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
iReg = 0;
else
throw std::runtime_error(
- "gradientswrtInviscidFlow: Unknown section name");
+ "gradientsInviscid: Unknown section name");
for (size_t iNode = 0; iNode < reg->getnNodes(); ++iNode)
{
tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[reg->nodes[iNode]->row];
@@ -562,722 +604,286 @@ void CoupledAdjoint::gradientswrtInviscidFlow()
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());
-}
-/**
- * @brief Compute all the gradients wrt the Mach number on the surface
- * Non-zero are: dRM_M, dRdStar_M
- * Zeros are: dRphi_M, dRrho_M, dRv_M, dRblw_M, dRx_M
- */
-void CoupledAdjoint::gradientswrtMachNumber()
-{
//**************************************************************************//
// dRM_M //
//**************************************************************************//
dRM_M.setIdentity();
-// //**************************************************************************//
-// // dRdStar_M //
-// //--------------------------------------------------------------------------//
-// // Central finite-difference //
-// //**************************************************************************//
-// std::deque<Eigen::Triplet<double>> T;
-// double delta = 1e-6;
-// size_t i = 0;
-// double saveM = 0.;
-// std::vector<Eigen::VectorXd> ddStar(2,
-// Eigen::VectorXd::Zero(dRdStar_M.cols()));
-// std::vector<double> dCdf(2, 0.);
-// for (auto &body : vsol->bodies)
-// {
-// for (auto §ion : body->sections)
-// {
-// for (auto ® : section->regions)
-// {
-// for (size_t j = 0; j < reg->getnNodes(); ++j)
-// {
-// saveM = reg->Me[j];
-
-// reg->Me[j] = saveM + delta;
-// ddStar[0] = this->runViscous();
-// dCdf[0] = vsol->Cdf;
-// reg->Me[j] = saveM - delta;
-// ddStar[1] = this->runViscous();
-// dCdf[1] = vsol->Cdf;
-
-// reg->Me[j] = saveM;
-
-// for (int k = 0; k < ddStar[0].size(); ++k)
-// T.push_back(Eigen::Triplet<double>(
-// k, i, -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
-// // Collect Cd gradients
-// dCd_M(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
-// ++i;
-// }
-// }
-// }
-// dRdStar_M.setFromTriplets(T.begin(), T.end());
-// dRdStar_M.prune(0.);
-// }
-}
-
-// /**
-// * @brief Compute all the gradients wrt the density on the surface
-// * Non-zero are: dRrho_rho, dRblw_rho
-// * Zeros are: dRphi_rho, dRM_rho, dRv_rho, dRdStar_rho, dRx_rho
-// */
-void CoupledAdjoint::gradientswrtDensity()
-{
-// size_t nNs = 0;
-// for (auto bBC : isol->pbl->bBCs)
-// nNs += bBC->nodes.size(); // Number of surface nodes
-// nNs += 1; // Duplicate stagnation point
-// size_t nEs = 0;
-// for (auto bBC : isol->pbl->bBCs)
-// nEs += bBC->uE.size(); // Number of blowing elements
-
-// //**************************************************************************//
-// // dRrho_rho //
-// //**************************************************************************//
-// dRrho_rho.setIdentity();
-
-// //**************************************************************************//
-// // dRblw_rho //
-// //--------------------------------------------------------------------------//
-// // Analytical //
-// //**************************************************************************//
-// std::deque<Eigen::Triplet<double>> T_blw;
-// int offSetElms = 0;
-// int offSetNodes = 0;
-// for (auto &body : vsol->bodies)
-// {
-// for (auto §ion : body->sections)
-// {
-// for (auto ® : section->regions)
-// {
-// std::vector<std::vector<double>> grad = reg->evalGradwrtRho();
-// for (size_t i = 0; i < grad.size(); ++i)
-// for (size_t j = 0; j < grad[i].size(); ++j)
-// T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
-// -grad[i][j]));
-// offSetElms += reg->getnElms();
-// offSetNodes += reg->getnNodes();
-// }
-// }
-// }
-// dRblw_rho.setFromTriplets(T_blw.begin(), T_blw.end());
-// dRblw_rho.prune(0.);
-}
+ //**************************************************************************//
+ // dRrho_rho //
+ //**************************************************************************//
+ dRrho_rho.setIdentity();
-/**
- * @brief Compute all the gradients wrt the velocity on the surface
- * Non-zeros are: dRv_v, dRdStar_v, dRblw_v
- * Zeros are: dRphi_v, dRM_v, dRrho_v, dRx_v
- */
-void CoupledAdjoint::gradientswrtVelocity()
-{
- // size_t nDim = isol->pbl->nDim; // Number of dimensions of the problem
- // size_t nNs = 0;
- // for (auto bBC : isol->pbl->bBCs)
- // nNs += bBC->nodes.size(); // Number of surface nodes
- // nNs += 1; // Duplicate stagnation point
- // size_t nEs = 0;
- // for (auto bBC : isol->pbl->bBCs)
- // nEs += bBC->uE.size(); // Number of blowing elements
-
- // //**************************************************************************//
- // // dRv_v //
- // //**************************************************************************//
- // dRv_v.setIdentity();
-
- // //**** Get velocity****/
- // Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nNs, nDim);
- // size_t i = 0;
- // int iReg = -1;
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // if (reg->getName() == "wake")
- // iReg = 1;
- // else if (reg->getName() == "upper" || reg->getName() == "lower")
- // iReg = 0;
- // else
- // throw std::runtime_error("Unknown section name");
- // for (size_t j = 0; j < reg->getnNodes(); ++j)
- // {
- // for (size_t idim = 0; idim < nDim; ++idim)
- // v(i, idim) =
- // isol->U[isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row](idim);
- // ++i;
- // }
- // }
- // }
- // }
-
- // //**************************************************************************//
- // // dvt_v //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nNs, nNs * nDim);
- // std::deque<Eigen::Triplet<double>> T_vt;
-
- // i = 0;
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // for (size_t j = 0; j < reg->getnNodes(); ++j)
- // {
- // for (size_t idim = 0; idim < nDim; ++idim)
- // {
- // T_vt.push_back(Eigen::Triplet<double>(
- // i, i * nDim + idim,
- // 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim)));
- // }
- // ++i;
- // }
- // }
- // }
- // }
- // dVt_v.setFromTriplets(T_vt.begin(), T_vt.end());
- // dVt_v.prune(0.);
-
- // //**************************************************************************//
- // // dRdStar_vt //
- // //--------------------------------------------------------------------------//
- // // Central finite-difference //
- // //**************************************************************************//
- // std::deque<Eigen::Triplet<double>> T;
- // double delta = 1e-6;
- // i = 0;
- // double saveM = 0.;
- // std::vector<Eigen::VectorXd> dvt(2, Eigen::VectorXd::Zero(dRdStar_M.cols()));
- // std::vector<double> dCdf(2, 0.);
- // Eigen::VectorXd dCd_vt = Eigen::VectorXd::Zero(dRM_M.cols());
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // for (size_t j = 0; j < reg->getnNodes(); ++j)
- // {
- // saveM = reg->vt[j];
-
- // reg->vt[j] = saveM + delta;
- // dvt[0] = this->runViscous();
- // dCdf[0] = vsol->Cdf;
- // reg->vt[j] = saveM - delta;
- // dvt[1] = this->runViscous();
- // dCdf[1] = vsol->Cdf;
-
- // reg->vt[j] = saveM;
-
- // for (int k = 0; k < dvt[0].size(); ++k)
- // T.push_back(Eigen::Triplet<double>(
- // k, i, -(dvt[0][k] - dvt[1][k]) / (2. * delta)));
- // dCd_vt(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
- // ++i;
- // }
- // }
- // }
- // }
- // Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt(dRdStar_M.rows(),
- // dRdStar_M.cols());
- // dRdStar_vt.setFromTriplets(T.begin(), T.end());
- // dRdStar_vt.prune(0.);
-
- // //**************************************************************************//
- // // dRdStar_v //
- // //**************************************************************************//
- // dRdStar_v = dRdStar_vt * dVt_v;
- // dRdStar_v.prune(0.);
-
- // //**************************************************************************//
- // // dCdf_v //
- // //**************************************************************************//
- // dCd_v =
- // dCd_vt.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim]
-
- // //**************************************************************************//
- // // dRblw_vt //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // std::deque<Eigen::Triplet<double>> T_blw;
- // int offSetElms = 0;
- // int offSetNodes = 0;
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // std::vector<std::vector<double>> grad = reg->evalGradwrtVt();
- // for (size_t i = 0; i < grad.size(); ++i)
- // for (size_t j = 0; j < grad[i].size(); ++j)
- // T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
- // -grad[i][j]));
- // offSetElms += reg->getnElms();
- // offSetNodes += reg->getnNodes();
- // }
- // }
- // }
- // Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt(nEs, nNs);
- // dRblw_vt.setFromTriplets(T_blw.begin(), T_blw.end());
- // dRblw_vt.prune(0.);
-
- // //**************************************************************************//
- // // dRblw_v //
- // //**************************************************************************//
- // dRblw_v = dRblw_vt * dVt_v;
- // dRblw_v.prune(0.);
-}
+ //**************************************************************************//
+ // dRv_v //
+ //**************************************************************************//
+ dRv_v.setIdentity();
-/**
- * @brief Compute all the gradients wrt the displacement thickness on the
- * surface Non-zero are: dRdStar_dStar, dRblw_dStar Zeros are: dRphi_dStar,
- * dRM_dStar, dRrho_dStar, dRv_dStar, dRx_dStar
- */
-void CoupledAdjoint::gradientswrtDeltaStar()
-{
//**************************************************************************//
- // dRdStar_dStar //
+ // dCl_x, dCd_x //
//--------------------------------------------------------------------------//
- // Central finite-difference //
+ // Analytical //
//**************************************************************************//
- // std::deque<Eigen::Triplet<double>> T;
- // double delta = 1e-6;
- // double saveDStar = 0.;
- // std::vector<Eigen::VectorXd> ddstar(
- // 2, Eigen::VectorXd::Zero(dRdStar_dStar.cols()));
- // std::vector<double> dCdf(2, 0.);
- // size_t i = 0;
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // for (size_t j = 0; j < reg->getnNodes(); ++j)
- // {
- // saveDStar = reg->deltaStarExt[j];
-
- // reg->deltaStarExt[j] = saveDStar + delta;
- // ddstar[0] = this->runViscous();
- // dCdf[0] = vsol->Cdf;
- // reg->deltaStarExt[j] = saveDStar - delta;
- // ddstar[1] = this->runViscous();
- // dCdf[1] = vsol->Cdf;
-
- // reg->deltaStarExt[j] = saveDStar;
-
- // for (int k = 0; k < ddstar[0].size(); ++k)
- // T.push_back(Eigen::Triplet<double>(
- // k, i, (ddstar[0][k] - ddstar[1][k]) / (2. * delta)));
- // dCd_dStar(i) = (dCdf[0] - dCdf[1]) / (2. * delta);
- // ++i;
- // }
- // }
- // }
- // }
- // Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar_dum(
- // dRdStar_dStar.rows(), dRdStar_dStar.cols());
- // // RdStar = dStar - f(dStar) = 0
- // // dRdStar_dStar = 1 - df_dStar (minus is 3 lines below and not in the
- // // triplets)
- // dRdStar_dStar_dum.setFromTriplets(T.begin(), T.end());
- // dRdStar_dStar.setIdentity();
- // dRdStar_dStar -= dRdStar_dStar_dum;
- // dRdStar_dStar.prune(0.);
-
- // //**************************************************************************//
- // // dRblw_dStar //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // std::deque<Eigen::Triplet<double>> T_blw;
- // int offSetElms = 0;
- // int offSetNodes = 0;
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // std::vector<std::vector<double>> grad = reg->evalGradwrtDeltaStar();
- // for (size_t i = 0; i < grad.size(); ++i)
- // for (size_t j = 0; j < grad[i].size(); ++j)
- // T_blw.push_back(Eigen::Triplet<double>(i + offSetElms, j + offSetNodes,
- // -grad[i][j]));
- // offSetElms += reg->getnElms();
- // offSetNodes += reg->getnNodes();
- // }
- // }
- // }
- // dRblw_dStar.setFromTriplets(T_blw.begin(), T_blw.end());
- // dRblw_dStar.prune(0.);
-}
+ std::vector<std::vector<double>> dCx_x =
+ iadjoint->computeGradientCoefficientsMesh();
+ dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size());
+ dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
-/**
- * @brief Compute all the gradients wrt the blowing velocity
- * Non-zero are: dRphi_blw, dRblw_blw
- * Zeros are: dRM_blw, dRrho_blw, dRv_blw, dRdStar_blw, dRx_blw
- */
-void CoupledAdjoint::gradientswrtBlowingVelocity()
-{
//**************************************************************************//
- // dRphi_blw //
+ // dRphi_x //
//--------------------------------------------------------------------------//
// Analytical //
//**************************************************************************//
- // std::deque<Eigen::Triplet<double>> T;
-
- // size_t jj = 0;
- // int iReg = 0;
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // if (reg->getName() == "wake")
- // iReg = 1;
- // else if (reg->getName() == "upper" || reg->getName() == "lower")
- // iReg = 0;
- // else
- // throw std::runtime_error("Unknown section name");
- // for (size_t iElm = 0; iElm < reg->getnElms(); ++iElm)
- // {
- // auto blowingElement = isol->pbl->bBCs[iReg]->uE[reg->no[iElm]].first;
- // Eigen::VectorXd be =
- // dart::BlowingResidual::buildGradientBlowing(*blowingElement);
- // for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++)
- // {
- // tbox::Node *nodi = blowingElement->nodes[ii];
- // T.push_back(
- // Eigen::Triplet<double>(isol->getRow(nodi->row), jj, be(ii)));
- // }
- // ++jj;
- // }
- // }
- // }
- // }
- // dRphi_blw.setFromTriplets(T.begin(), T.end());
- // dRphi_blw.prune(0.);
- // dRphi_blw.makeCompressed();
-
- // /**** dRblw_blw ****/
- // dRblw_blw.setIdentity();
-}
+ dRphi_x = iadjoint->getRu_X();
+ //**************************************************************************//
+ // dRx_x //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ dRx_x = iadjoint->getRx_X();
+
+ //**************************************************************************//
+ // dRM_x, dRrho_x, dRv_x //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+
+ std::deque<Eigen::Triplet<double>> T_dM_x;
+ std::deque<Eigen::Triplet<double>> T_drho_x;
+ std::deque<Eigen::Triplet<double>> T_dv_x;
+
+ jj = 0;
+ iReg = 0;
+
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ if (reg->getName() == "wake")
+ iReg = 1;
+ else if (reg->getName() == "upper" || reg->getName() == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error(
+ "gradientsInviscid: Unknown section name");
+ for (size_t iNode = 0; iNode < reg->getnNodes(); ++iNode)
+ {
+ tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[reg->nodes[iNode]->row];
+ double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size());
+ for (auto e : pbl->fluid->neMap[srfNode])
+ {
+ // dv/dx = dv/dgradPhi * dgradPhi/dx = dgradPhi/dx = (-J^(-1) *
+ // dJ/dx_j)*grad_xPhi dM/dx = dM/dgradPhi * dgradPhi/dx = dM/dgradPhi *
+ // gradPhi * dv/dx drho/dx = drho/dgradPhi * dgradPhi/dx = drho/dgradPhi
+ // * gradPhi * dv/dx
+ for (size_t inodOnElm = 0; inodOnElm < e->nodes.size(); inodOnElm++)
+ for (int idim = 0; idim < ndim; ++idim)
+ {
+ Eigen::VectorXd gradV = (-e->getJinv(0)) *
+ e->getGradJ(0)[inodOnElm * ndim + idim] *
+ e->computeGradient(isol->phi, 0);
+ double gradM = pbl->fluid->mach->evalGrad(*e, isol->phi, 0) *
+ e->computeGradient(isol->phi, 0).transpose() * gradV;
+ double gradrho = pbl->fluid->rho->evalGrad(*e, isol->phi, 0) *
+ e->computeGradient(isol->phi, 0).transpose() *
+ gradV;
+ T_dM_x.push_back(Eigen::Triplet<double>(
+ jj, e->nodes[inodOnElm]->row * ndim + idim, -gradM / nAdjElms));
+ T_drho_x.push_back(Eigen::Triplet<double>(
+ jj, e->nodes[inodOnElm]->row * ndim + idim,
+ -gradrho / nAdjElms));
+ for (int idim2 = 0; idim2 < ndim; idim2++)
+ T_dv_x.push_back(Eigen::Triplet<double>(
+ jj * ndim + idim2, e->nodes[inodOnElm]->row * ndim + idim,
+ -gradV(idim2) / nAdjElms));
+ }
+ }
+ ++jj;
+ }
+ }
+ }
+ }
+ dRM_x.setFromTriplets(T_dM_x.begin(), T_dM_x.end());
+ dRrho_x.setFromTriplets(T_drho_x.begin(), T_drho_x.end());
+ dRv_x.setFromTriplets(T_dv_x.begin(), T_dv_x.end());
+ dRM_x.prune(0.);
+ dRrho_x.prune(0.);
+ dRv_x.prune(0.);
-/**
- * @brief Compute all the gradients wrt the angle of attack
- * Non-zero are: dRphi_AoA
- * Zeros are: dRM_AoA, dRrho_AoA, dRv_AoA, dRdStar_AoA, dRblw_AoA, dRx_AoA
- */
-void CoupledAdjoint::gradientwrtAoA()
-{
//**************************************************************************//
// dCl_Aoa dCd_Aoa //
//--------------------------------------------------------------------------//
// Analytical //
//**************************************************************************//
- // std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa();
- // dCl_AoA = dCx_AoA[0];
- // dCd_AoA = dCx_AoA[1];
-
- // //**************************************************************************//
- // // dRphi_Aoa //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // dRphi_AoA = iadjoint->getRu_A();
-}
+ std::vector<double> dCx_AoA = iadjoint->computeGradientCoefficientsAoa();
+ dCl_AoA = dCx_AoA[0];
+ dCd_AoA = dCx_AoA[1];
-/**
- * @brief Compute all the gradients wrt the mesh coordinates
- * Non-zero are: dRphi_x, dRM_x, dRrho_x, dRv_x, dRdStar_x, dRblw_x, dRx_x
- * Zeros are: -
- */
-void CoupledAdjoint::gradientwrtMesh()
-{
+ //**************************************************************************//
+ // dRphi_Aoa //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ dRphi_AoA = iadjoint->getRu_A();
+
+ //**************************************************************************//
+ // dRphi_blw //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ std::deque<Eigen::Triplet<double>> T;
- // int nDim = isol->pbl->nDim;
-
- // //**************************************************************************//
- // // dCl_x, dCd_x //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // std::vector<std::vector<double>> dCx_x =
- // iadjoint->computeGradientCoefficientsMesh();
- // dCl_x = Eigen::Map<Eigen::VectorXd>(dCx_x[0].data(), dCx_x[0].size());
- // dCdp_x = Eigen::Map<Eigen::VectorXd>(dCx_x[1].data(), dCx_x[1].size());
-
- // //**************************************************************************//
- // // dRphi_x //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // dRphi_x = iadjoint->getRu_X();
- // //**************************************************************************//
- // // dRx_x //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // dRx_x = iadjoint->getRx_X();
-
- // //**************************************************************************//
- // // dRdStar_x, dCdf_dx //
- // //--------------------------------------------------------------------------//
- // // Central finite-difference //
- // //--------------------------------------------------------------------------//
- // // We don't do dRdStar_x = dRdStar_loc * dloc_dx because Cdf = f(xi(x), x) //
- // // and we would have to compute dCdf/dx and dCdf/dloc. Same for deltaStar. //
- // // It is more computationally expansive to compute //
- // // ddeltaStar/dx, ddeltaStar/dloc than ddeltaStar_dx, ddeltaStar_dy //
- // //**************************************************************************//
- // std::deque<Eigen::Triplet<double>> T_dStar_x;
- // double delta = 1e-8;
- // int iReg = 0;
- // std::vector<Eigen::VectorXd> ddStar(2,
- // Eigen::VectorXd::Zero(dRdStar_M.cols()));
- // std::vector<double> dCdf(2, 0.);
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // if (reg->getName() == "wake")
- // iReg = 1;
- // else if (reg->getName() == "upper" || reg->getName() == "lower")
- // iReg = 0;
- // else
- // throw std::runtime_error("Unknown section name");
- // for (size_t j = 0; j < reg->getnNodes(); ++j)
- // {
- // int rowj = isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row;
- // // x
- // double saveX = reg->nodes[j]->pos[0];
-
- // reg->nodes[j]->pos[0] = saveX + delta;
- // section->computeLocalCoordinates(reg);
- // for (auto &nod : reg->nodes)
- // nod->xiExt = nod->xi;
- // ddStar[0] = this->runViscous();
- // dCdf[0] = vsol->Cdf;
- // reg->nodes[j]->pos[0] = saveX - delta;
- // section->computeLocalCoordinates(reg);
- // for (auto &nod : reg->nodes)
- // nod->xiExt = nod->xi;
- // ddStar[1] = this->runViscous();
- // dCdf[1] = vsol->Cdf;
-
- // reg->nodes[j]->pos[0] = saveX;
- // section->computeLocalCoordinates(reg);
- // for (auto &nod : reg->nodes)
- // nod->xiExt = nod->xi;
-
- // for (int k = 0; k < ddStar[0].size(); ++k)
- // T_dStar_x.push_back(Eigen::Triplet<double>(
- // k, rowj * isol->pbl->nDim + 0,
- // -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
- // dCdf_x(rowj * isol->pbl->nDim + 0) += (dCdf[0] - dCdf[1]) / (2. * delta);
- // // y
- // double saveY = reg->nodes[j]->pos[1];
-
- // reg->nodes[j]->pos[1] = saveY + delta;
- // section->computeLocalCoordinates(reg);
- // for (auto &nod : reg->nodes)
- // nod->xiExt = nod->xi;
- // ddStar[0] = this->runViscous();
- // dCdf[0] = vsol->Cdf;
- // reg->nodes[j]->pos[1] = saveY - delta;
- // section->computeLocalCoordinates(reg);
- // for (auto &nod : reg->nodes)
- // nod->xiExt = nod->xi;
- // ddStar[1] = this->runViscous();
- // dCdf[1] = vsol->Cdf;
-
- // reg->nodes[j]->pos[1] = saveY;
- // section->computeLocalCoordinates(reg);
- // for (auto &nod : reg->nodes)
- // nod->xiExt = nod->xi;
-
- // for (int k = 0; k < ddStar[0].size(); ++k)
- // T_dStar_x.push_back(Eigen::Triplet<double>(
- // k, rowj * isol->pbl->nDim + 1,
- // -(ddStar[0][k] - ddStar[1][k]) / (2. * delta)));
- // dCdf_x(rowj * isol->pbl->nDim + 1) += (dCdf[0] - dCdf[1]) / (2. * delta);
- // }
- // }
- // }
- // }
- // dRdStar_x.setFromTriplets(T_dStar_x.begin(), T_dStar_x.end());
- // dRdStar_x.prune(0.);
-
- // //**************************************************************************//
- // // dRblw_x //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // std::deque<Eigen::Triplet<double>> T_blw_x;
- // int offSetElms = 0;
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // if (reg->getName() == "wake")
- // iReg = 1;
- // else if (reg->getName() == "upper" || reg->getName() == "lower")
- // iReg = 0;
- // else
- // throw std::runtime_error("Unknown section name");
- // std::vector<std::vector<double>> gradX = reg->evalGradwrtX();
- // std::vector<std::vector<double>> gradY = reg->evalGradwrtY();
- // for (size_t i = 0; i < gradX.size(); ++i)
- // for (size_t j = 0; j < gradX[i].size(); ++j)
- // {
- // int rowj = isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row;
- // T_blw_x.push_back(Eigen::Triplet<double>(
- // i + offSetElms, rowj * isol->pbl->nDim + 0, -gradX[i][j]));
- // T_blw_x.push_back(Eigen::Triplet<double>(
- // i + offSetElms, rowj * isol->pbl->nDim + 1, -gradY[i][j]));
- // }
- // offSetElms += reg->getnElms();
- // }
- // }
- // }
- // dRblw_x.setFromTriplets(T_blw_x.begin(), T_blw_x.end());
- // dRblw_x.prune(0.);
-
- // //**************************************************************************//
- // // dRM_x, dRrho_x, dRv_x //
- // //--------------------------------------------------------------------------//
- // // Analytical //
- // //**************************************************************************//
- // auto pbl = isol->pbl;
-
- // std::deque<Eigen::Triplet<double>> T_dM_x;
- // std::deque<Eigen::Triplet<double>> T_drho_x;
- // std::deque<Eigen::Triplet<double>> T_dv_x;
-
- // int jj = 0;
- // iReg = 0;
-
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (auto ® : section->regions)
- // {
- // if (reg->getName() == "wake")
- // iReg = 1;
- // else if (reg->getName() == "upper" || reg->getName() == "lower")
- // iReg = 0;
- // else
- // throw std::runtime_error(
- // "gradientswrtInviscidFlow: Unknown section name");
- // for (size_t iNode = 0; iNode < reg->getnNodes(); ++iNode)
- // {
- // tbox::Node *srfNode = pbl->bBCs[iReg]->nodes[reg->nodes[iNode]->row];
- // double nAdjElms = static_cast<double>(pbl->fluid->neMap[srfNode].size());
- // for (auto e : pbl->fluid->neMap[srfNode])
- // {
- // // dv/dx = dv/dgradPhi * dgradPhi/dx = dgradPhi/dx = (-J^(-1) *
- // // dJ/dx_j)*grad_xPhi dM/dx = dM/dgradPhi * dgradPhi/dx = dM/dgradPhi *
- // // gradPhi * dv/dx drho/dx = drho/dgradPhi * dgradPhi/dx = drho/dgradPhi
- // // * gradPhi * dv/dx
- // for (size_t inodOnElm = 0; inodOnElm < e->nodes.size(); inodOnElm++)
- // for (int idim = 0; idim < nDim; ++idim)
- // {
- // Eigen::VectorXd gradV = (-e->getJinv(0)) *
- // e->getGradJ(0)[inodOnElm * nDim + idim] *
- // e->computeGradient(isol->phi, 0);
- // double gradM = pbl->fluid->mach->evalGrad(*e, isol->phi, 0) *
- // e->computeGradient(isol->phi, 0).transpose() * gradV;
- // double gradrho = pbl->fluid->rho->evalGrad(*e, isol->phi, 0) *
- // e->computeGradient(isol->phi, 0).transpose() *
- // gradV;
- // T_dM_x.push_back(Eigen::Triplet<double>(
- // jj, e->nodes[inodOnElm]->row * nDim + idim, -gradM / nAdjElms));
- // T_drho_x.push_back(Eigen::Triplet<double>(
- // jj, e->nodes[inodOnElm]->row * nDim + idim,
- // -gradrho / nAdjElms));
- // for (int idim2 = 0; idim2 < nDim; idim2++)
- // T_dv_x.push_back(Eigen::Triplet<double>(
- // jj * nDim + idim2, e->nodes[inodOnElm]->row * nDim + idim,
- // -gradV(idim2) / nAdjElms));
- // }
- // }
- // ++jj;
- // }
- // }
- // }
- // }
- // dRM_x.setFromTriplets(T_dM_x.begin(), T_dM_x.end());
- // dRrho_x.setFromTriplets(T_drho_x.begin(), T_drho_x.end());
- // dRv_x.setFromTriplets(T_dv_x.begin(), T_dv_x.end());
- // dRM_x.prune(0.);
- // dRrho_x.prune(0.);
- // dRv_x.prune(0.);
+ jj = 0;
+ iReg = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ if (reg->getName() == "wake")
+ iReg = 1;
+ else if (reg->getName() == "upper" || reg->getName() == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ for (size_t iElm = 0; iElm < reg->getnElms(); ++iElm)
+ {
+ auto blowingElement = isol->pbl->bBCs[iReg]->uE[reg->no[iElm]].first;
+ Eigen::VectorXd be =
+ dart::BlowingResidual::buildGradientBlowing(*blowingElement);
+ for (size_t ii = 0; ii < blowingElement->nodes.size(); ii++)
+ {
+ tbox::Node *nodi = blowingElement->nodes[ii];
+ T.push_back(
+ Eigen::Triplet<double>(isol->getRow(nodi->row), jj, be(ii)));
+ }
+ ++jj;
+ }
+ }
+ }
+ }
+ dRphi_blw.setFromTriplets(T.begin(), T.end());
+ dRphi_blw.prune(0.);
+ dRphi_blw.makeCompressed();
}
-/**
- * @brief Run the viscous solver and extract the deltaStar
- * on all the stations
- *
- * @return Eigen::VectorXd deltaStar
- */
-Eigen::VectorXd CoupledAdjoint::runViscous()
+void CoupledAdjoint::gradientsViscous()
{
- // int exitCode = vsol->run();
- // if (exitCode != 0 && vsol->verbose > 0)
- // std::cout << "vsol terminated with exit code " << exitCode << std::endl;
-
- // // Extract deltaStar
- // std::vector<Eigen::VectorXd> dStar_p;
-
- // // Extract parts in a loop
- // for (auto &body : vsol->bodies)
- // {
- // for (auto §ion : body->sections)
- // {
- // for (size_t i = 0; i < section->regions.size(); ++i)
- // {
- // std::vector<double> deltaStar_vec = section->regions[i]->deltaStar;
- // Eigen::VectorXd deltaStar_eigen =
- // Eigen::Map<Eigen::VectorXd>(deltaStar_vec.data(), deltaStar_vec.size());
- // dStar_p.push_back(deltaStar_eigen);
- // }
- // }
- // }
- // // Calculate the total size
- // int nNs = 0;
- // for (const auto &p : dStar_p)
- // nNs += p.size();
-
- // // Concatenate the vectors
- // Eigen::VectorXd deltaStar(nNs);
- // int idx = 0;
- // for (const auto &p : dStar_p)
- // {
- // deltaStar.segment(idx, p.size()) = p;
- // idx += p.size();
- // }
- // return deltaStar;
-}
+ size_t ndim = isol->pbl->nDim; // Number of dimensions of the problem
+ size_t nnodes = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nnodes += bBC->nodes.size(); // Number of surface nodes
+ nnodes += 1; // Duplicate stagnation point
+ size_t nelms = 0;
+ for (auto bBC : isol->pbl->bBCs)
+ nelms += bBC->uE.size(); // Number of blowing elements
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_vt_dummy(nnodes, nnodes);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_vt_dummy(nelms, nnodes);
+ Eigen::VectorXd dCd_vt_dummy(nnodes);
+
+ std::map<std::string, Eigen::SparseMatrix<double, Eigen::RowMajor> *> adjointMatrixMap;
+ std::map<std::string, Eigen::VectorXd *> adjointVectorMap;
+ adjointMatrixMap["dRdStar_M"] = &dRdStar_M;
+ adjointMatrixMap["dRdStar_vt"] = &dRdStar_vt_dummy;
+ adjointMatrixMap["dRdStar_rho"] = &dRdStar_rho;
+ adjointMatrixMap["dRdStar_dStar"] = &dRdStar_dStar;
+ adjointMatrixMap["dRdStar_blw"] = &dRdStar_blw;
+ adjointMatrixMap["dRdStar_x"] = &dRdStar_x;
+ adjointMatrixMap["dRblw_M"] = &dRblw_M;
+ adjointMatrixMap["dRblw_vt"] = &dRblw_vt_dummy;
+ adjointMatrixMap["dRblw_rho"] = &dRblw_rho;
+ adjointMatrixMap["dRblw_dStar"] = &dRblw_dStar;
+ adjointMatrixMap["dRblw_blw"] = &dRblw_blw;
+ adjointMatrixMap["dRblw_x"] = &dRblw_x;
+ adjointVectorMap["dCdf_M"] = &dCd_M;
+ adjointVectorMap["dCdf_rho"] = &dCd_rho;
+ adjointVectorMap["dCdf_vt"] = &dCd_vt_dummy;
+ adjointVectorMap["dCdf_dStar"] = &dCd_dStar;
+ adjointVectorMap["dCdf_blw"] = &dCd_dStar;
+ adjointVectorMap["dCdf_x"] = &dCdf_x;
+
+ // Run viscous adjoint
+ vadjoint->solve();
+ vadjoint->evaluate(adjointMatrixMap, adjointVectorMap);
+
+ for (auto &mat : adjointMatrixMap)
+ {
+ mat.second->prune(0.);
+ mat.second->makeCompressed();
+ }
+
+ //**** Get velocity****/
+ Eigen::MatrixXd v = Eigen::MatrixXd::Zero(nnodes, ndim);
+ size_t i = 0;
+ int iReg = -1;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ if (reg->getName() == "wake")
+ iReg = 1;
+ else if (reg->getName() == "upper" || reg->getName() == "lower")
+ iReg = 0;
+ else
+ throw std::runtime_error("Unknown section name");
+ for (size_t j = 0; j < reg->getnNodes(); ++j)
+ {
+ for (size_t idim = 0; idim < ndim; ++idim)
+ v(i, idim) =
+ isol->U[isol->pbl->bBCs[iReg]->nodes[reg->nodes[j]->row]->row](idim);
+ ++i;
+ }
+ }
+ }
+ }
+ //**************************************************************************//
+ // dvt_v //
+ //--------------------------------------------------------------------------//
+ // Analytical //
+ //**************************************************************************//
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dVt_v(nnodes, nnodes * ndim);
+ std::deque<Eigen::Triplet<double>> T_vt;
+
+ i = 0;
+ for (auto &body : vsol->bodies)
+ {
+ for (auto §ion : body->sections)
+ {
+ for (auto ® : section->regions)
+ {
+ for (size_t j = 0; j < reg->getnNodes(); ++j)
+ {
+ for (size_t idim = 0; idim < ndim; ++idim)
+ {
+ T_vt.push_back(Eigen::Triplet<double>(
+ i, i * ndim + idim,
+ 1 / std::sqrt(v(i, 0) * v(i, 0) + v(i, 1) * v(i, 1)) * v(i, idim)));
+ }
+ ++i;
+ }
+ }
+ }
+ }
+ dVt_v.setFromTriplets(T_vt.begin(), T_vt.end());
+ dVt_v.prune(0.);
+
+ dRblw_v = dRblw_vt_dummy * dVt_v;
+ dRblw_v.prune(0.);
+ dRblw_v.makeCompressed();
+ dRdStar_v = dRdStar_vt_dummy * dVt_v;
+ dRdStar_v.prune(0.);
+ dRdStar_v.makeCompressed();
+
+ dCd_v = dCd_vt_dummy.transpose() * dVt_v; // [1xnNs] x [nNs x nNs*nDim] = [1 x nNs*nDim]
+}
/**
* @brief Transpose all matrices of the adjoint problem included in the system.
* Not the others.
@@ -1336,7 +942,7 @@ void CoupledAdjoint::transposeMatrices()
*
* The matrices are written in a file with the following format:
* - The first line contains the column headers
- * - The following lines contain the matrix data
+ * - The following linelms contain the matrix data
* - The matrix data is written with a precision of 12 digits
*/
void CoupledAdjoint::writeMatrixToFile(const Eigen::MatrixXd &matrix,
diff --git a/blast/src/blCoupledAdjoint.h b/blast/src/blCoupledAdjoint.h
index e4e3165..35a8c09 100644
--- a/blast/src/blCoupledAdjoint.h
+++ b/blast/src/blCoupledAdjoint.h
@@ -18,6 +18,7 @@
#define BLCOUPLEDADJOINT_H
#include "blast.h"
+#include "blAdjoint.h"
#include "dart.h"
#include "wObject.h"
#include "wTimers.h"
@@ -38,7 +39,9 @@ class BLAST_API CoupledAdjoint : public fwk::wSharedObject
{
public:
- std::shared_ptr<dart::Adjoint> iadjoint; ///< Adjoint solver
+ std::shared_ptr<dart::Adjoint> iadjoint; ///< Inviscid adjoint solver
+ std::shared_ptr<blast::blAdjoint> vadjoint; ///< Viscous adjoint solver
+
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)
@@ -48,155 +51,110 @@ public:
double tdCd_AoA; ///< Total erivative of the drag coefficient with respect to
///< the angle of attack
- std::vector<Eigen::Vector3d>
- tdCl_x; ///< Total derivative of the lift coefficient with respect to the
- ///< mesh coordinates
- std::vector<Eigen::Vector3d>
- tdCd_x; ///< Total derivative of the drag coefficient with respect to the
- ///< mesh coordinates
+ std::vector<Eigen::Vector3d> tdCl_x; ///< Total derivative of the lift coefficient with respect to the
+ ///< mesh coordinates
+ std::vector<Eigen::Vector3d> tdCd_x; ///< Total derivative of the drag coefficient with respect to the
+ ///< mesh coordinates
private:
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRphi_phi; ///< Inviscid flow Jacobian
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRphi_M; ///< Derivative of the inviscid flow residual with respect to the
- ///< Mach number
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRphi_v; ///< Derivative of the inviscid flow residual with respect to the
- ///< velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRphi_rho; ///< Derivative of the inviscid flow residual with respect to
- ///< the density
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRphi_dStar; ///< Derivative of the inviscid flow residual with respect to
- ///< delta*
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRphi_blw; ///< Derivative of the inviscid flow residual with respect to
- ///< the blowing velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRphi_x; ///< Derivative of the potential residual with respect to the
- ///< mesh coordinates
-
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRM_phi; ///< Derivative of the Mach number residual with respect to the
- ///< inviscid flow
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRM_M; ///< Derivative of the Mach number residual with respect to the
- ///< Mach number
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRM_v; ///< Derivative of the Mach number residual with respect to the
- ///< velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRM_rho; ///< Derivative of the Mach number residual with respect to the
- ///< density
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRM_dStar; ///< Derivative of the Mach number residual with respect to
- ///< delta*
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRM_blw; ///< Derivative of the Mach number residual with respect to the
- ///< blowing velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRM_x; ///< Derivative of the Mach number residual with respect to the
- ///< mesh coordinates
-
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRrho_phi; ///< Derivative of the density residual with respect to the
- ///< inviscid flow
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRrho_M; ///< Derivative of the density residual with respect to the Mach
- ///< number
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRrho_v; ///< Derivative of the density residual with respect to the
- ///< velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRrho_rho; ///< Derivative of the density residual with respect to the
- ///< density
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRrho_dStar; ///< Derivative of the density residual with respect to
- ///< delta*
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRrho_blw; ///< Derivative of the density residual with respect to the
- ///< blowing velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRrho_x; ///< Derivative of the density residual with respect to the mesh
- ///< coordinates
-
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRv_phi; ///< Derivative of the velocity residual with respect to the
- ///< inviscid flow
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRv_M; ///< Derivative of the velocity residual with respect to the Mach
- ///< number
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRv_v; ///< Derivative of the velocity residual with respect to the
- ///< velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRv_rho; ///< Derivative of the velocity residual with respect to the
- ///< density
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRv_dStar; ///< Derivative of the velocity residual with respect to delta*
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRv_blw; ///< Derivative of the velocity residual with respect to the
- ///< blowing velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRv_x; ///< Derivative of the velocity residual with respect to the mesh
- ///< coordinates
-
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRdStar_phi; ///< Derivative of delta* residual with respect to the
- ///< inviscid flow
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRdStar_M; ///< Derivative of delta* residual with respect to the Mach
- ///< number
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRdStar_v; ///< Derivative of delta* residual with respect to the velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRdStar_rho; ///< Derivative of delta* residual with respect to the
- ///< density
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRdStar_dStar; ///< Derivative of delta* residual with respect to delta*
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRdStar_blw; ///< Derivative of delta* residual with respect to the
- ///< blowing velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRdStar_x; ///< Derivative of delta* residual with respect to the mesh
- ///< coordinates
-
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRblw_phi; ///< Derivative of the blowing velocity residual with respect
- ///< to the inviscid flow
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRblw_M; ///< Derivative of the blowing velocity residual with respect to
- ///< the Mach number
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRblw_v; ///< Derivative of the blowing velocity residual with respect to
- ///< the velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRblw_rho; ///< Derivative of the blowing velocity residual with respect
- ///< to the density
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRblw_dStar; ///< Derivative of the blowing velocity residual with respect
- ///< to delta*
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRblw_blw; ///< Derivative of the blowing velocity residual with respect
- ///< to the blowing velocity
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRblw_x; ///< Derivative of the blowing velocity residual with respect to
- ///< the mesh coordinates
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_phi; ///< Inviscid flow Jacobian
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_M; ///< Derivative of the inviscid flow residual with respect to the
+ ///< Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_v; ///< Derivative of the inviscid flow residual with respect to the
+ ///< velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_rho; ///< Derivative of the inviscid flow residual with respect to
+ ///< the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_dStar; ///< Derivative of the inviscid flow residual with respect to
+ ///< delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_blw; ///< Derivative of the inviscid flow residual with respect to
+ ///< the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRphi_x; ///< Derivative of the potential residual with respect to the
+ ///< mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_phi; ///< Derivative of the Mach number residual with respect to the
+ ///< inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_M; ///< Derivative of the Mach number residual with respect to the
+ ///< Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_v; ///< Derivative of the Mach number residual with respect to the
+ ///< velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_rho; ///< Derivative of the Mach number residual with respect to the
+ ///< density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_dStar; ///< Derivative of the Mach number residual with respect to
+ ///< delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_blw; ///< Derivative of the Mach number residual with respect to the
+ ///< blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRM_x; ///< Derivative of the Mach number residual with respect to the
+ ///< mesh coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_phi; ///< Derivative of the density residual with respect to the
+ ///< inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_M; ///< Derivative of the density residual with respect to the Mach
+ ///< number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_v; ///< Derivative of the density residual with respect to the
+ ///< velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_rho; ///< Derivative of the density residual with respect to the
+ ///< density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_dStar; ///< Derivative of the density residual with respect to
+ ///< delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_blw; ///< Derivative of the density residual with respect to the
+ ///< blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRrho_x; ///< Derivative of the density residual with respect to the mesh
+ ///< coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_phi; ///< Derivative of the velocity residual with respect to the
+ ///< inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_M; ///< Derivative of the velocity residual with respect to the Mach
+ ///< number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_v; ///< Derivative of the velocity residual with respect to the
+ ///< velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_rho; ///< Derivative of the velocity residual with respect to the
+ ///< density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_dStar; ///< Derivative of the velocity residual with respect to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_blw; ///< Derivative of the velocity residual with respect to the
+ ///< blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRv_x; ///< Derivative of the velocity residual with respect to the mesh
+ ///< coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_phi; ///< Derivative of delta* residual with respect to the
+ ///< inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_M; ///< Derivative of delta* residual with respect to the Mach
+ ///< number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_v; ///< Derivative of delta* residual with respect to the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_rho; ///< Derivative of delta* residual with respect to the
+ ///< density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_dStar; ///< Derivative of delta* residual with respect to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_blw; ///< Derivative of delta* residual with respect to the
+ ///< blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRdStar_x; ///< Derivative of delta* residual with respect to the mesh
+ ///< coordinates
+
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_phi; ///< Derivative of the blowing velocity residual with respect
+ ///< to the inviscid flow
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_M; ///< Derivative of the blowing velocity residual with respect to
+ ///< the Mach number
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_v; ///< Derivative of the blowing velocity residual with respect to
+ ///< the velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_rho; ///< Derivative of the blowing velocity residual with respect
+ ///< to the density
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_dStar; ///< Derivative of the blowing velocity residual with respect
+ ///< to delta*
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_blw; ///< Derivative of the blowing velocity residual with respect
+ ///< to the blowing velocity
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRblw_x; ///< Derivative of the blowing velocity residual with respect to
+ ///< the mesh coordinates
// Objective functions
- 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
- Eigen::VectorXd dCd_M; ///< Derivative of the drag coefficient with respect to
- ///< the Mach number
- Eigen::VectorXd dCd_rho; ///< Derivative of the drag coefficient with respect
- ///< to the density
- Eigen::VectorXd dCd_v; ///< Derivative of the drag coefficient with respect to
- ///< the velocity
- Eigen::VectorXd
- dCd_dStar; ///< Derivative of the drag coefficient with respect to delta*
+ 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
+ Eigen::VectorXd dCd_M; ///< Derivative of the drag coefficient with respect to
+ ///< the Mach number
+ Eigen::VectorXd dCd_rho; ///< Derivative of the drag coefficient with respect
+ ///< to the density
+ Eigen::VectorXd dCd_v; ///< Derivative of the drag coefficient with respect to
+ ///< the velocity
+ Eigen::VectorXd dCd_dStar; ///< Derivative of the drag coefficient with respect to delta*
// Angle of attack
Eigen::VectorXd dRphi_AoA; ///< Derivative of the inviscid flow residual with
@@ -207,44 +165,31 @@ private:
///< respect to the angle of attack
// Mesh
- Eigen::SparseMatrix<double, Eigen::RowMajor>
- dRx_x; ///< Derivative of the mesh residual with respect to the mesh
- ///< coordinates
- Eigen::VectorXd dCl_x; ///< Partial derivative of the lift coefficient with
- ///< respect to the mesh coordinates
- Eigen::VectorXd dCdp_x; ///< Partial derivative of the pressure drag
- ///< coefficient with respect to the mesh coordinates
- Eigen::VectorXd dCdf_x; ///< Partial derivative of the friction drag
- ///< coefficient with respect to the mesh coordinates
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_x; ///< Derivative of the mesh residual with respect to the mesh
+ ///< coordinates
+ Eigen::VectorXd dCl_x; ///< Partial derivative of the lift coefficient with
+ ///< respect to the mesh coordinates
+ Eigen::VectorXd dCdp_x; ///< Partial derivative of the pressure drag
+ ///< coefficient with respect to the mesh coordinates
+ Eigen::VectorXd dCdf_x; ///< Partial derivative of the friction drag
+ ///< coefficient with respect to the mesh coordinates
fwk::Timers tms; //< internal timers
fwk::Timers tms2; //< internal precise timers
public:
CoupledAdjoint(std::shared_ptr<dart::Adjoint> iAdjoint,
- std::shared_ptr<blast::Driver> vSolver);
+ std::shared_ptr<blast::blAdjoint> vAdjoint);
virtual ~CoupledAdjoint() { std::cout << "~blast::CoupledAdjoint()\n"; }
void run();
void reset();
private:
- void buildAdjointMatrix(
- std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>>
- &matrices,
- Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
-
- Eigen::VectorXd runViscous();
-
- void gradientswrtInviscidFlow();
- void gradientswrtMachNumber();
- void gradientswrtDensity();
- void gradientswrtVelocity();
- void gradientswrtDeltaStar();
- void gradientswrtBlowingVelocity();
-
- void gradientwrtAoA();
- void gradientwrtMesh();
+ void buildAdjointMatrix(std::vector<std::vector<Eigen::SparseMatrix<double, Eigen::RowMajor> *>> &matrices, Eigen::SparseMatrix<double, Eigen::RowMajor> &matrixAdjoint);
+
+ void gradientsInviscid();
+ void gradientsViscous();
void transposeMatrices();
void writeMatrixToFile(const Eigen::MatrixXd &matrix,
diff --git a/blast/src/blDriver.h b/blast/src/blDriver.h
index 79356e9..61dffa7 100644
--- a/blast/src/blDriver.h
+++ b/blast/src/blDriver.h
@@ -19,21 +19,23 @@ namespace blast
*/
class BLAST_API Driver : public fwk::wSharedObject
{
+
+ friend class blAdjoint;
+
public:
std::vector<std::shared_ptr<Body>> bodies; ///< Boundary layer regions sorted by body and section.
unsigned int verbose; ///< Verbosity level of the solver 0<=verbose<=1.
private:
- bldouble Re; ///< Freestream Reynolds number.
- bldouble CFL0; ///< Initial CFL number for pseudo time integration.
- Solver *tSolver; ///< Pseudo-time solver.
- int solverExitCode; ///< Exit code of viscous calculations.
+ bldouble Re; ///< Freestream Reynolds number.
+ bldouble CFL0; ///< Initial CFL number for pseudo time integration.
+ Solver *tSolver; ///< Pseudo-time solver.
+ int solverExitCode; ///< Exit code of viscous calculations.
+protected:
bldouble Cdt;
bldouble Cdf;
bldouble Cdp;
-
-protected:
fwk::Timers tms; ///< Internal timers
public:
diff --git a/blast/src/blNode.cpp b/blast/src/blNode.cpp
index 2628aa7..57d0c18 100644
--- a/blast/src/blNode.cpp
+++ b/blast/src/blNode.cpp
@@ -27,7 +27,7 @@ void blNode::write(std::ostream &out) const
out << "node #" << no << " = " << pos.transpose() << " row=" << row;
}
-blNode::blNode(int n, blVector3d const &p, int r) : wObject(), no(n), pos(p), row(r)
+blNode::blNode(int n, blVector3d const &pos_, int row_, int adjrow_) : no(n), pos(pos_), row(row_), adjrow(adjrow_)
{
xoc = 0.0;
xi = 0.0;
diff --git a/blast/src/blNode.h b/blast/src/blNode.h
index 84d1604..9a777f0 100644
--- a/blast/src/blNode.h
+++ b/blast/src/blNode.h
@@ -31,12 +31,13 @@ namespace blast
/**
* @brief Boundary layer node data
*/
-class BLAST_API blNode : public fwk::wObject
+class BLAST_API blNode : public fwk::wSharedObject
{
public:
int no;
blVector3d pos;
int row;
+ int adjrow;
bldouble xoc; ///< x/c coordinate in the global frame of reference.
bldouble xi; ///< coordinate in the local frame of reference.
@@ -46,8 +47,8 @@ private:
int regime; ///< Regime of the boundary layer (0=laminar, 1=turbulent).
public:
- blNode(int n, blVector3d const &p, int r);
- void setRegime(int const r) { regime = r; }
+ blNode(int n, blVector3d const &pos_, int row_, int adjrow_);
+ void setRegime(int const reg) { regime = reg; }
int getRegime() const { return regime; }
#ifndef SWIG
diff --git a/blast/src/blast.h b/blast/src/blast.h
index 239c65f..c7f421b 100644
--- a/blast/src/blast.h
+++ b/blast/src/blast.h
@@ -32,6 +32,7 @@
#include "tbox.h"
#include "blaster_def.h"
#include <Eigen/Dense>
+#include <Eigen/Sparse>
#ifndef SWIG
#ifndef USE_CODI
@@ -43,6 +44,7 @@ using bltape = bldouble::Tape;
#endif // USE_CODI
using blMatrixXd = Eigen::Matrix<bldouble, Eigen::Dynamic, Eigen::Dynamic>;
+using blSparseMatrix = Eigen::SparseMatrix<bldouble, Eigen::RowMajor>;
using blVectorXd = Eigen::Matrix<bldouble, Eigen::Dynamic, 1>;
using blMatrix3d = Eigen::Matrix<bldouble, 3, 3>;
using blVector3d = Eigen::Matrix<bldouble, 3, 1>;
@@ -59,6 +61,9 @@ class Driver;
class Solver;
// Adjoint
+#ifdef USE_CODI
+class blAdjoint;
+#endif // USE_CODI
class CoupledAdjoint;
// Data structures
diff --git a/blast/tests/dart/t_adjoint_2D.py b/blast/tests/dart/t_adjoint_2D.py
index d04cf01..c205189 100644
--- a/blast/tests/dart/t_adjoint_2D.py
+++ b/blast/tests/dart/t_adjoint_2D.py
@@ -104,9 +104,10 @@ def main():
coupler, isol, vsol = vutils.initBlast(icfg, vcfg, task='optimization')
morpher = isol.iobj['mrf']
-
+
# Adjoint objects
- cAdj = blast.CoupledAdjoint(isol.adjointSolver, vsol)
+ vadj = blast.blAdjoint(vsol, isol.getnDim())
+ cAdj = blast.CoupledAdjoint(isol.adjointSolver, vadj)
# Run coupled case
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
@@ -137,7 +138,7 @@ def main():
coupler.reset()
aeroCoeffs = coupler.run(write=False)
dcl.append(isol.getCl())
- dcd.append(isol.getCd()+vsol.Cdf)
+ dcd.append(isol.getCd()+vsol.getCdf())
dCl_AoA_FD = (dcl[-1] - dcl[-2])/(2.*dalpha)
dCd_AoA_FD = (dcd[-1] - dcd[-2])/(2.*dalpha)
tms['fd_blaster_aoa'].stop()
@@ -187,7 +188,7 @@ def main():
if cmp % 20 == 0:
print('F-D opers', cmp, '/', isol.blw['airfoil']['wing'].nodes.size() * isol.getnDim() * 2, '(node '+str(n.row)+')', 'time', str(round(time.time() - start_time, 3))+'s')
cl[isgn] = isol.getCl()
- cd[isgn] = isol.getCd() + vsol.Cdf
+ cd[isgn] = isol.getCd() + vsol.getCdf()
# Reset mesh
for node in isol.iobj['sol'].pbl.msh.nodes:
--
GitLab