dart/api/core.py

@@ -249,6 +249,8 @@ def init_dart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
         grst = cfg['G_restart'] if 'G_restart' in cfg else 50
         gtol = cfg['G_tol'] if 'G_tol' in cfg else 1e-5
         linsol = tbox.Gmres(gfil, 1e-6, grst, gtol, 200)
+    elif cfg['LSolver'] == 'SparseLu':
+        linsol = tbox.SparseLu()
     else:
         raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
     # initialize the nonlinear (outer) solver

dart/src/wAdjoint.cpp

@@ -256,12 +256,6 @@ void Adjoint::linearizeFlow()
     // Partials wrt AoA
     dRu_A = Eigen::VectorXd::Zero(sol->pbl->msh->nodes.size());
     this->buildGradientResidualAoa(dRu_A);
-    // Partials wrt Blw
-    if (sol->pbl->bBCs.size() > 0)
-    {
-        dRu_Blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(sol->pbl->msh->nodes.size(), sol->pbl->bBCs[0]->uE.size());
-        this->buildGradientResidualBlowing(dRu_Blw);
-    }
     // Partials wrt mesh
     dRu_X = Eigen::SparseMatrix<double, Eigen::RowMajor>(sol->pbl->msh->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
     this->buildGradientResidualMesh(dRu_X);
@@ -316,28 +310,6 @@ std::vector<double> Adjoint::computeJacVecFlowMesh(std::vector<double> const &dR
     return std::vector<double>(dX.data(), dX.data() + dX.size());
 }
 
-/**
- * @brief Compute the matrix-vector product of the gradients of the flow residuals (wrt blowing velocity) with the flow residuals
- *
- * dX = dRu/dX^T * dRu
- */
-std::vector<double> Adjoint::computeJacVecFlowBlowing(std::vector<double> const &dR)
-{
-    Eigen::VectorXd dX = dRu_Blw.transpose() * Eigen::Map<const Eigen::VectorXd>(dR.data(), dR.size());
-    return std::vector<double>(dX.data(), dX.data() + dX.size());
-}
-
-/**
- * @brief Compute the matrix-matrix product of the gradients of the flow residuals (wrt blowing velocity) with the flow residuals
- *
- * dX = dRu/dX^T * dRu
- */
-Eigen::SparseMatrix<double, Eigen::RowMajor> Adjoint::computeJacVecFlowBlowing(Eigen::SparseMatrix<double, Eigen::RowMajor> const &dR)
-{
-    Eigen::SparseMatrix<double, Eigen::RowMajor> dX = dRu_Blw.transpose() * dR;
-    return dX;
-}
-
 /**
  * @brief Compute the partial gradients of the loads, on a given body
  *
@@ -523,34 +495,6 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
     }
 }
 
-/**
- * @brief Build gradients of the residual with respect to blowing velocity
- */
-void Adjoint::buildGradientResidualBlowing(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
-{
-    // Multithread
-    tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
-    tbb::spin_mutex mutex;
-
-    std::deque<Eigen::Triplet<double>> T;
-
-    size_t jj = 0;
-    tbb::parallel_for_each(sol->pbl->bBCs[0]->uE.begin(), sol->pbl->bBCs[0]->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
-        Eigen::VectorXd be = BlowingResidual::buildGradientBlowing(*p.first);
-
-        tbb::spin_mutex::scoped_lock lock(mutex);
-        for (size_t ii = 0; ii < p.first->nodes.size(); ii++)
-        {
-            Node *nodi = p.first->nodes[ii];
-            T.push_back(Eigen::Triplet<double>(sol->rows[nodi->row], jj, be(ii)));
-        }
-        ++jj;
-    });
-    dR.setFromTriplets(T.begin(), T.end());
-    dR.prune(0.);
-    dR.makeCompressed();
-}
-
 /**
  * @brief Build gradients of the lift and drag with respect to angle of attack
  */
@@ -622,17 +566,20 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
     // Blowing boundary
     if (sol->pbl->bBCs.size() > 0)
     {
-        for (auto p : sol->pbl->bBCs[0]->uE)
+        for (auto b : sol->pbl->bBCs)
         {
-            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)
+            for (auto p : b->uE)
             {
-                size_t rowi = p.first->nodes[i]->row;
-                for (size_t j = 0; j < p.first->nodes.size(); ++j)
+                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)
                 {
-                    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)));
+                    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)));
+                    }
                 }
             }
         }

dart/src/wAdjoint.h

@@ -87,13 +87,9 @@ public:
     std::vector<std::vector<double>> computeGradientCoefficientsFlow();
     std::vector<double> computeGradientCoefficientsAoa();
     std::vector<std::vector<double>> computeGradientCoefficientsMesh();
-    // Gradient of the blowing velocity
-    Eigen::SparseMatrix<double, Eigen::RowMajor> computeJacVecFlowBlowing(Eigen::SparseMatrix<double, Eigen::RowMajor> const &dR);
-    std::vector<double> computeJacVecFlowBlowing(std::vector<double> const &dR);
 
     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_Blw() const { return dRu_Blw; }
     Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_X() const { return dRu_X; }
     Eigen::SparseMatrix<double, Eigen::RowMajor> getRx_X() const { return dRx_X; }
 

dart/src/wSolver.h

@@ -67,10 +67,10 @@ public:
     double Cd;                      ///< drag coefficient
     double Cs;                      ///< sideforce coefficient
     double Cm;                      ///< pitch moment coefficient (positive nose-up)
+    std::vector<int> rows;          ///< unknown nodal index
 
 protected:
     fwk::Timers tms;       ///< internal timers
-    std::vector<int> rows; ///< unknown nodal index
 
 public:
     Solver(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol, double rTol, double aTol, int mIt, int nthrds, int vrb);