Add tbb loop in Adjoint. Add doc. Refactoring.

dart/api/core.py

@@ -219,12 +219,12 @@ def init_dart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
         linsol = tbox.Pardiso()
     elif cfg['LSolver'] == 'MUMPS':
         linsol = tbox.Mumps()
+    elif cfg['LSolver'] == 'SparseLU':
+        linsol = tbox.SparseLu()
     elif cfg['LSolver'] == 'GMRES':
         linsol = tbox.Gmres(cfg.get('G_fill', 2), 1e-6, cfg.get('G_restart', 50), cfg.get('G_tol', 1e-5), 200)
-    elif cfg['LSolver'] == 'SparseLu':
-        linsol = tbox.SparseLu()
     else:
-        raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
+        raise RuntimeError('Available linear solvers: PARDISO, MUMPS, SparseLU or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
     # initialize the nonlinear (outer) solver
     _sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb)
 

dart/src/wAdjoint.cpp

@@ -564,25 +564,24 @@ 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 blow : sol->pbl->bBCs)
     {
-        for (auto b : sol->pbl->bBCs)
-        {
-            for (auto p : b->uE)
+        tbb::parallel_for_each(blow->uE.begin(), blow->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
+            // Build elementary matrices
+            Eigen::MatrixXd A = BlowingResidual::buildGradientMesh(*p.first, sol->phi, *p.second, sol->pbl->nDim);
+            // Assembly
+            tbb::spin_mutex::scoped_lock lock(mutex);
+            for (size_t i = 0; i < p.first->nodes.size(); ++i)
             {
-                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)
                 {
-                    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)));
-                    }
+                    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

dart/src/wAdjoint.h

@@ -87,12 +87,13 @@ public:
     std::vector<double> computeGradientCoefficientsAoa();
     std::vector<std::vector<double>> computeGradientCoefficientsMesh();
 
+#ifndef SWIG
+    // Getters for partials (required by VII-BLASTER)
     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

@@ -59,25 +59,6 @@ void Blowing::create()
     }
 }
 
-/**
- * @brief Set blowing velocity
- */
-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

@@ -41,8 +41,8 @@ public:
     Blowing(std::shared_ptr<tbox::MshData> _msh, std::string const &names);
     virtual ~Blowing();
 
-    void setU(size_t i, double ue);
-    double getU(size_t i) const;
+    void setU(size_t i, double ue) { uE[i].second->update(ue); }
+    double getU(size_t i) const { return uE[i].second->eval(*uE[i].first, std::vector<double>(), 0); }
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;

dart/src/wBlowingResidual.cpp

@@ -27,7 +27,7 @@ using namespace dart;
 /**
  * @brief Build the residual vector, on one boundary element
  *
- * b = \int psi * vn dV
+ * b = \int psi * vn dS
  */
 Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, F0El const &f)
 {
@@ -47,11 +47,16 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
 
 /**
  * @brief Build the gradient of the residual vector with respect to the blowing velocity, on one blowing boundary element
+ *
+ * b = \int psi dS
  */
 Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
 {
+    // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
+
+    // Build velocity gradient
     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);
@@ -60,19 +65,23 @@ Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
 
 /**
  * @brief Build the gradient of the residual vector with respect to the mesh, on one blowing boundary element
+ *
+ * A = \int psi * vn ddS
  */
 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());
+    // Get pre-computed values
     Cache &cache = e.getVCache();
     Gauss &gauss = cache.getVGauss();
-    // Blowing velocity
+
+    // Build velocity gradient
     double vn = f.eval(e, phi, 0);
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
     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];
+        std::vector<double> dDetJ = e.getGradDetJ(k);
+        for (size_t i = 0; i < dDetJ.size(); i++)
+            A.col(i) += cache.getSf(k) * vn * gauss.getW(k) * dDetJ[i];
     }
     return A;
 }

dart/src/wSolver.h

@@ -80,7 +80,7 @@ public:
     virtual Status run();
     void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
 
-    int getRow(int i) const { return rows[i]; };
+    int getRow(size_t i) const { return rows[i]; };
 
 protected:
     void computeFlow();