(feat) Added adjoint for viscous cases

Computation of blowing velocity residuals gradients and blowing contribution to mesh gradients

dart/src/wAdjoint.cpp

@@ -27,6 +27,8 @@
 #include "wWakeResidual.h"
 #include "wKuttaResidual.h"
 #include "wLoadFunctional.h"
+#include "wBlowing.h"
+#include "wBlowingResidual.h"
 
 #include "wMshData.h"
 #include "wNode.h"
@@ -561,6 +563,28 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
     });
     // Farfield B.C. are not taken into account because they have no influence on the solution
 
+    // Blowing boundary
+    if (sol->pbl->bBCs.size() > 0)
+    {
+        for (auto b : sol->pbl->bBCs)
+        {
+            for (auto p : b->uE)
+            {
+                Eigen::MatrixXd A = BlowingResidual::buildGradientMesh(*p.first, sol->phi, *p.second, sol->pbl->nDim);
+                for (size_t i = 0; i < p.first->nodes.size(); ++i)
+                {
+                    size_t rowi = p.first->nodes[i]->row;
+                    for (size_t j = 0; j < p.first->nodes.size(); ++j)
+                    {
+                        int rowj = p.first->nodes[j]->row;
+                        for (int n = 0; n < sol->pbl->nDim; ++n)
+                            T.push_back(Eigen::Triplet<double>(sol->rows[rowi], sol->pbl->nDim * rowj + n, A(i, sol->pbl->nDim * j + n)));
+                    }
+                }
+            }
+        }
+    }
+
     // Wake
     for (auto wake : sol->pbl->wBCs)
     {

dart/src/wAdjoint.h

@@ -87,6 +87,11 @@ public:
     std::vector<double> computeGradientCoefficientsAoa();
     std::vector<std::vector<double>> computeGradientCoefficientsMesh();
 
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_U() const { return dRu_U; }
+    Eigen::VectorXd getRu_A() const { return dRu_A; }
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_X() const { return dRu_X; }
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRx_X() const { return dRx_X; }
+
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
 #endif

dart/src/wBlowing.cpp

@@ -67,6 +67,17 @@ void Blowing::setU(size_t i, double ue)
     uE[i].second->update(ue);
 }
 
+/**
+ * @brief Get blowing velocity
+ */
+double Blowing::getU(size_t i) const
+{
+    // Return a constant value
+    tbox::Element *e = uE[i].first;
+    std::vector<double> phi(0, 0.);
+    return uE[i].second->eval(*e, phi, 0);
+}
+
 void Blowing::write(std::ostream &out) const
 {
     out << " Blowing boundary condition on " << *tag << std::endl;

dart/src/wBlowing.h

@@ -42,6 +42,7 @@ public:
     virtual ~Blowing();
 
     void setU(size_t i, double ue);
+    double getU(size_t i) const;
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;

dart/src/wBlowingResidual.cpp

@@ -44,3 +44,35 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
         b += cache.getSf(k) * vn * gauss.getW(k) * e.getDetJ(k);
     return b;
 }
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the blowing velocity, on one blowing boundary element
+ */
+Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
+{
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+    Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
+    for (size_t k = 0; k < gauss.getN(); ++k)
+        b += cache.getSf(k) * gauss.getW(k) * e.getDetJ(k);
+    return b;
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the mesh, on one blowing boundary element
+ */
+Eigen::MatrixXd BlowingResidual::buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim)
+{
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+    // Blowing velocity
+    double vn = f.eval(e, phi, 0);
+    for (size_t k = 0; k < gauss.getN(); ++k)
+    {
+        std::vector<double> gradJ = e.getGradDetJ(k);
+        for (size_t i = 0; i < gradJ.size(); i++)
+            A.col(i) += cache.getSf(k) * vn * gauss.getW(k) * gradJ[i];
+    }
+    return A;
+}

dart/src/wBlowingResidual.h

@@ -33,6 +33,8 @@ class DART_API BlowingResidual
 public:
     // Newton
     static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, F0El const &f);
+    static Eigen::VectorXd buildGradientBlowing(tbox::Element const &e);
+    static Eigen::MatrixXd buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim);
 };
 
 } // namespace dart