(feat) Added auto-diff adjoint solver for boundary layer

blast/_src/blastw.h

@@ -21,5 +21,6 @@
 #include "blNode.h"
 
 // Solvers
+#include "blAdjoint.h"
 #include "blCoupledAdjoint.h"
 #include "blDriver.h"

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;

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:

blast/src/blAdjoint.h

+/*
+ * 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

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

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:

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;

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

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

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: