Update waves

Created by: acrovato

Context

This pull request brings the last developments from the flow module as well as some upgrade to waves.

Major changes in flow

• flow is now able to compute 3D flows
• The Jacobian (tangent) matrix for the Netwton method is now fully analytical
• The adjacent elements identification is now based on ```std::unordered_set````

Major changes in waves

• tbox and fwk are now included in a namespace
• waves is now fully included in a namespace
• All elements now have a private memory (Mem) which can be publicly accessed
• All elements now have a cache including Gauss points (Gauss) which can be publicly accessed
• MshDeformation now uses a linear elasticity solver
• MshCrack class implemented (creates duplicate internal boundaries in meshes)

Test

The ctest battery has been run on gaston, 100% tests passed under 110s.

Future work

• Use automated indentation rules for the whole code
• Move methods using VTK outside of MshData
• Clean tbox (mainly MshData, Tri3, Tera4 and Hex8)

flow/src/wSolver.cpp

@@ -14,6 +14,7 @@
 #include "wElement.h"
 #include "wLinesearch.h"
 #include "wMem.h"
+#include "wCache.h"
 
 #include <gmm/gmm.h>
 #include <gmm/gmm_MUMPS_interface.h>
@@ -174,8 +175,10 @@ bool Solver::run()
         }
 
         // Build Jacobian matrix
+        tms["Jac"].start();
         gmm::csr_matrix<double> J;
         buildJac(J);
+        tms["Jac"].stop();
 
         // Solve linear SoE J \DeltaPhi = R
         //gmm::iteration iter(1e-8);
@@ -226,6 +229,8 @@ bool Solver::run()
         else continue;
     } while(itC < maxIt);
 
+    std::cout << "TMS" << tms << std::endl;
+
     // Compute field variables
     computeFlow();
 
@@ -288,6 +293,85 @@ void Solver::buildJac(gmm::csr_matrix<double> &J)
 	        }
         }
         else {
+            /* Analytical NR try */
+            /* @todo clean, make mu function, change to new buildNK and buildNs, rework grad2, rm timers */
+
+            // switching function
+            double mu = std::max(muCv*(1-(mCOv*mCOv)/(mach*mach)) ,0.);
+            // to merge with upper part
+            gmm::dense_matrix<double> Ae(e->nodes.size(), e->nodes.size());
+            e->buildNK(Ae, phi, fluid->rho);
+            // scale - term 1
+            gmm::scale(Ae, 1-mu);
+
+            // term 3 (mu*rhoU*dkNj*dkNi)
+            gmm::dense_matrix<double> Ae_(e->nodes.size(), e->nodes.size());
+            e->buildK(Ae_, phi, Fct0C(1.));
+            gmm::scale(Ae_, mu*fluid->rho.eval(*eU, phi, 0));
+
+            // term 2 (mu*drhoU*dkNi)
+            std::vector<double> bei(e->nodes.size());
+            std::vector<double> bejU(eU->nodes.size());
+            e->buildNs(bei, phi, Fct0C(1.));
+            // warning
+            std::vector<double> gradU(pbl->nDim);
+            eU->computeGrad(phi, gradU, 0);
+            std::vector<double> dNU(pbl->nDim);
+            Cache &cacheU = eU->getVCache();
+            Mem &memU = eU->getVMem();
+            for (size_t j = 0; j < eU->nodes.size(); ++j) {
+                gmm::mult(memU.getJinv(0), gmm::col_vector(cacheU.dff[0][j]), gmm::col_vector(dNU));
+                for (size_t i = 0; i < pbl->nDim; ++i)
+                    bejU[j] += gradU[i]*dNU[i];
+            }
+            gmm::dense_matrix<double> Ae__(e->nodes.size(), eU->nodes.size());
+            for (size_t i = 0; i < e->nodes.size(); ++i)
+                for (size_t j = 0; j < eU->nodes.size(); ++j)
+                    Ae__(i,j) = mu*fluid->rho.evalD(*eU, phi, 0) * bei[i]*bejU[j];
+
+            // term 4 -(rho-rhoU)*dmu*dkNi
+            // make function mu and use evalD...
+            std::vector<double> bej(e->nodes.size());
+            std::vector<double> grad(pbl->nDim);
+            e->computeGrad(phi, grad, 0);
+            std::vector<double> dN(pbl->nDim);
+            Cache &cache = e->getVCache();
+            Mem &mem = e->getVMem();
+            for (size_t j = 0; j < e->nodes.size(); ++j) {
+                gmm::mult(mem.getJinv(0), gmm::col_vector(cache.dff[0][j]), gmm::col_vector(dN));
+                for (size_t i = 0; i < pbl->nDim; ++i)
+                    bej[j] += grad[i]*dN[i];
+            }
+            double grad2 = 0;
+            e->computeGrad2(phi, grad2, 0);
+            double a2 = 1/(0.7*0.7) + 0.5*(1.4-1) * (1 - grad2);
+            double dm = 2*muCv*mCOv*mCOv/pow(mach,3.)*(1/sqrt(grad2*a2)+ 0.5*(1.4-1)*sqrt(grad2)/pow(a2, 3./2.));
+            double deltaRho = (fluid->rho.eval(*eU, phi, 0)-fluid->rho.eval(*e, phi, 0));
+            gmm::dense_matrix<double> Ae___(e->nodes.size(), e->nodes.size());
+            for (size_t i = 0; i < e->nodes.size(); ++i)
+                for (size_t j = 0; j < e->nodes.size(); ++j)
+                    Ae___(i,j) = deltaRho*dm * bei[i]*bej[j];
+
+
+            // Assembly
+            tbb::spin_mutex::scoped_lock lock(mutex);
+	        for(size_t ii=0; ii<e->nodes.size(); ++ii)
+	        {
+	        	Node *nodi = e->nodes[ii];
+	        	for(size_t jj=0; jj<e->nodes.size(); ++jj)
+	        	{
+	        		Node *nodj = e->nodes[jj];
+	        		A(nodi->row,nodj->row) += Ae(ii,jj) + Ae_(ii,jj) + Ae___(ii,jj);
+	        	}
+                for(size_t jj=0; jj<eU->nodes.size(); ++jj)
+	        	{
+                    Node *nodj = eU->nodes[jj];
+                    A(nodi->row,nodj->row) += Ae__(ii,jj);
+                }
+	        }
+
+            /* @todo numerical Jac*/
+/*
             // Merge nodes of current and upwind elements (and remove duplicate)
             std::vector<Node *> nodeU;
             for (auto n : e->nodes)
@@ -335,6 +419,7 @@ void Solver::buildJac(gmm::csr_matrix<double> &J)
                 phiPlu[nodj->row] = phiTmp;
                 phiMin[nodj->row] = phiTmp;
             }
+*/
         }
 	});