Feature eigen

Context

This MR generalizes the use of Eigen in waves.

Tasks

Following MR !50 (merged), here is the updated to-do list:

• replace gmm::matrix<double> by Eigen::Matrix<double>
• replace std::vector<double> by Eigen::(Row)Vector<double>
• refactor elementary matrices computation
• change method return type, ie void build(Matrix const &A, Matrix &K) should become Matrix K build(Matrix const &A)
• replace tbox::Pt by Eigen::VectorXd
• investigate Pardiso interface
• investigate Eigen multithreading feature
• clean the thermo-mechanical matrices

Notes

• I added a general interface tbox::LinearSolver from which derive tbox::Mumps, tbox::Pardiso, tbox::Gmres and tbox::SparseLu. The first class is a direct interface to the MUMPS solver, while the other three merely call Eigen interface or implementation of these solvers. The linear solvers are interfaced in python and the nonlinear solvers must now be initialized with a problem and a linear solver.
• I did not clean the thermo-mechanical matrices computation in tbox::Tetra4 and tbox::Hex8. This could be further addressed in the branch clean_mirrors, which brings LFS support for large geometry files and improved tests.
• Regarding the multi-threading in Eigen, it does not make sense to use it if the top level functions in waves are multi-threaded, which should be the case. However, multi-threading can be used for linear solvers. For now, only Pardiso can use multi-threading (GMRES and SparseLU are not multi-threaded and I did not look into MUMPS).

Additional changes

• I replaced the std::vector holding the shape functions by Eigen::Matrices, and consequently updated the Cache, ShapeFunction and Gauss classes. Additionally, the members of Cache and Gauss are now private and can only be accessed by getters (which I inlined).
• I replaced the std::vector by std::deque to hold the list of Eigen::Triplet for the global matrices assembly, as it proved to be faster.
• I changed the assembly procedure of periodic BCs (tbox::MshDeform) and Kutta condition (flow::Newton and flow::Picard) so that the contributions are directly assembled on the right rows, thus avoiding the costly manipulation of the rows after the matrix has been built.

Tests

Passed on ubuntu18.04 (python 2.7.15), debian4.9 (python 2.7.13), msys2 (python 3.8.1).

flow/src/wNewton.cpp

@@ -266,20 +266,6 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
     // Multithread
     tbb::spin_mutex mutex;
 
-    // List of rows to build Jacobian matrix
-    std::vector<int> rows(pbl->msh->nodes.size());
-    for (size_t i = 0; i < rows.size(); ++i)
-        rows[i] = pbl->msh->nodes[i]->row;
-    // Kutta condition step 1: equality of normal flux through the wake
-    // -> add the upper wake node rows to lower wake nodes rows and zero the upper wake node rows
-    // Nishida Thesis, pp. 29, Eq 2.17
-    // -> directly assemble upper rows on lower rows
-    // Galbraith paper
-    for (auto wake : pbl->wBCs)
-        tbb::parallel_do(wake->nodMap.begin(), wake->nodMap.end(), [&](std::pair<Node *, Node *> p) {
-            rows[p.first->row] = p.second->row;
-        });
-
     // List of triplets to build Jacobian matrix
     std::deque<Eigen::Triplet<double>> T;
 
@@ -470,7 +456,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
             for (size_t ii = 0; ii < e->nodes.size(); ++ii)
             {
                 Node *nodi = e->nodes[ii];
-                R(nodi->row) -= beU(ii);
+                R(rows[nodi->row]) -= beU(ii);
             }
         }
         // Assembly (subsonic)
@@ -478,7 +464,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
         for (size_t ii = 0; ii < e->nodes.size(); ++ii)
         {
             Node *nodi = e->nodes[ii];
-            R(nodi->row) -= be(ii);
+            R(rows[nodi->row]) -= be(ii);
         }
     });
     // Apply freestream (Neumann) BCs
@@ -493,7 +479,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
             for (size_t ii = 0; ii < e->nodes.size(); ++ii)
             {
                 Node *nodi = e->nodes[ii];
-                R(nodi->row) += be(ii);
+                R(rows[nodi->row]) += be(ii);
             }
         });
     }
@@ -509,25 +495,13 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
             for (size_t ii = 0; ii < p.first->nodes.size(); ++ii)
             {
                 Node *nodi = p.first->nodes[ii];
-                R(nodi->row) += be(ii);
+                R(rows[nodi->row]) += be(ii);
             }
         });
     }
     // Slip BCs are naturally enforced in FEM, hence not applied explicitely
 
-    // Kutta condition (3 steps)
-    // Equality of normal flux through the wake
-    // -> add the upper wake node rows to lower wake nodes rows and zero the upper wake node rows
-    // Nishida Thesis, pp. 29, Eq 2.17
-    // here we do not use Galbraith method (as in J) as it would be slower (rows[] needs to be rebuilt!)
-    for (auto wake : pbl->wBCs)
-    {
-        tbb::parallel_do(wake->nodMap.begin(), wake->nodMap.end(), [&](std::pair<Node *, Node *> p) {
-            R(p.second->row) += R[p.first->row];
-            R(p.first->row) = 0.;
-        });
-    }
-    // Equality of velocity norm through the wake
+    // Kutta condition step 2: equality of velocity norm through the wake
     // -> add the Kutta equation on the upper wake node rows
     // Nishida Thesis, pp. 29-30, Eq. 2.18
     for (auto wake : pbl->wBCs)
@@ -546,7 +520,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
             }
         });
     }
-    // Equality of velocity norm on the Trailing edge (Kutta)
+    // Kutta condition step 3: equality of velocity norm on the Trailing edge (Kutta)
     // -> add the Kutta equation on the upper TE node rows (only works for 2D cases)
     // Galbraith paper, pp. 7, Eq. 34
     for (auto kutta : pbl->kSCs)
@@ -567,10 +541,8 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
 
     // Apply Dirichlet BCs
     for (auto dBC : pbl->dBCs)
-    {
         for (auto nod : dBC->nodes)
             R(nod->row) = 0.;
-    }
 
     if (verbose > 1)
         std::cout << "R (" << R.size() << ")\n";