diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 96bc4057ca80d1ef413ea396522dd8179a46ccc9..a200fa2df84bcfeedef912f5b40e5f4e60b37468 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,7 +1,7 @@
# gitlab-ci file for waves
default:
- image: rboman/waves
+ image: rboman/waves:2020.0
before_script:
- source /opt/intel/mkl/bin/mklvars.sh intel64
- source /opt/intel/tbb/bin/tbbvars.sh intel64
@@ -35,7 +35,7 @@ build-py2-no-tlnos:
- make -j 4
build-py3:
- image: rboman/waves-py3
+ image: rboman/waves-py3:2020.0
stage: build
script:
- printenv | sort
@@ -71,7 +71,7 @@ ctest-py2:
- build-py2
ctest-py3:
- image: rboman/waves-py3
+ image: rboman/waves-py3:2020.0
stage: test
script:
- cd build
diff --git a/flow/src/wAdjoint.cpp b/flow/src/wAdjoint.cpp
index caf2b63a8663cda55c089db85f51a9ec08fea9bb..cbf88d7618f37d4383ca514bdf5681049a6281c5 100644
--- a/flow/src/wAdjoint.cpp
+++ b/flow/src/wAdjoint.cpp
@@ -32,11 +32,9 @@
#include "wMem.h"
#include "wResults.h"
#include "wMshExport.h"
+#include "wMumpsInterface.h"
-#include "gmm_gmm.h"
-#include <gmm/gmm_MUMPS_interface.h>
-#include "gmm_extra.h"
-#include "std_extra.h"
+#include <Eigen/Sparse>
#include <tbb/task_scheduler_init.h>
#include <tbb/tbb.h>
@@ -82,26 +80,27 @@ void Adjoint::run()
std::cout << "--- Adjoint solver ---\n"
<< "Number of threads: " << nthreads << "\n"
<< std::endl;
+ // Map solution vectors
+ Eigen::Map<Eigen::VectorXd> lambdaL_(lambdaL.data(), lambdaL.size()), lambdaD_(lambdaD.data(), lambdaD.size());
// build SoE
- gmm::csr_matrix<double> dR;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> dR(sol->pbl->msh->nodes.size(), sol->pbl->msh->nodes.size());
buildDR(dR);
- std::vector<double> dL(sol->pbl->msh->nodes.size()), dD(sol->pbl->msh->nodes.size());
- buildDI(dL, dD);
-
- // Solve linear SoE dR {lambdaL, lambdaD} = {dL, dD}
- bool solve_success = gmm::MUMPS_solve(dR, lambdaL, dL);
- if (!solve_success)
- throw std::runtime_error("Problem occured in gmm::MUMPS_solve while solving for lift sensitivities!");
- solve_success = gmm::MUMPS_solve(dR, lambdaD, dD);
- if (!solve_success)
- throw std::runtime_error("Problem occured in gmm::MUMPS_solve while solving for drag sensitivities!");
+ std::vector<double> dummyL(sol->pbl->msh->nodes.size()), dummyD(sol->pbl->msh->nodes.size()); // dummy bs vector
+ Eigen::Map<Eigen::VectorXd> dL_(dummyL.data(), dummyL.size()), dD_(dummyD.data(), dummyD.size());
+ buildDI(dL_, dD_);
+
+ // Setup linear solver and solve linear SoE dR {lambdaL, lambdaD} = {dL, dD}
+ std::map<std::string, int> mumpsOpt;
+ MumpsInterface mumps(mumpsOpt);
+ mumps.setSoE(dR, dL_);
+ mumps.compute(lambdaL_);
+ mumps.setSoE(dR, dD_);
+ mumps.compute(lambdaD_);
+ mumps.finalize();
// Check residual
- std::vector<double> rLambdaL(sol->pbl->msh->nodes.size()), rLambdaD(sol->pbl->msh->nodes.size()), dR_lambdaL(sol->pbl->msh->nodes.size()), dR_lambdaD(sol->pbl->msh->nodes.size());
- gmm::mult(dR, lambdaL, dR_lambdaL);
- gmm::mult(dR, lambdaD, dR_lambdaD);
- gmm::add(dR_lambdaL, gmm::scaled(dL, -1.0), rLambdaL);
- gmm::add(dR_lambdaD, gmm::scaled(dD, -1.0), rLambdaD);
+ Eigen::VectorXd rLambdaL = dR * lambdaL_ - dL_;
+ Eigen::VectorXd rLambdaD = dR * lambdaD_ - dD_;
std::cout << std::setw(8) << "L_Iter"
<< std::setw(15) << "Sens[lambdaL]"
@@ -110,10 +109,10 @@ void Adjoint::run()
<< std::setw(15) << "Res[lambdaD]" << std::endl;
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8) << 0
- << std::setw(15) << gmm::vect_norm2(lambdaL)
- << std::setw(15) << gmm::vect_norm2(lambdaD)
- << std::setw(15) << log10(gmm::vect_norm2(rLambdaL))
- << std::setw(15) << log10(gmm::vect_norm2(rLambdaD)) << std::endl;
+ << std::setw(15) << lambdaL_.norm()
+ << std::setw(15) << lambdaD_.norm()
+ << std::setw(15) << log10(rLambdaL.norm())
+ << std::setw(15) << log10(rLambdaD.norm()) << std::endl;
std::cout << ANSI_COLOR_YELLOW << "Warning: the adjoint solver is experimental and has not been validated!\n"
<< ANSI_COLOR_RESET << std::endl;
@@ -124,13 +123,14 @@ void Adjoint::run()
* @authors Adrien Crovato
* @todo Refactor, this is heavily copy-pasted from Newton::buildJ
*/
-void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
+void Adjoint::buildDR(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
{
// Multithread
tbb::spin_mutex mutex;
- // Temporary Jacobian matrix
- gmm::row_matrix<gmm::wsvector<double>> A(sol->pbl->msh->nodes.size(), sol->pbl->msh->nodes.size());
+ // List of triplets to build Jacobian matrix
+ std::vector<Eigen::Triplet<double>> T;
+ T.reserve(sol->pbl->msh->nodes.size());
// Full Potential Equation with upwind bias: analytical derivatives
auto fluid = sol->pbl->medium;
@@ -160,40 +160,28 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
// 3rd term (mu*drhoU*grad_phi*grad_psi)
Eigen::VectorXd bei;
- std::vector<double> bejU(eU->nodes.size());
+ Eigen::RowVectorXd bejU(eU->nodes.size());
e->buildDs(bei, sol->phi, Fct0C(1.));
std::vector<double> gradU(sol->pbl->nDim);
eU->computeGrad(sol->phi, gradU, 0);
- std::vector<double> dNU(sol->pbl->nDim);
+ Eigen::VectorXd dNU(sol->pbl->nDim);
Cache &cacheU = eU->getVCache();
-
- /* TEMPO EIGEN/GMM */
- Eigen::MatrixXd JJinvU = eU->getVMem().getJinv(0);
- gmm::dense_matrix<double> JinvU(JJinvU.rows(), JJinvU.cols());
- for (size_t i = 0; i < JinvU.nrows(); ++i)
- for (size_t j = 0; j < JinvU.ncols(); ++j)
- JinvU(i, j) = JJinvU(i, j);
-
for (size_t j = 0; j < eU->nodes.size(); ++j)
{
- gmm::mult(JinvU, gmm::col_vector(cacheU.dff[0][j]), gmm::col_vector(dNU));
- for (size_t i = 0; i < sol->pbl->nDim; ++i)
- bejU[j] += gradU[i] * dNU[i];
+ dNU = eU->getVMem().getJinv(0) * Eigen::Map<Eigen::VectorXd>(cacheU.dff[0][j].data(), sol->pbl->nDim);
+ bejU(j) = Eigen::Map<Eigen::VectorXd>(gradU.data(), sol->pbl->nDim).dot(dNU);
}
- gmm::dense_matrix<double> Ae3(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)
- Ae3(i, j) = mu * fluid->rho.evalD(*eU, sol->phi, 0) * bei[i] * bejU[j];
+ Eigen::MatrixXd Ae3(e->nodes.size(), eU->nodes.size());
+ Ae3 = mu * fluid->rho.evalD(*eU, sol->phi, 0) * bei * bejU;
// 4th term (mu*rhoU*grad_phi*grad_psi)
Eigen::MatrixXd Ae4;
e->buildK(Ae4, sol->phi, Fct0C(1.));
Ae4 *= mu * fluid->rho.eval(*eU, sol->phi, 0);
// 5th term dmu*(rhoU-rho)*grad_phi*grad_psi
- double deltaRho = fluid->rho.eval(*eU, sol->phi, 0) - fluid->rho.eval(*e, sol->phi, 0);
Eigen::MatrixXd Ae5;
e->buildDK(Ae5, sol->phi, fluid->mach);
- Ae5 *= deltaRho * dmu;
+ Ae5 *= (fluid->rho.eval(*eU, sol->phi, 0) - fluid->rho.eval(*e, sol->phi, 0)) * dmu;
// Assembly (supersonic)
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -203,12 +191,12 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- A(nodi->row, nodj->row) += Ae4(ii, jj) + Ae5(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae4(ii, jj) + Ae5(ii, jj)));
}
for (size_t jj = 0; jj < eU->nodes.size(); ++jj)
{
Node *nodj = eU->nodes[jj];
- A(nodi->row, nodj->row) += Ae3(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae3(ii, jj)));
}
}
}
@@ -220,10 +208,12 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- A(nodi->row, nodj->row) += Ae1(ii, jj) + Ae2(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae1(ii, jj) + Ae2(ii, jj)));
}
}
});
+ // Build Jacobian matrix without BCs
+ dR.setFromTriplets(T.begin(), T.end());
// Apply wake BCs for adjoint problem
for (auto wake : sol->pbl->wBCs)
{
@@ -241,9 +231,9 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
for (size_t jj = 0; jj < we->nRow; ++jj)
{
Node *nodj = we->surUpE->nodes[jj];
- A(nodi->row, nodj->row) += 2 * Kupup(ii, jj);
+ dR.coeffRef(nodi->row, nodj->row) += 2 * Kupup(ii, jj);
nodj = we->surLwE->nodes[jj];
- A(nodi->row, nodj->row) += Kuplw(ii, jj);
+ dR.coeffRef(nodi->row, nodj->row) += Kuplw(ii, jj);
}
}
for (size_t ii = 0; ii < we->nColLw; ++ii)
@@ -252,9 +242,9 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
for (size_t jj = 0; jj < we->nRow; ++jj)
{
Node *nodj = we->surUpE->nodes[jj];
- A(nodi->row, nodj->row) -= 2 * Klwup(ii, jj);
+ dR.coeffRef(nodi->row, nodj->row) -= 2 * Klwup(ii, jj);
nodj = we->surLwE->nodes[jj];
- A(nodi->row, nodj->row) -= Klwlw(ii, jj);
+ dR.coeffRef(nodi->row, nodj->row) -= Klwlw(ii, jj);
}
}
});
@@ -264,8 +254,13 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
{
for (auto nod : dBC->nodes)
{
- A.row(nod->row) = gmm::wsvector<double>(sol->pbl->msh->nodes.size());
- A(nod->row, nod->row) = 1.0;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(dR, nod->row); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
}
}
for (auto fBC : sol->pbl->fBCs)
@@ -274,17 +269,24 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
{
for (auto nod : e->nodes)
{
- A.row(nod->row) = gmm::wsvector<double>(sol->pbl->msh->nodes.size());
- A(nod->row, nod->row) = 1.0;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(dR, nod->row); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
}
}
}
// Transpose
- gmm::transposed(A);
- gmm::copy(A, dR);
+ dR.transpose();
+ // Clean matrix and turn to compressed row format
+ dR.prune(0.);
+ dR.makeCompressed();
if (verbose)
- std::cout << "dR (" << gmm::mat_nrows(dR) << "," << gmm::mat_ncols(dR) << ") nnz=" << gmm::nnz(dR) << "\n";
+ std::cout << "J (" << dR.rows() << "," << dR.cols() << ") nnz=" << dR.nonZeros() << "\n";
}
/**
@@ -292,7 +294,7 @@ void Adjoint::buildDR(gmm::csr_matrix<double> &dR)
* @authors Adrien Crovato
* @todo Refactor, to be checked
*/
-void Adjoint::buildDI(std::vector<double> &dL, std::vector<double> &dD)
+void Adjoint::buildDI(Eigen::Map<Eigen::VectorXd> &dL, Eigen::Map<Eigen::VectorXd> &dD)
{
// Multithread
tbb::spin_mutex mutex;
@@ -319,8 +321,8 @@ void Adjoint::buildDI(std::vector<double> &dL, std::vector<double> &dD)
for (size_t ii = 0; ii < eV->nodes.size(); ++ii)
{
Node *nodi = eV->nodes[ii];
- dL[nodi->row] += factorL * be(ii);
- dD[nodi->row] += factorD * be(ii);
+ dL(nodi->row) += factorL * be(ii);
+ dD(nodi->row) += factorD * be(ii);
}
});
}
@@ -329,8 +331,8 @@ void Adjoint::buildDI(std::vector<double> &dL, std::vector<double> &dD)
{
for (auto nod : dBC->nodes)
{
- dL[nod->row] = 0.;
- dD[nod->row] = 0.;
+ dL(nod->row) = 0.;
+ dD(nod->row) = 0.;
}
}
for (auto fBC : sol->pbl->fBCs)
@@ -339,16 +341,16 @@ void Adjoint::buildDI(std::vector<double> &dL, std::vector<double> &dD)
{
for (auto nod : e->nodes)
{
- dL[nod->row] = 0.;
- dD[nod->row] = 0.;
+ dL(nod->row) = 0.;
+ dD(nod->row) = 0.;
}
}
}
if (verbose)
{
- std::cout << "dL (" << gmm::vect_size(dL) << ")\n";
- std::cout << "dD (" << gmm::vect_size(dD) << ")\n";
+ std::cout << "dL (" << dL.size() << ")\n";
+ std::cout << "dD (" << dD.size() << ")\n";
}
}
diff --git a/flow/src/wAdjoint.h b/flow/src/wAdjoint.h
index cf752783fd039dd1ad44c2037b4d2cfc390982fd..8cecd60761c272b615e24e9fd9e0f49531c0bf14 100644
--- a/flow/src/wAdjoint.h
+++ b/flow/src/wAdjoint.h
@@ -20,12 +20,10 @@
#include "flow.h"
#include "wObject.h"
-#include "gmm_gmm.h"
-
#include <iostream>
#include <vector>
#include <memory>
-#include <Eigen/Dense>
+#include <Eigen/Sparse>
namespace flow
{
@@ -58,8 +56,8 @@ public:
#endif
private:
- void buildDR(gmm::csr_matrix<double> &dR);
- void buildDI(std::vector<double> &dL, std::vector<double> &dD);
+ void buildDR(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR);
+ void buildDI(Eigen::Map<Eigen::VectorXd> &dL, Eigen::Map<Eigen::VectorXd> &dD);
void buildWake(WakeElement *&we, Eigen::MatrixXd &Kupup, Eigen::MatrixXd &Kuplw, Eigen::MatrixXd &Klwup, Eigen::MatrixXd &Klwlw);
void buildBoundary(tbox::Element *&e, tbox::Element *&eV, Eigen::RowVectorXd &be);
};
diff --git a/flow/src/wAssign.h b/flow/src/wAssign.h
index 4ce9492e7b7c0a774d8c7e04e77b3259cfe3488d..81b8ab0440d96b0627dd8465ff010974000256cd 100644
--- a/flow/src/wAssign.h
+++ b/flow/src/wAssign.h
@@ -22,6 +22,7 @@
#include "wF0Ps.h"
#include <vector>
#include <string>
+#include <Eigen/Dense>
namespace flow
{
diff --git a/flow/src/wFpeLSFunction.cpp b/flow/src/wFpeLSFunction.cpp
index b13cf9e4774bdb9fc79dbdd82f0dcb06706512d7..061ba4c016b814293e19263354ee63126cfc0776 100644
--- a/flow/src/wFpeLSFunction.cpp
+++ b/flow/src/wFpeLSFunction.cpp
@@ -19,9 +19,8 @@ using namespace tbox;
using namespace flow;
FpeLSFunction::FpeLSFunction(flow::Newton &_solver, std::vector<double> &_u,
- std::vector<double> &_u0, std::vector<double> &_du,
- std::vector<double> &_fRes) : LSFunction(),
- solver(_solver), u(_u), u0(_u0), du(_du), fRes(_fRes)
+ std::vector<double> &_u0, Eigen::Map<Eigen::VectorXd> &_du,
+ Eigen::Map<Eigen::VectorXd> &_fRes) : LSFunction(), solver(_solver), u(_u), u0(_u0), du(_du), fRes(_fRes)
{
}
@@ -29,7 +28,7 @@ FpeLSFunction::FpeLSFunction(flow::Newton &_solver, std::vector<double> &_u,
* @brief Update initial and change in solution
* @authors Adrien Crovato
*/
-void FpeLSFunction::update(std::vector<double> &_u0, std::vector<double> &_du)
+void FpeLSFunction::update(std::vector<double> &_u0, Eigen::Map<Eigen::VectorXd> &_du)
{
u0 = _u0;
du = _du;
@@ -41,12 +40,10 @@ void FpeLSFunction::update(std::vector<double> &_u0, std::vector<double> &_du)
*/
double FpeLSFunction::eval(double alpha)
{
- double resNorm = 0.;
- gmm::add(u0, gmm::scaled(du, alpha), u);
- gmm::clear(fRes);
+ for (size_t i = 0; i < u.size(); ++i)
+ u[i] = u0[i] + alpha * du(i);
solver.buildRes(fRes);
- resNorm = gmm::vect_norm2(fRes);
- return resNorm;
+ return fRes.norm();
}
void FpeLSFunction::write(std::ostream &out) const
diff --git a/flow/src/wFpeLSFunction.h b/flow/src/wFpeLSFunction.h
index 4261bd348409c84db9d606980a8c5306e9c06175..4b947e79788da12eeda76c6aa81cdbc1a49ec22e 100644
--- a/flow/src/wFpeLSFunction.h
+++ b/flow/src/wFpeLSFunction.h
@@ -34,15 +34,15 @@ namespace flow
class FLOW_API FpeLSFunction : public tbox::LSFunction
{
private:
- flow::Newton &solver; ///< solver object (must be of Newton type)
- std::vector<double> &u; ///< solution
- std::vector<double> u0; ///< initial solution
- std::vector<double> du; ///< change in solution
- std::vector<double> &fRes; ///< residual vector
+ flow::Newton &solver; ///< solver object (must be of Newton type)
+ std::vector<double> &u; ///< solution
+ std::vector<double> u0; ///< initial solution
+ Eigen::Map<Eigen::VectorXd> du; ///< change in solution
+ Eigen::Map<Eigen::VectorXd> &fRes; ///< residual vector
public:
- FpeLSFunction(flow::Newton &_solver, std::vector<double> &_u, std::vector<double> &_u0, std::vector<double> &_du, std::vector<double> &_fRes);
- void update(std::vector<double> &_u0, std::vector<double> &_du);
+ FpeLSFunction(flow::Newton &_solver, std::vector<double> &_u, std::vector<double> &_u0, Eigen::Map<Eigen::VectorXd> &_du, Eigen::Map<Eigen::VectorXd> &_fRes);
+ void update(std::vector<double> &_u0, Eigen::Map<Eigen::VectorXd> &_du);
virtual double eval(double alpha) override;
#ifndef SWIG
diff --git a/flow/src/wNewton.cpp b/flow/src/wNewton.cpp
index fef2dcb1d832cba156fa6df1dd87d6c0fb3c1a4a..4a77d5516ddaf8d92f648beb449f36fc4944d325 100644
--- a/flow/src/wNewton.cpp
+++ b/flow/src/wNewton.cpp
@@ -31,11 +31,9 @@
#include "wLinesearch.h"
#include "wMem.h"
#include "wCache.h"
+#include "wMumpsInterface.h"
-#include "gmm_gmm.h"
-#include <gmm/gmm_MUMPS_interface.h>
-#include "gmm_extra.h"
-#include "std_extra.h"
+#include <Eigen/Sparse>
#include <tbb/task_scheduler_init.h>
#include <tbb/tbb.h>
@@ -101,19 +99,24 @@ bool Newton::run()
int avC = 0; // adaptive viscosity counter
bool rUpwd = true; // update best upwind element
+ // Setup linear (inner) solver MUMPS
+ std::map<std::string, int> mumpsOpt;
+ MumpsInterface mumps(mumpsOpt);
+
+ // Residual and change in solution
+ std::vector<double> dPhi(phi.size()); // dummy dPhi
+ Eigen::Map<Eigen::VectorXd> rPhi_(rPhi.data(), rPhi.size()), dPhi_(dPhi.data(), dPhi.size());
+
// Find initial best upwind element
findUpwd();
// Compute initial residual
- buildRes(rPhi);
- double resInit = gmm::vect_norm2(rPhi);
- double absRes = gmm::vect_norm2(rPhi);
+ buildRes(rPhi_);
+ double resInit = rPhi_.norm();
+ double absRes = rPhi_.norm();
double relRes = absRes / resInit;
- // Change in solution
- std::vector<double> deltaPhi(phi.size(), 0.);
-
// Linesearch
- FpeLSFunction fpe(*this, phi, phi, deltaPhi, rPhi);
+ FpeLSFunction fpe(*this, phi, phi, dPhi_, rPhi_);
BankRoseLS ls(fpe);
//FleuryLS ls(fpe, 1., 1e-2, false);
ls.set(maxLsIt, lsTol, verbose);
@@ -166,17 +169,16 @@ bool Newton::run()
}
// Build Jacobian matrix
- gmm::csr_matrix<double> J;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> J(pbl->msh->nodes.size(), pbl->msh->nodes.size());
buildJac(J);
// Solve linear SoE J \DeltaPhi = R
- bool solve_success = gmm::MUMPS_solve(J, deltaPhi, rPhi);
- if (!solve_success)
- throw std::runtime_error("Problem occured in gmm::MUMPS_solve!");
+ mumps.setSoE(J, rPhi_);
+ mumps.compute(dPhi_);
itC++;
// Compute new solution and residual with linesearch
- fpe.update(phi, deltaPhi);
+ fpe.update(phi, dPhi_);
try
{
ls.run();
@@ -184,7 +186,7 @@ bool Newton::run()
catch (...)
{
}
- absRes = gmm::vect_norm2(rPhi);
+ absRes = rPhi_.norm();
relRes = absRes / resInit;
// Re-find best upwind element
if (rUpwd)
@@ -215,6 +217,9 @@ bool Newton::run()
continue;
} while (itC < maxIt);
+ // Terminate MUMPS
+ mumps.finalize();
+
// Compute field variables
Solver::computeFlow();
@@ -243,13 +248,14 @@ bool Newton::run()
* @brief Build the Jacobian (tangent) matrix
* @authors Adrien Crovato
*/
-void Newton::buildJac(gmm::csr_matrix<double> &J)
+void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
{
// Multithread
tbb::spin_mutex mutex;
- // Temporary Jacobian matrix
- gmm::row_matrix<gmm::wsvector<double>> A(pbl->msh->nodes.size(), pbl->msh->nodes.size());
+ // List of triplets to build Jacobian matrix
+ std::vector<Eigen::Triplet<double>> T;
+ T.reserve(pbl->msh->nodes.size());
// Full Potential Equation with upwind bias: analytical derivatives
auto fluid = pbl->medium;
@@ -278,30 +284,19 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
// 3rd term (mu*drhoU*grad_phi*grad_psi)
Eigen::VectorXd bei;
- std::vector<double> bejU(eU->nodes.size());
+ Eigen::RowVectorXd bejU(eU->nodes.size());
e->buildDs(bei, phi, Fct0C(1.));
std::vector<double> gradU(pbl->nDim);
eU->computeGrad(phi, gradU, 0);
- std::vector<double> dNU(pbl->nDim);
+ Eigen::VectorXd dNU(pbl->nDim);
Cache &cacheU = eU->getVCache();
-
- /* TEMPO EIGEN/GMM */
- Eigen::MatrixXd JJinvU = eU->getVMem().getJinv(0);
- gmm::dense_matrix<double> JinvU(JJinvU.rows(), JJinvU.cols());
- for (size_t i = 0; i < JinvU.nrows(); ++i)
- for (size_t j = 0; j < JinvU.ncols(); ++j)
- JinvU(i, j) = JJinvU(i, j);
-
for (size_t j = 0; j < eU->nodes.size(); ++j)
{
- gmm::mult(JinvU, 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];
+ dNU = eU->getVMem().getJinv(0) * Eigen::Map<Eigen::VectorXd>(cacheU.dff[0][j].data(), pbl->nDim);
+ bejU(j) = Eigen::Map<Eigen::VectorXd>(gradU.data(), pbl->nDim).dot(dNU);
}
- gmm::dense_matrix<double> Ae3(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)
- Ae3(i, j) = mu * fluid->rho.evalD(*eU, phi, 0) * bei[i] * bejU[j];
+ Eigen::MatrixXd Ae3(e->nodes.size(), eU->nodes.size());
+ Ae3 = mu * fluid->rho.evalD(*eU, phi, 0) * bei * bejU;
// 4th term (mu*rhoU*grad_phi*grad_psi)
Eigen::MatrixXd Ae4;
@@ -309,10 +304,9 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
Ae4 *= mu * fluid->rho.eval(*eU, phi, 0);
// 5th term dmu*(rhoU-rho)*grad_phi*grad_psi
- double deltaRho = fluid->rho.eval(*eU, phi, 0) - fluid->rho.eval(*e, phi, 0);
Eigen::MatrixXd Ae5;
e->buildDK(Ae5, phi, fluid->mach);
- Ae5 *= deltaRho * dmu;
+ Ae5 *= (fluid->rho.eval(*eU, phi, 0) - fluid->rho.eval(*e, phi, 0)) * dmu;
// Assembly (supersonic)
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -322,12 +316,12 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- A(nodi->row, nodj->row) += Ae4(ii, jj) + Ae5(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae4(ii, jj) + Ae5(ii, jj)));
}
for (size_t jj = 0; jj < eU->nodes.size(); ++jj)
{
Node *nodj = eU->nodes[jj];
- A(nodi->row, nodj->row) += Ae3(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae3(ii, jj)));
}
}
}
@@ -339,10 +333,12 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- A(nodi->row, nodj->row) += Ae1(ii, jj) + Ae2(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae1(ii, jj) + Ae2(ii, jj)));
}
}
});
+ // Build Jacobian matrix without BCs
+ J.setFromTriplets(T.begin(), T.end());
// Kutta condition (3 steps): analytical derivatives
// Egality of normal flux through the wake
@@ -353,11 +349,17 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
tbb::parallel_do(wake->nodMap.begin(), wake->nodMap.end(), [&](std::pair<Node *, Node *> p) {
int idxUp = p.first->row;
int idxLw = p.second->row;
- for (auto node : pbl->msh->nodes)
+ // first grab upper row values...
+ std::vector<Eigen::Triplet<double>> tUp;
+ tUp.reserve(pbl->msh->nodes.size() / 10);
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(J, idxUp); it; ++it)
{
- A(idxLw, node->row) = A(idxLw, node->row) + A(idxUp, node->row);
- A(idxUp, node->row) = 0.;
+ tUp.push_back(Eigen::Triplet<double>(it.row(), it.col(), it.value()));
+ it.valueRef() = 0.;
}
+ // ...then add them to lower row values
+ for (auto t : tUp)
+ J.coeffRef(idxLw, t.col()) += t.value();
});
}
// Egality of velocity norm through the wake
@@ -377,12 +379,12 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
for (size_t jj = 0; jj < we->nColUp; ++jj)
{
Node *nodj = we->volUpE->nodes[jj];
- A(nodi->row, nodj->row) += Aue(ii, jj);
+ J.coeffRef(nodi->row, nodj->row) += Aue(ii, jj);
}
for (size_t jj = 0; jj < we->nColLw; ++jj)
{
Node *nodj = we->volLwE->nodes[jj];
- A(nodi->row, nodj->row) -= Ale(ii, jj);
+ J.coeffRef(nodi->row, nodj->row) -= Ale(ii, jj);
}
}
});
@@ -404,7 +406,7 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
for (size_t jj = 0; jj < ke->nCol; ++jj)
{
Node *nodj = ke->volE->nodes[jj];
- A(nodi->row, nodj->row) += Ae(ii, jj);
+ J.coeffRef(nodi->row, nodj->row) += Ae(ii, jj);
}
}
});
@@ -415,25 +417,36 @@ void Newton::buildJac(gmm::csr_matrix<double> &J)
{
for (auto nod : dBC->nodes)
{
- A.row(nod->row) = gmm::wsvector<double>(pbl->msh->nodes.size());
- A(nod->row, nod->row) = 1.0;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(J, nod->row); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
}
}
- gmm::copy(A, J);
+ // Clean matrix and turn to compressed row format
+ J.prune(0.);
+ J.makeCompressed();
+
if (verbose)
- std::cout << "J (" << gmm::mat_nrows(J) << "," << gmm::mat_ncols(J) << ") nnz=" << gmm::nnz(J) << "\n";
+ std::cout << "J (" << J.rows() << "," << J.cols() << ") nnz=" << J.nonZeros() << "\n";
}
/**
* @brief Build the residual vector
* @authors Adrien Crovato
*/
-void Newton::buildRes(std::vector<double> &R)
+void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
{
// Multithread
tbb::spin_mutex mutex;
+ // Reset residual
+ R.setZero();
+
// Full Potential Equation with upwind bias
auto fluid = pbl->medium;
tbb::parallel_do(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
@@ -462,7 +475,7 @@ void Newton::buildRes(std::vector<double> &R)
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- R[nodi->row] -= beU(ii);
+ R(nodi->row) -= beU(ii);
}
}
// Assembly (subsonic)
@@ -470,7 +483,7 @@ void Newton::buildRes(std::vector<double> &R)
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- R[nodi->row] -= be(ii);
+ R(nodi->row) -= be(ii);
}
});
// Apply freestream (Neumann) BCs
@@ -486,7 +499,7 @@ void Newton::buildRes(std::vector<double> &R)
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- R[nodi->row] += be(ii);
+ R(nodi->row) += be(ii);
}
});
}
@@ -503,7 +516,7 @@ void Newton::buildRes(std::vector<double> &R)
for (size_t ii = 0; ii < p.first->nodes.size(); ++ii)
{
Node *nodi = p.first->nodes[ii];
- R[nodi->row] += be(ii);
+ R(nodi->row) += be(ii);
}
});
}
@@ -518,8 +531,8 @@ void Newton::buildRes(std::vector<double> &R)
tbb::parallel_do(wake->nodMap.begin(), wake->nodMap.end(), [&](std::pair<Node *, Node *> p) {
int idxUp = p.first->row;
int idxLw = p.second->row;
- R[idxLw] = R[idxLw] + R[idxUp];
- R[idxUp] = 0.;
+ R(idxLw) += R[idxUp];
+ R(idxUp) = 0.;
});
}
// Egality of velocity norm through the wake
@@ -536,8 +549,8 @@ void Newton::buildRes(std::vector<double> &R)
for (size_t ii = 0; ii < we->nRow; ++ii)
{
Node *nodi = we->wkN[ii].second;
- R[nodi->row] -= bue(ii);
- R[nodi->row] += ble(ii);
+ R(nodi->row) -= bue(ii);
+ R(nodi->row) += ble(ii);
}
});
}
@@ -556,7 +569,7 @@ void Newton::buildRes(std::vector<double> &R)
for (size_t ii = 0; ii < ke->nRow; ++ii)
{
Node *nodi = ke->teN[ii].second;
- R[nodi->row] -= be(ii);
+ R(nodi->row) -= be(ii);
}
});
}
@@ -565,11 +578,11 @@ void Newton::buildRes(std::vector<double> &R)
for (auto dBC : pbl->dBCs)
{
for (auto nod : dBC->nodes)
- R[nod->row] = 0.;
+ R(nod->row) = 0.;
}
if (verbose)
- std::cout << "R (" << gmm::vect_size(R) << ")\n";
+ std::cout << "R (" << R.size() << ")\n";
}
/**
diff --git a/flow/src/wNewton.h b/flow/src/wNewton.h
index 86f287f11cd9d541b995561db206abcf22cfd9ba..bfa19aa6cdeac258230610d2d26cbe567f655c18 100644
--- a/flow/src/wNewton.h
+++ b/flow/src/wNewton.h
@@ -23,7 +23,7 @@
#include <iostream>
#include <vector>
#include <memory>
-#include "gmm_gmm.h"
+#include <Eigen/Sparse>
namespace flow
{
@@ -56,8 +56,8 @@ public:
#endif
private:
- void buildJac(gmm::csr_matrix<double> &J);
- void buildRes(std::vector<double> &R);
+ void buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J);
+ void buildRes(Eigen::Map<Eigen::VectorXd> &R);
void findUpwd();
};
diff --git a/flow/src/wPicard.cpp b/flow/src/wPicard.cpp
index f472a7df272f1181d9c312d8a4d552e8051eed93..bbee1333e0b0ed894cff1d9c70ecdf320e9a5840 100644
--- a/flow/src/wPicard.cpp
+++ b/flow/src/wPicard.cpp
@@ -27,12 +27,10 @@
#include "wTag.h"
#include "wNode.h"
#include "wElement.h"
-
-#include "gmm_gmm.h"
-#include "gmm_extra.h"
-#include "std_extra.h"
#include "wMumpsInterface.h"
+#include <Eigen/Sparse>
+
#include <tbb/task_scheduler_init.h>
#include <tbb/tbb.h>
#include <tbb/parallel_do.h>
@@ -92,7 +90,9 @@ bool Picard::run()
int it = 0;
double resInit = 1.;
double absRes = 1., relRes = 1.;
- std::vector<double> phiOld(phi.size()); // old solution (for relaxation)
+ // Solution, residual and old solution
+ Eigen::Map<Eigen::VectorXd> phi_(phi.data(), phi.size()), rPhi_(rPhi.data(), rPhi.size());
+ Eigen::VectorXd phiOld(phi.size());
// Setup linear (inner) solver MUMPS
std::map<std::string, int> mumpsOpt;
@@ -107,19 +107,18 @@ bool Picard::run()
do
{
// Reinitialize SoE
- gmm::csr_matrix<double> A;
- std::vector<double> b(pbl->msh->nodes.size());
+ Eigen::SparseMatrix<double, Eigen::RowMajor> A(pbl->msh->nodes.size(), pbl->msh->nodes.size());
+ std::vector<double> b(pbl->msh->nodes.size()); // dummy b vector
+ Eigen::Map<Eigen::VectorXd> b_ = Eigen::Map<Eigen::VectorXd>(b.data(), b.size());
// Build system of equations A \phi = b
build(A, b);
// Compute residual
- std::vector<double> A_phi(pbl->msh->nodes.size());
- gmm::mult(A, phi, A_phi);
- gmm::add(A_phi, gmm::scaled(b, -1.0), rPhi);
+ rPhi_ = A * phi_ - b_;
if (it == 0)
- resInit = gmm::vect_norm2(rPhi);
- absRes = gmm::vect_norm2(rPhi);
+ resInit = rPhi_.norm();
+ absRes = rPhi_.norm();
relRes = absRes / resInit;
// Compute aerodynamic coefficient on boundaries
@@ -141,17 +140,14 @@ bool Picard::run()
it++;
// Transfer old solution
- phiOld = phi;
+ phiOld = phi_;
// Solve linear SoE
- mumps.setSoE(A, b);
- mumps.analyze(); // could be done once at first iteration if structure does not change
- mumps.factorize();
- mumps.solve(phi);
+ mumps.setSoE(A, b_);
+ mumps.compute(phi_);
// Relaxation
- for (size_t i = 0; i < phi.size(); ++i)
- phi[i] = (1 - relax) * phiOld[i] + relax * phi[i];
+ phi_ = (1 - relax) * phiOld + relax * phi_;
} while (it < maxIt);
// Terminate MUMPS
@@ -185,13 +181,14 @@ bool Picard::run()
* @brief Build LHS matrix and RHS vector for Picard iteration
* @authors Adrien Crovato
*/
-void Picard::build(gmm::csr_matrix<double> &A, std::vector<double> &b)
+void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<double> &b)
{
// Multithread
tbb::spin_mutex mutex;
- // Matrix to be built
- gmm::row_matrix<gmm::wsvector<double>> localA(pbl->msh->nodes.size(), pbl->msh->nodes.size());
+ // List of triplets to build matrix
+ std::vector<Eigen::Triplet<double>> T;
+ T.reserve(pbl->msh->nodes.size());
auto fluid = pbl->medium;
tbb::parallel_do(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
@@ -210,18 +207,25 @@ void Picard::build(gmm::csr_matrix<double> &A, std::vector<double> &b)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- localA(nodi->row, nodj->row) += Ae(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae(ii, jj)));
}
}
});
+ // Build matrix without BCs
+ A.setFromTriplets(T.begin(), T.end());
// Apply Dirichlet BCs to A matrix
for (auto dBC : pbl->dBCs)
{
for (auto nod : dBC->nodes)
{
- localA.row(nod->row) = gmm::wsvector<double>(pbl->msh->nodes.size());
- localA(nod->row, nod->row) = 1.0;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, nod->row); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
}
}
@@ -273,12 +277,19 @@ void Picard::build(gmm::csr_matrix<double> &A, std::vector<double> &b)
tbb::parallel_do(wake->nodMap.begin(), wake->nodMap.end(), [&](std::pair<Node *, Node *> p) {
int idxUp = p.first->row;
int idxLw = p.second->row;
- for (auto node : pbl->msh->nodes)
+ // first grab upper row values...
+ std::vector<Eigen::Triplet<double>> tUp;
+ tUp.reserve(pbl->msh->nodes.size() / 10);
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(A, idxUp); it; ++it)
{
- localA(idxLw, node->row) = localA(idxLw, node->row) + localA(idxUp, node->row);
- localA(idxUp, node->row) = 0.;
+ tUp.push_back(Eigen::Triplet<double>(it.row(), it.col(), it.value()));
+ it.valueRef() = 0.;
}
- b[idxLw] = b[idxLw] + b[idxUp];
+ // ...then add them to lower row values
+ for (auto t : tUp)
+ A.coeffRef(idxLw, t.col()) += t.value();
+ // same for rhs
+ b[idxLw] += b[idxUp];
b[idxUp] = 0.;
});
}
@@ -299,12 +310,12 @@ void Picard::build(gmm::csr_matrix<double> &A, std::vector<double> &b)
for (size_t jj = 0; jj < we->volUpE->nodes.size(); ++jj)
{
Node *nodj = we->volUpE->nodes[jj];
- localA(nodi->row, nodj->row) += Aue(ii, jj);
+ A.coeffRef(nodi->row, nodj->row) += Aue(ii, jj);
}
for (size_t jj = 0; jj < we->volLwE->nodes.size(); ++jj)
{
Node *nodj = we->volLwE->nodes[jj];
- localA(nodi->row, nodj->row) -= Ale(ii, jj);
+ A.coeffRef(nodi->row, nodj->row) -= Ale(ii, jj);
}
}
});
@@ -326,17 +337,20 @@ void Picard::build(gmm::csr_matrix<double> &A, std::vector<double> &b)
for (size_t jj = 0; jj < ke->nCol; ++jj)
{
Node *nodj = ke->volE->nodes[jj];
- localA(nodi->row, nodj->row) += Ae(ii, jj);
+ A.coeffRef(nodi->row, nodj->row) += Ae(ii, jj);
}
}
});
}
- gmm::copy(localA, A);
+ // Clean matrix and turn to compressed row format
+ A.prune(0.);
+ A.makeCompressed();
+
if (verbose)
{
- std::cout << "A (" << gmm::mat_nrows(A) << "," << gmm::mat_ncols(A) << ") nnz=" << gmm::nnz(A) << "\n";
- std::cout << "b (" << gmm::vect_size(b) << ")\n";
+ std::cout << "A (" << A.rows() << "," << A.cols() << ") nnz=" << A.nonZeros() << "\n";
+ std::cout << "b (" << b.size() << ")\n";
}
}
diff --git a/flow/src/wPicard.h b/flow/src/wPicard.h
index 525d236596f0221149a84da0d34a21852ea79b0c..4237f8509aa70c2431ffb70baec1782029b5940d 100644
--- a/flow/src/wPicard.h
+++ b/flow/src/wPicard.h
@@ -23,7 +23,7 @@
#include <iostream>
#include <vector>
#include <memory>
-#include "gmm_gmm.h"
+#include <Eigen/Sparse>
namespace flow
{
@@ -48,7 +48,7 @@ public:
#endif
private:
- void build(gmm::csr_matrix<double> &A, std::vector<double> &b);
+ void build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<double> &b);
};
} // namespace flow
diff --git a/heat/src/wSolver.cpp b/heat/src/wSolver.cpp
index 88e802ada22f80e6c0828e8b411069942525e1b5..7fb301cf62d2fd40fafae92f21134c97caca13ed 100644
--- a/heat/src/wSolver.cpp
+++ b/heat/src/wSolver.cpp
@@ -31,10 +31,9 @@
#include "wPeriodic.h"
#include "wResults.h"
#include "wMshExport.h"
-#include "gmm_gmm.h"
-#include <gmm/gmm_MUMPS_interface.h>
-#include "gmm_extra.h"
-#include "std_extra.h"
+#include "wMumpsInterface.h"
+
+#include <Eigen/Sparse>
#include <typeinfo>
#include <algorithm>
#include <iomanip>
@@ -95,8 +94,8 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
n->row = i++;
// setup results
- T1.resize(msh->nodes.size());
- gmm::clear(T1); // in case of restart
+ T1.resize(msh->nodes.size(), 0.); // in case of restart
+ Eigen::Map<Eigen::VectorXd> T1_(T1.data(), T1.size());
Results results;
results.scalars_at_nodes["T"] = &T1;
@@ -138,6 +137,10 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
}
tms["prepro"].stop(); // ----- end of preprocessing phase
+ // Setup linear (inner) solver MUMPS
+ std::map<std::string, int> mumpsOpt;
+ MumpsInterface mumps(mumpsOpt);
+
// --------------------------------------------------------------------
// Iteration LOOP
// --------------------------------------------------------------------
@@ -145,12 +148,13 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
int it = 0;
while (it <= itmax)
{
- std::vector<double> rhs(msh->nodes.size()); // RHS
+ std::vector<double> rhs(msh->nodes.size()); // dummy RHS
+ Eigen::Map<Eigen::VectorXd> rhs_ = Eigen::Map<Eigen::VectorXd>(rhs.data(), rhs.size()); // actual RHS
tms["Kass"].start();
// K matrix assembly
- gmm::csr_matrix<double> K;
- buildK(K, rhs);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> K(msh->nodes.size(), msh->nodes.size());
+ buildK(K, rhs_);
//tomatlab("K.txt", K);
//throw std::runtime_error("stop!");
tms["Kass"].stop();
@@ -159,8 +163,8 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
if (!pbl->pdic)
{
- builds(rhs); // sources
- buildq(rhs); // prescribed fluxes (=sources over lines)
+ builds(rhs_); // sources
+ buildq(rhs_); // prescribed fluxes (=sources over lines)
// BCs
for (auto dBC : pbl->dBCs) // apply Dirichlet BCs to RHS
@@ -173,18 +177,14 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
tms["residual"].start();
// method #1
-
- std::vector<double> qint(msh->nodes.size());
- gmm::mult(K, T1, qint); // calcul de q_int = K T1:
+ Eigen::VectorXd qint = K * T1_; // calcul de q_int = K T1:
// method #2
if (0)
{
- std::vector<double> qint2(msh->nodes.size());
+ Eigen::VectorXd qint2;
buildqint(qint2); // calcul de q_int par integration du flux sur l'elem
- std::vector<double> diff(msh->nodes.size());
- gmm::add(qint, gmm::scaled(qint2, -1.0), diff);
- double normdiff = gmm::vect_norm2(diff);
+ double normdiff = (qint - qint2).norm();
//if(verbose>1)
std::cout << "norm diff=" << normdiff << '\n';
}
@@ -194,10 +194,9 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
std::cout << i << ":" << qint[i] << "," << qint2[i] << "\n";
*/
- std::vector<double> res(msh->nodes.size()); // residual
- gmm::add(qint, gmm::scaled(rhs, -1.0), res);
- double rnorm2 = gmm::vect_norm2(res);
- double rnorminf = gmm::vect_norminf(res);
+ Eigen::VectorXd res = qint - rhs_; // residual
+ double rnorm2 = res.norm();
+ double rnorminf = res.lpNorm<Eigen::Infinity>();
tms["residual"].stop();
if (verbose > 0)
@@ -216,9 +215,13 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
tms["solve"].start();
if (verbose > 1)
std::cout << "solving...\n";
- bool ok = gmm::MUMPS_solve(K, T1, rhs);
- if (!ok)
- throw std::runtime_error("problem while solving system with MUMPS");
+ mumps.setSoE(K, rhs_);
+ if (it == 1)
+ {
+ mumps.analyze();
+ }
+ mumps.factorize();
+ mumps.solve(T1_);
tms["solve"].stop();
}
if (it > itmax)
@@ -227,6 +230,9 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
std::cout << " in less than " << itmax << " iterations!\n\n"; // TODO: generer une erreur!
}
+ // Terminate MUMPS
+ mumps.finalize();
+
// --------------------------------------------------------------------
// Postprocessing
// --------------------------------------------------------------------
@@ -260,7 +266,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
results.tensors_at_elems["k [tensor]"] = &kmoy;
double Vtot = 0.0;
- std::vector<double> qM(2);
+ Eigen::Vector2d qM = Eigen::Vector2d::Zero();
for (auto mat : pbl->media)
{
tbb::parallel_do(mat->tag->elems.begin(), mat->tag->elems.end(),
@@ -270,11 +276,8 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
//std::cout << "processing element #" << e->no << "\n";
// (flux moyen . volume) sur l'element
- Eigen::VectorXd qqV;
- e->computeqV(T1, qqV, mat->k);
- std::vector<double> qV(qqV.rows());
- for (size_t ii = 0; ii < qqV.rows(); ++ii)
- qV[ii] = qqV(ii);
+ Eigen::VectorXd qV;
+ e->computeqV(T1, qV, mat->k);
double V = e->computeV();
@@ -289,33 +292,31 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
gradT[i] = Pt(grad[0], grad[1]);
// mean q over the element
- kgradT[i] = Pt(qV[0] / V, qV[1] / V);
+ kgradT[i] = Pt(qV(0) / V, qV(1) / V);
// mean k over the element
- Eigen::MatrixXd kkmean;
- e->computeV(T1, kkmean, mat->k);
- gmm::dense_matrix<double> kmean(kkmean.rows(), kkmean.cols());
- for (size_t ii = 0; ii < kkmean.rows(); ++ii)
- for (size_t jj = 0; jj < kkmean.cols(); ++jj)
- kmean(ii, jj) = kkmean(ii, jj);
- gmm::scale(kmean, 1. / V);
+ Eigen::MatrixXd kmean;
+ e->computeV(T1, kmean, mat->k);
+ kmean /= V;
kmoy11[i] = kmean(0, 0);
kmoy22[i] = kmean(1, 1);
kmoy12[i] = kmean(0, 1);
kmoy21[i] = kmean(1, 0);
// tensor
- gmm::resize(kmean, 3, 3);
- gmm::copy(kmean, kmoy[i]);
+ kmean.conservativeResize(3, 3);
+ for (size_t ii = 0; ii < kmean.rows(); ++ii)
+ for (size_t jj = 0; jj < kmean.cols(); ++jj)
+ kmoy[i](ii, jj) = kmean(ii, jj);
//std::cout << kmean << '\n';
}
tbb::spin_mutex::scoped_lock lock(mutex);
Vtot += V;
- gmm::add(qV, qM); // qM = qM + qV
+ qM += qV; // qM = qM + qV
});
}
- gmm::scale(qM, 1. / Vtot);
+ qM /= Vtot;
if (verbose > 1)
{
@@ -326,11 +327,11 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
// periodic BCs
if (pbl->pdic)
{
- pbl->pdic->qM = Pt(qM[0], qM[1]);
+ pbl->pdic->qM = Pt(qM(0), qM(1));
if (verbose > 1)
{
- std::cout << "=> k_x = " << ((pbl->pdic->dTxM == 0.0) ? 0 : qM[0] / pbl->pdic->dTxM) << '\n';
- std::cout << "=> k_y = " << ((pbl->pdic->dTyM == 0.0) ? 0 : qM[1] / pbl->pdic->dTyM) << '\n';
+ std::cout << "=> k_x = " << ((pbl->pdic->dTxM == 0.0) ? 0 : qM(0) / pbl->pdic->dTxM) << '\n';
+ std::cout << "=> k_y = " << ((pbl->pdic->dTyM == 0.0) ? 0 : qM(1) / pbl->pdic->dTyM) << '\n';
}
}
@@ -352,12 +353,15 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
}
}
-void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
+void Solver::buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map<Eigen::VectorXd> &rhs)
{
tbb::spin_mutex mutex;
std::shared_ptr<MshData> msh = pbl->msh;
- gmm::row_matrix<gmm::wsvector<double>> K(msh->nodes.size(), msh->nodes.size());
+
+ // List of triplets to build matrix
+ std::vector<Eigen::Triplet<double>> T;
+ T.reserve(msh->nodes.size());
// periodic BCs
// (on va sommer les lignes des noeuds "d" sur celles des noeuds "n")
@@ -392,12 +396,8 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::MatrixXd KKe;
- e->buildK(KKe, T1, mat->k, false);
- gmm::dense_matrix<double> Ke(KKe.rows(), KKe.cols());
- for (size_t ii = 0; ii < KKe.rows(); ++ii)
- for (size_t jj = 0; jj < KKe.cols(); ++jj)
- Ke(ii, jj) = KKe(ii, jj);
+ Eigen::MatrixXd Ke;
+ e->buildK(Ke, T1, mat->k, false);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -407,7 +407,7 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- K(rows[nodi->row], nodj->row) += Ke(ii, jj);
+ T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ke(ii, jj)));
}
}
});
@@ -426,16 +426,12 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
{
if (verbose > 1)
std::cout << "\tprocessing " << *mat << '\n';
- Eigen::MatrixXd KKe; // bidon
+ Eigen::MatrixXd Ke; // bidon
std::for_each(mat->tag->elems.begin(), mat->tag->elems.end(),
[&](Element *e) {
if (e->type() != ELTYPE::TRI3)
return;
- e->buildK(KKe, T1, mat->k, true);
- gmm::dense_matrix<double> Ke(KKe.rows(), KKe.cols());
- for (size_t ii = 0; ii < KKe.rows(); ++ii)
- for (size_t jj = 0; jj < KKe.cols(); ++jj)
- Ke(ii, jj) = KKe(ii, jj);
+ e->buildK(Ke, T1, mat->k, true);
});
}
@@ -460,12 +456,8 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
[&](Element *e) {
if (e->type() != ELTYPE::TRI3)
return;
- Eigen::MatrixXd KKe;
- e->buildK(KKe, T1, mat->k, false);
- gmm::dense_matrix<double> Ke(KKe.rows(), KKe.cols());
- for (size_t ii = 0; ii < KKe.rows(); ++ii)
- for (size_t jj = 0; jj < KKe.cols(); ++jj)
- Ke(ii, jj) = KKe(ii, jj);
+ Eigen::MatrixXd Ke;
+ e->buildK(Ke, T1, mat->k, false);
// assembly
//tbb::spin_mutex::scoped_lock lock(mutex);
@@ -475,7 +467,7 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- K(rows[nodi->row], nodj->row) += Ke(ii, jj);
+ T.push_back(Eigen::Triplet<double>(rows[nodi->row], nodj->row, Ke(ii, jj)));
}
}
});
@@ -483,6 +475,8 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
if (verbose > 1)
std::cout << "assembly done.\n";
}
+ // Build matrix without BCs
+ K2.setFromTriplets(T.begin(), T.end());
// periodic BCs
if (pbl->pdic)
@@ -500,20 +494,24 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
for (auto np : pbl->pdic->LR)
{
- K.row(np.d->row) = gmm::wsvector<double>(msh->nodes.size()); // clears related row
+ // clears related row
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K2, np.d->row); it; ++it)
+ it.valueRef() = 0.;
- K(np.d->row, np.d->row) = 1.0;
- K(np.d->row, np.n->row) = -1.0;
- rhs[np.d->row] = lx * dTxM;
+ K2.coeffRef(np.d->row, np.d->row) = 1.0;
+ K2.coeffRef(np.d->row, np.n->row) = -1.0;
+ rhs(np.d->row) = lx * dTxM;
}
for (auto np : pbl->pdic->BT)
{
- K.row(np.d->row) = gmm::wsvector<double>(msh->nodes.size()); // clears related row
+ // clears related row
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K2, np.d->row); it; ++it)
+ it.valueRef() = 0.;
- K(np.d->row, np.d->row) = 1.0;
- K(np.d->row, np.n->row) = -1.0;
- rhs[np.d->row] = ly * dTyM;
+ K2.coeffRef(np.d->row, np.d->row) = 1.0;
+ K2.coeffRef(np.d->row, np.n->row) = -1.0;
+ rhs(np.d->row) = ly * dTyM;
}
// corners \theta_2,3,4 fct \theta_1
@@ -521,22 +519,24 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
for (auto i = 0; i < 4; ++i)
{
Node *nod = pbl->pdic->corners[i];
- K.row(nod->row) = gmm::wsvector<double>(msh->nodes.size()); // clears related row
+ // clears related row
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K2, nod->row); it; ++it)
+ it.valueRef() = 0.;
}
Node *n1 = pbl->pdic->corners[0];
for (auto i = 1; i < 4; ++i)
{
Node *ni = pbl->pdic->corners[i];
- K(ni->row, ni->row) = 1.0;
- K(ni->row, n1->row) = -1.0;
+ K2.coeffRef(ni->row, ni->row) = 1.0;
+ K2.coeffRef(ni->row, n1->row) = -1.0;
}
Node *n2 = pbl->pdic->corners[1];
Node *n3 = pbl->pdic->corners[2];
Node *n4 = pbl->pdic->corners[3];
- rhs[n2->row] = lx * dTxM;
- rhs[n3->row] = lx * dTxM + ly * dTyM;
- rhs[n4->row] = ly * dTyM;
+ rhs(n2->row) = lx * dTxM;
+ rhs(n3->row) = lx * dTxM + ly * dTyM;
+ rhs(n4->row) = ly * dTyM;
// thermal energy consistency condition (n1->row)
@@ -551,29 +551,25 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::VectorXd CCe;
+ Eigen::VectorXd Ce;
// ** Ce matrix => K matrix
Fct0C f(1.0);
- e->builds(CCe, T1, f);
- std::vector<double> Ce(CCe.rows());
- for (size_t ii = 0; ii < CCe.rows(); ++ii)
- Ce[ii] = CCe(ii);
-
+ e->builds(Ce, T1, f);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- K(n1->row, nodj->row) += (mat->rhoc) * Ce[jj];
+ K2.coeffRef(n1->row, nodj->row) += (mat->rhoc) * Ce(jj);
- Vtot += Ce[jj];
- rhocM_Vtot += (mat->rhoc) * Ce[jj];
+ Vtot += Ce(jj);
+ rhocM_Vtot += (mat->rhoc) * Ce(jj);
}
});
}
- rhs[n1->row] = rhocM_Vtot * TM;
+ rhs(n1->row) = rhocM_Vtot * TM;
if (verbose > 1)
std::cout << "Vtot=" << Vtot << "\n";
if (verbose > 1)
@@ -591,28 +587,35 @@ void Solver::buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs)
for (auto nod : dBC->nodes)
{
//std::cout << "processing " << *nod << "\n";
- K.row(nod->row) = gmm::wsvector<double>(msh->nodes.size()); // clears related row
- K(nod->row, nod->row) = 1.0;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K2, nod->row); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
}
}
if (verbose > 1)
std::cout << "done." << std::endl;
}
- // copy in a read-only CSR
- gmm::copy(K, K2);
+ // Clean matrix and turn to compressed row format
+ K2.prune(0.);
+ K2.makeCompressed();
if (verbose > 1)
{
- std::cout << "K (" << gmm::mat_nrows(K2) << "," << gmm::mat_ncols(K2) << ") nnz=" << gmm::nnz(K2) << "\n";
+ std::cout << "K (" << K2.rows() << "," << K2.cols() << ") nnz=" << K2.nonZeros() << "\n";
}
- //gmm::cout << K2 << '\n';
}
-void Solver::buildqint(std::vector<double> &qint)
+void Solver::buildqint(Eigen::VectorXd &qint)
{
tbb::spin_mutex mutex;
+ qint = Eigen::VectorXd::Zero(pbl->msh->nodes.size());
+
for (auto mat : pbl->media)
{
if (verbose > 1)
@@ -623,18 +626,15 @@ void Solver::buildqint(std::vector<double> &qint)
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::VectorXd qqe;
- e->computeqint(T1, qqe, mat->k);
- std::vector<double> qe(qqe.rows());
- for (size_t ii = 0; ii < qqe.rows(); ++ii)
- qe[ii] = qqe(ii);
+ Eigen::VectorXd qe;
+ e->computeqint(T1, qe, mat->k);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- qint[nodi->row] += qe[ii];
+ qint(nodi->row) += qe(ii);
}
});
}
@@ -649,7 +649,7 @@ void Solver::buildqint(std::vector<double> &qint)
{
for (auto nod : dBC->nodes)
{
- qint[nod->row] = T1[nod->row]; // pour comparer avec K.T1
+ qint(nod->row) = T1[nod->row]; // pour comparer avec K.T1
}
}
if (verbose > 1)
@@ -661,7 +661,7 @@ void Solver::buildqint(std::vector<double> &qint)
* @brief volume sources
*/
-void Solver::builds(std::vector<double> &s)
+void Solver::builds(Eigen::Map<Eigen::VectorXd> &s)
{
tbb::spin_mutex mutex;
@@ -677,13 +677,10 @@ void Solver::builds(std::vector<double> &s)
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::VectorXd sse;
+ Eigen::VectorXd se;
// ** se vector => s vector
- e->builds(sse, T1, src->f);
- std::vector<double> se(sse.rows());
- for (size_t ii = 0; ii < sse.rows(); ++ii)
- se[ii] = sse(ii);
+ e->builds(se, T1, src->f);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -691,14 +688,14 @@ void Solver::builds(std::vector<double> &s)
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- s[nodi->row] += se[ii];
+ s(nodi->row) += se(ii);
}
});
}
if (verbose > 1)
{
- std::cout << "s (" << gmm::vect_size(s) << ")\n";
+ std::cout << "s (" << s.size() << ")\n";
}
}
@@ -706,7 +703,7 @@ void Solver::builds(std::vector<double> &s)
* @brief terme BC Neumann (pourrait etre mergé avec "builds"... c'est idem)
* @todo merge this routine with builds
*/
-void Solver::buildq(std::vector<double> &q)
+void Solver::buildq(Eigen::Map<Eigen::VectorXd> &q)
{
tbb::spin_mutex mutex;
@@ -722,13 +719,10 @@ void Solver::buildq(std::vector<double> &q)
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::VectorXd qqe;
+ Eigen::VectorXd qe;
// ** qe vector => q vector
- e->builds(qqe, T1, bnd->f);
- std::vector<double> qe(qqe.rows());
- for (size_t ii = 0; ii < qqe.rows(); ++ii)
- qe[ii] = qqe(ii);
+ e->builds(qe, T1, bnd->f);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -736,14 +730,14 @@ void Solver::buildq(std::vector<double> &q)
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- q[nodi->row] += qe[ii];
+ q(nodi->row) += qe(ii);
}
});
}
if (verbose > 1)
{
- std::cout << "q (" << gmm::vect_size(q) << ")\n";
+ std::cout << "q (" << q.size() << ")\n";
}
}
diff --git a/heat/src/wSolver.h b/heat/src/wSolver.h
index d8cf7b8eb902b93bbdac23ec59086a9f0243be3d..d91f0f36813b0eb09138ecd4e54474d46183efa7 100644
--- a/heat/src/wSolver.h
+++ b/heat/src/wSolver.h
@@ -25,7 +25,7 @@
#include <iostream>
#include <vector>
#include <memory>
-#include "gmm_gmm.h"
+#include <Eigen/Sparse>
using namespace tbox;
@@ -67,10 +67,10 @@ public:
#endif
private:
- void buildK(gmm::csr_matrix<double> &K2, std::vector<double> &rhs);
- void builds(std::vector<double> &s);
- void buildq(std::vector<double> &s);
- void buildqint(std::vector<double> &qint);
+ void buildK(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, Eigen::Map<Eigen::VectorXd> &rhs);
+ void builds(Eigen::Map<Eigen::VectorXd> &s);
+ void buildq(Eigen::Map<Eigen::VectorXd> &s);
+ void buildqint(Eigen::VectorXd &qint);
};
} // namespace heat
diff --git a/heat/tests/fe2/Lmpi.py b/heat/tests/fe2/Lmpi.py
index fdb8c21aaf542cf39d6327d2d84046f40cad9127..baf2439b92cfffceb02638f74771c78414fbc7b7 100644
--- a/heat/tests/fe2/Lmpi.py
+++ b/heat/tests/fe2/Lmpi.py
@@ -57,7 +57,7 @@ class Macro(object):
self.solver = h.Solver(self.pbl)
self.solver.nthreads = 1
- self.solver.restol = 1e-11
+ self.solver.restol = 1e-10
self.solver.itmax = 1000
if __name__ == "__main__":
diff --git a/heat/tests/periodic/honeycomb.py b/heat/tests/periodic/honeycomb.py
index c0836368d4d9439a0acc4c03b5e325949ae136f1..1464d321263bdc7ca3e100126fec3c6ef209d32b 100755
--- a/heat/tests/periodic/honeycomb.py
+++ b/heat/tests/periodic/honeycomb.py
@@ -36,6 +36,7 @@ class MyMaterial(tbox.Fct2U):
k2 = 100.
k.resize(2, 2)
k[0,0] = k[1,1] = k1+(k2-k1)/(t2-t1)*(u-t1)
+ k[0,1] = k[1,0] = 0.
#k[0,0] = k[1,1] = k1 # simple case
#print "returning k=%f" %k
return k
diff --git a/heat/tests/periodic/lincomb.py b/heat/tests/periodic/lincomb.py
index a833b73fbb345cc5276f20786cfca70276506738..246dfcf65dc2cfde0a23d5bce4f6971c56d87916 100755
--- a/heat/tests/periodic/lincomb.py
+++ b/heat/tests/periodic/lincomb.py
@@ -45,7 +45,8 @@ class MyMaterial(tbox.Fct2XYZ):
import math
v = old_div((self.k1+self.k2),2)+ abs(self.k2-self.k1)/2*math.sin(2*math.pi/L*4*(x+old_div(y,2)))
k.resize(2, 2)
- k[0,0] = k[1,1] = v;
+ k[0,0] = k[1,1] = v
+ k[0,1] = k[1,0] = 0.
#k[0,0] = self.k1 # constant
#k[1,1] = self.k2
diff --git a/mirrors/src/wSolver.cpp b/mirrors/src/wSolver.cpp
index 69f0274516939b826dc969cc7d3ee6aad54471bf..2ab298b05f55f718267d2da65fd65c1e78c057fe 100644
--- a/mirrors/src/wSolver.cpp
+++ b/mirrors/src/wSolver.cpp
@@ -30,9 +30,6 @@
#include "wNode.h"
#include "wPt.h"
#include "wANSYSSolution.h"
-#include "std_extra.h"
-#include "gmm_gmm.h"
-#include <gmm/gmm_MUMPS_interface.h>
#include <map>
#include "wTag.h"
#include <math.h>
@@ -43,6 +40,9 @@
#include "wGaussHex8.h"
#include "wResults.h"
#include "wMshExport.h"
+#include "wMumpsInterface.h"
+
+#include <Eigen/Sparse>
#include <math.h>
@@ -142,10 +142,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
std::cout << "\tAssembly...\n";
// memory allocation
- gmm::row_matrix<gmm::wsvector<double>> Tmp_1(pbl_size, pbl_size);
- gmm::csc_matrix<double> K(pbl_size, pbl_size);
+ std::vector<Eigen::Triplet<double>> Trip; // list of triplets to build stiffness matrix
+ Trip.reserve(pbl_size);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> K(pbl_size, pbl_size);
std::vector<double> rhs(pbl_size);
std::vector<double> u(pbl_size);
+ Eigen::Map<Eigen::VectorXd> u_(u.data(), u.size()), rhs_(rhs.data(), rhs.size());
// loop on the media
for (auto m : pbl.media)
@@ -185,16 +187,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
{
// Element *e = current_tag.elems[index];
- Eigen::MatrixXd KK_e;
+ Eigen::MatrixXd K_e;
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
return; //if the element is not a tetrahedron or a hexahedron, it is not treated
- e->buildK_mechanic(KK_e, H);
- gmm::dense_matrix<double> K_e(KK_e.rows(), KK_e.cols());
- for (size_t ii = 0; ii < KK_e.rows(); ++ii)
- for (size_t jj = 0; jj < KK_e.cols(); ++jj)
- K_e(ii, jj) = KK_e(ii, jj);
+ e->buildK_mechanic(K_e, H);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -208,7 +206,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
for (auto iii = 0; iii < 3; ++iii)
for (auto jjj = 0; jjj < 3; ++jjj)
if (K_e(3 * ii + iii, 3 * jj + jjj) != 0)
- Tmp_1(3 * (nodi->row) + iii, 3 * (nodj->row) + jjj) += K_e(3 * ii + iii, 3 * jj + jjj);
+ Trip.push_back(Eigen::Triplet<double>(3 * (nodi->row) + iii, 3 * (nodj->row) + jjj, K_e(3 * ii + iii, 3 * jj + jjj)));
}
}
});
@@ -226,16 +224,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
// L matrix
//=============================tbb==================================
tbb::parallel_do(current_tag.elems.begin(), current_tag.elems.end(), [&](Element *e) {
- Eigen::MatrixXd LL_e;
+ Eigen::MatrixXd L_e;
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
return; //if the element is not a tetrahedron or a hexahedron, it is not treated
- e->buildL_thermic(LL_e, K_conductivity);
- gmm::dense_matrix<double> L_e(LL_e.rows(), LL_e.cols());
- for (size_t ii = 0; ii < LL_e.rows(); ++ii)
- for (size_t jj = 0; jj < LL_e.cols(); ++jj)
- L_e(ii, jj) = LL_e(ii, jj);
+ e->buildL_thermic(L_e, K_conductivity);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -247,7 +241,7 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
tbb::spin_mutex::scoped_lock lock(mutex);
//=====================================//
if (L_e(ii, jj) != 0)
- Tmp_1(T_first_index + (nodi->row), T_first_index + (nodj->row)) += L_e(ii, jj);
+ Trip.push_back(Eigen::Triplet<double>(T_first_index + (nodi->row), T_first_index + (nodj->row), L_e(ii, jj)));
}
}
});
@@ -272,16 +266,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
// S matrix
//=============================tbb==================================
tbb::parallel_do(current_tag.elems.begin(), current_tag.elems.end(), [&](Element *e) {
- Eigen::MatrixXd SS_e;
+ Eigen::MatrixXd S_e;
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
return; //if the element is not a tetrahedron or a hexahedron, it is not treated
- e->buildS_thermomechanic(SS_e, D);
- gmm::dense_matrix<double> S_e(SS_e.rows(), SS_e.cols());
- for (size_t ii = 0; ii < SS_e.rows(); ++ii)
- for (size_t jj = 0; jj < SS_e.cols(); ++jj)
- S_e(ii, jj) = SS_e(ii, jj);
+ e->buildS_thermomechanic(S_e, D);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -294,13 +284,15 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
//=====================================//
for (size_t iii = 0; iii < 3; ++iii)
if (S_e(3 * ii + iii, jj) != 0)
- Tmp_1(3 * (nodi->row) + iii, T_first_index + (nodj->row)) += S_e(3 * ii + iii, jj);
+ Trip.push_back(Eigen::Triplet<double>(3 * (nodi->row) + iii, T_first_index + (nodj->row), S_e(3 * ii + iii, jj)));
}
}
});
//=============================tbb==================================
}
}
+ // Build matrix without BCs
+ K.setFromTriplets(Trip.begin(), Trip.end());
tms["Kass"].stop();
// ------------------------------------------
@@ -321,16 +313,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
//loop on the elements
for (auto e : current_tag.elems)
{
- Eigen::MatrixXd NN_e;
+ Eigen::MatrixXd N_e;
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a tetrahedron or a hexahedron, it is not treated
- e->buildM(NN_e);
- gmm::dense_matrix<double> N_e(NN_e.rows(), NN_e.cols());
- for (size_t ii = 0; ii < NN_e.rows(); ++ii)
- for (size_t jj = 0; jj < NN_e.cols(); ++jj)
- N_e(ii, jj) = NN_e(ii, jj);
+ e->buildM(N_e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -370,16 +358,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
//loop on the elements
for (auto e : current_tag.elems)
{
- Eigen::MatrixXd NN_e;
+ Eigen::MatrixXd N_e;
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a triangle or a quadrangle, it is not treated
- e->buildS(NN_e);
- gmm::dense_matrix<double> N_e(NN_e.rows(), NN_e.cols());
- for (size_t ii = 0; ii < NN_e.rows(); ++ii)
- for (size_t jj = 0; jj < NN_e.cols(); ++jj)
- N_e(ii, jj) = NN_e(ii, jj);
+ e->buildS(N_e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -403,16 +387,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
//loop on the elements
for (auto e : current_tag.elems)
{
- Eigen::MatrixXd NN_e;
+ Eigen::MatrixXd N_e;
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a triangle or a quadrangle, it is not treated
- e->buildS(NN_e);
- gmm::dense_matrix<double> N_e(NN_e.rows(), NN_e.cols());
- for (size_t ii = 0; ii < NN_e.rows(); ++ii)
- for (size_t jj = 0; jj < NN_e.cols(); ++jj)
- N_e(ii, jj) = NN_e(ii, jj);
+ e->buildS(N_e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -460,8 +440,13 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
if (b->which_dof.x[iii] == 1)
{
Node *nodi = e->nodes[ii];
- Tmp_1.row(3 * (nodi->row) + iii) = gmm::wsvector<double>(pbl_size); // clears related row
- Tmp_1(3 * (nodi->row) + iii, 3 * (nodi->row) + iii) = 1.0;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K, 3 * nodi->row + iii); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
rhs[3 * (nodi->row) + iii] = b->values.x[iii];
}
}
@@ -487,9 +472,14 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
Node *nodi = e->nodes[ii];
- Tmp_1.row(T_first_index + (nodi->row)) = gmm::wsvector<double>(pbl_size); // clears related row
- Tmp_1(T_first_index + (nodi->row), T_first_index + (nodi->row)) = 1.0;
- rhs[T_first_index + (nodi->row)] = b->values - T_ref;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K, T_first_index + nodi->row); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
+ rhs[T_first_index + nodi->row] = b->values - T_ref;
}
}
}
@@ -503,9 +493,9 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
// ------------------------------------------
tms["Kclean"].start();
std::cout << "\tCleaning matrices...\n";
- // clean K
- gmm::clean(Tmp_1, 1E-12);
- gmm::copy(Tmp_1, K);
+ // Clean matrix and turn to compressed row format
+ K.prune(0.);
+ K.makeCompressed();
tms["Kclean"].stop();
// ------------------------------------------
@@ -515,9 +505,10 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
// ------------------------------------------
tms["solve"].start();
std::cout << "\tSolving linear system...\n";
- bool ok = gmm::MUMPS_solve(K, u, rhs);
- if (!ok)
- throw std::runtime_error("problem while solving system with MUMPS");
+ std::map<std::string, int> mumpsOpt;
+ MumpsInterface mumps(mumpsOpt);
+ mumps.setSoE(K, rhs_);
+ mumps.compute(u_);
tms["solve"].stop();
// ------------------------------------------
@@ -665,16 +656,12 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
for (auto e : current_tag.elems)
{
- Eigen::MatrixXd NN_e;
+ Eigen::MatrixXd N_e;
if (e->type() != ELTYPE::HEX8 && e->type() != ELTYPE::TETRA4)
continue; //if the element is not a triangle or a quadrangle, it is not treated
- e->buildS(NN_e);
- gmm::dense_matrix<double> N_e(NN_e.rows(), NN_e.cols());
- for (size_t ii = 0; ii < NN_e.rows(); ++ii)
- for (size_t jj = 0; jj < NN_e.cols(); ++jj)
- N_e(ii, jj) = NN_e(ii, jj);
+ e->buildS(N_e);
for (size_t ii = 0; ii < e->nodes.size(); ++ii)
{
@@ -722,12 +709,10 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
Z_10[ii] = sqrt(8) * pow(rho[ii], 3.) * cos(3 * theta[ii]);
}
- gmm::row_matrix<gmm::wsvector<double>> W(msh->nodes.size(), msh->nodes.size());
- gmm::dense_matrix<double> Z(msh->nodes.size(), 10);
- gmm::row_matrix<gmm::wsvector<double>> TMP_Z(10, msh->nodes.size());
- gmm::dense_matrix<double> A(10, 10);
- std::vector<double> rhs_Z(10); // RHS
- std::vector<double> c_Z(10);
+ Eigen::MatrixXd W = Eigen::MatrixXd::Zero(msh->nodes.size(), msh->nodes.size());
+ Eigen::MatrixXd Z(msh->nodes.size(), 10);
+ Eigen::MatrixXd A(10, 10);
+ Eigen::VectorXd rhs_Z(10); // RHS
for (size_t ii = 0; ii < msh->nodes.size(); ++ii)
{
@@ -745,11 +730,9 @@ void Solver::start(std::shared_ptr<MshExport> mshWriter)
Z(ii, 9) = Z_10[ii];
}
- gmm::mult(gmm::transposed(Z), W, TMP_Z);
- gmm::mult(TMP_Z, gmm::col_vector(n_displacement), gmm::col_vector(rhs_Z));
- gmm::mult(TMP_Z, Z, A);
-
- gmm::lu_solve(A, c_Z, rhs_Z);
+ A = Z.transpose() * W * Z;
+ rhs_Z = Z.transpose() * W * Eigen::Map<Eigen::VectorXd>(n_displacement.data(), msh->nodes.size());
+ Eigen::VectorXd c_Z = A.partialPivLu().solve(rhs_Z);
std::cout << c_Z;
diff --git a/mrstlnos/src/Results_def.hpp b/mrstlnos/src/Results_def.hpp
index c4b34a8d2f870054c2190cd44cf8ee827639791e..3179124aae80da02c74cc18128838db7cb4c67ab 100644
--- a/mrstlnos/src/Results_def.hpp
+++ b/mrstlnos/src/Results_def.hpp
@@ -155,9 +155,6 @@ void Write_Results(
for (auto &k : strain_at_nodes)
gmm::resize(k, 3, 3);
- gmm::dense_matrix<double> zero_3x3(3, 3);
- gmm::clear(zero_3x3);
-
// Fill the result vectors
for (size_t ell = 0; ell < ensemble_size; ++ell)
@@ -208,8 +205,8 @@ void Write_Results(
displacement_y[i] = dy;
displacement_z[i] = dz;
displacement_u[i] = sqrt(dx * dx + dy * dy + dz * dz);
- gmm::copy(zero_3x3, strain_at_nodes[i]);
- gmm::copy(zero_3x3, stress_at_nodes[i]);
+ gmm::clear(strain_at_nodes[i]);
+ gmm::clear(stress_at_nodes[i]);
shear_modulus_per_node[i] = ET::coeff(random_field->operator()(pbl.msh->nodes[nodesWO_LO_to_GO[i]]->pos.x[0], pbl.msh->nodes[nodesWO_LO_to_GO[i]]->pos.x[1], pbl.msh->nodes[nodesWO_LO_to_GO[i]]->pos.x[2]), ell);
displacement[i] = Pt(dx, dy, dz);
if (numPrimalDPN == 4)
@@ -259,8 +256,8 @@ void Write_Results(
}
G = ET::coeff(random_field->operator()(x, y, z), ell);
}
- gmm::dense_matrix<double> H(9, 9);
- gmm::dense_matrix<double> D(3, 3);
+ Eigen::MatrixXd H = Eigen::MatrixXd::Zero(9, 9);
+ Eigen::MatrixXd D = Eigen::MatrixXd::Zero(3, 3);
for (auto i = 0; i < 3; ++i)
for (auto j = 0; j < 3; ++j)
@@ -283,17 +280,17 @@ void Write_Results(
D(i, i) = -3. * kappa * beta;
}
std::vector<double> x_tmp(8), y_tmp(8), z_tmp(8), T_tmp(8);
- gmm::dense_matrix<double> Strain(3, 3);
- gmm::dense_matrix<double> Sigma(3, 3);
+ Eigen::MatrixXd Strain(3, 3);
+ Eigen::MatrixXd Sigma(3, 3);
auto e = current_tag->elems[i_elem];
- std::vector<gmm::dense_matrix<double>> Strain_at_node(e->nodes.size());
- std::vector<gmm::dense_matrix<double>> Sigma_at_node(e->nodes.size());
+ std::vector<Eigen::MatrixXd> Strain_at_node(e->nodes.size());
+ std::vector<Eigen::MatrixXd> Sigma_at_node(e->nodes.size());
for (auto &k : Strain_at_node)
- gmm::resize(k, 3, 3);
+ k = Eigen::MatrixXd::Zero(3, 3);
for (auto &k : Sigma_at_node)
- gmm::resize(k, 3, 3);
+ k = Eigen::MatrixXd::Zero(3, 3);
if (e->type() == ELTYPE::HEX8 || e->type() == ELTYPE::TETRA4)
{
@@ -307,17 +304,12 @@ void Write_Results(
T_tmp[i] = T[node_i];
else
T_tmp[i] = 0.;
-
- gmm::dense_matrix<double> Strain_tmp(3, 3);
- gmm::clear(Strain_tmp);
- gmm::copy(Strain_tmp, Strain_at_node[i]);
- gmm::copy(Strain_tmp, Sigma_at_node[i]);
}
strain_stress_x_y_z(e, Strain, Sigma, H, D, x_tmp, y_tmp, z_tmp, T_tmp, numPrimalDPN == 4);
- gmm::copy(Strain, strain[elems_LO_to_GO[e->no - 1]]);
- gmm::copy(Sigma, stress[elems_LO_to_GO[e->no - 1]]);
+ strain[elems_LO_to_GO[e->no - 1]] = Strain;
+ stress[elems_LO_to_GO[e->no - 1]] = Sigma;
G_vector[elems_LO_to_GO[e->no - 1]] = G;
@@ -330,8 +322,13 @@ void Write_Results(
LocalOrdinal node_i = nodesWO_GO_to_LO[e->nodes[i]->row];
if (ell == 0)
++adjacent_elems[node_i]; // Increases the number of adjacent elements by one only for the first sample
- gmm::add(Strain_at_node[i], strain_at_nodes[node_i]);
- gmm::add(Sigma_at_node[i], stress_at_nodes[node_i]);
+
+ for (size_t iii = 0; iii < 3; ++iii)
+ for (size_t jjj = 0; jjj < 3; ++jjj)
+ {
+ strain_at_nodes[node_i](iii, jjj) += Strain_at_node[i](iii, jjj);
+ stress_at_nodes[node_i](iii, jjj) += Sigma_at_node[i](iii, jjj);
+ }
}
}
}
diff --git a/mrstlnos/src/Stress_computation.h b/mrstlnos/src/Stress_computation.h
index 463ab1fd4e2aa4f44e9d5304f0d8d5a95113d1d2..906c927f4842acd80afbd77eb0889e9c278759c0 100644
--- a/mrstlnos/src/Stress_computation.h
+++ b/mrstlnos/src/Stress_computation.h
@@ -12,19 +12,15 @@
/**
* @brief Those functions compute the stress and strain at nodes and in average per element.
- * They are currently implemented using gmm data structure but should be refactored using Kokkos:view
+ * They are currently implemented using Eigen data structure but should be refactored using Kokkos:view
* such that it will be possible to use parallel for in the future.
* @authors Kim Liegeois
*/
-void strain_at_gp(size_t n_n, Mem &mem, Cache &cache, size_t a, gmm::dense_matrix<double> u_e, gmm::dense_matrix<double> &Strain)
+void strain_at_gp(size_t n_n, Mem &mem, Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
{
- Eigen::MatrixXd const &JJ = mem.getJinv(a);
-
- gmm::dense_matrix<double> du_d_Phi(3, 3);
- gmm::dense_matrix<double> Tmp(3, 3);
- gmm::dense_matrix<double> Tmp2(3, 3);
-
+ Eigen::MatrixXd const &J = mem.getJinv(a);
+ Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
for (size_t i = 0; i < n_n; ++i)
{
std::vector<double> &dffi = cache.dff[a][i];
@@ -32,23 +28,12 @@ void strain_at_gp(size_t n_n, Mem &mem, Cache &cache, size_t a, gmm::dense_matri
for (auto j = 0; j < 3; ++j)
du_d_Phi(j, k) += dffi[j] * u_e(k, i);
}
-
- /* TEMPO EIGEN/GMM */
- gmm::dense_matrix<double> J(JJ.rows(), JJ.cols());
- for (size_t i = 0; i < J.nrows(); ++i)
- for (size_t j = 0; j < J.ncols(); ++j)
- J(i, j) = JJ(i, j);
-
- gmm::mult(J, du_d_Phi, Tmp);
- gmm::add(gmm::transposed(Tmp), Tmp, Tmp2);
- gmm::copy(gmm::scaled(Tmp2, 0.5), Tmp);
- gmm::add(Tmp, Strain);
+ Strain += (J * du_d_Phi).transpose() * (J * du_d_Phi) * 0.5;
}
-void compute_stress(size_t n_n, gmm::dense_matrix<double> Strain, gmm::dense_matrix<double> &Stress, const gmm::dense_matrix<double> &H, const gmm::dense_matrix<double> &D, const std::vector<double> &T, bool thermal)
+void compute_stress(size_t n_n, const Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &T, bool thermal)
{
- gmm::clear(Stress);
-
+ Stress = Eigen::Matrix3d::Zero();
for (size_t i = 0; i < 3; ++i)
for (size_t j = 0; j < 3; ++j)
for (size_t k = 0; k < 3; ++k)
@@ -66,7 +51,7 @@ void compute_stress(size_t n_n, gmm::dense_matrix<double> Strain, gmm::dense_mat
}
}
-void strain_stress_x_y_z(Element *elem, gmm::dense_matrix<double> &Strain, gmm::dense_matrix<double> &Stress, const gmm::dense_matrix<double> &H, const gmm::dense_matrix<double> &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
+void strain_stress_x_y_z(Element *elem, Eigen::MatrixXd &Strain, Eigen::MatrixXd &Stress, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
{
Cache &cache = elem->getVCache();
Mem &mem = elem->getVMem();
@@ -74,10 +59,7 @@ void strain_stress_x_y_z(Element *elem, gmm::dense_matrix<double> &Strain, gmm::
size_t n_gp = cache.dff.size();
size_t n_n = elem->nodes.size();
- gmm::dense_matrix<double> u_e(3, n_n);
- gmm::dense_matrix<double> Tmp(3, 3);
-
- gmm::clear(Strain);
+ Eigen::MatrixXd u_e(3, n_n);
for (size_t i = 0; i < n_n; ++i)
{
@@ -86,18 +68,17 @@ void strain_stress_x_y_z(Element *elem, gmm::dense_matrix<double> &Strain, gmm::
u_e(2, i) = z[i];
}
+ Strain = Eigen::Matrix3d::Zero();
for (size_t a = 0; a < n_gp; ++a)
{
strain_at_gp(n_n, mem, cache, a, u_e, Strain);
}
-
- gmm::copy(gmm::scaled(Strain, 1. / n_n), Tmp);
- gmm::copy(Tmp, Strain); // Prevent Level 2 warning: a conflict is possible in copy
+ Strain /= n_n;
compute_stress(n_n, Strain, Stress, H, D, T, thermal);
}
-void strain_stress_x_y_z_at_node(Element *elem, std::vector<gmm::dense_matrix<double>> &Strain_at_node, std::vector<gmm::dense_matrix<double>> &Stress_at_node, const gmm::dense_matrix<double> &H, const gmm::dense_matrix<double> &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
+void strain_stress_x_y_z_at_node(Element *elem, std::vector<Eigen::MatrixXd> &Strain_at_node, std::vector<Eigen::MatrixXd> &Stress_at_node, const Eigen::MatrixXd &H, const Eigen::MatrixXd &D, const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &T, bool thermal = false)
{
Cache &cache = elem->getVCache();
Mem &mem = elem->getVMem();
@@ -105,15 +86,15 @@ void strain_stress_x_y_z_at_node(Element *elem, std::vector<gmm::dense_matrix<do
size_t n_gp = cache.dff.size();
size_t n_n = elem->nodes.size();
- std::vector<gmm::dense_matrix<double>> Strain_at_gp(n_gp);
- std::vector<gmm::dense_matrix<double>> Stress_at_gp(n_gp);
- gmm::dense_matrix<double> u_e(3, n_n);
- gmm::dense_matrix<double> Strain(3, 3);
+ std::vector<Eigen::MatrixXd> Strain_at_gp(n_gp);
+ std::vector<Eigen::MatrixXd> Stress_at_gp(n_gp);
+ Eigen::MatrixXd u_e(3, n_n);
+ Eigen::MatrixXd Strain(3, 3);
for (auto &k : Strain_at_gp)
- gmm::resize(k, 3, 3);
+ k.resize(3, 3);
for (auto &k : Stress_at_gp)
- gmm::resize(k, 3, 3);
+ k.resize(3, 3);
for (size_t i = 0; i < n_n; ++i)
{
@@ -124,9 +105,9 @@ void strain_stress_x_y_z_at_node(Element *elem, std::vector<gmm::dense_matrix<do
for (size_t a = 0; a < n_gp; ++a)
{
- gmm::clear(Strain);
+ Strain = Eigen::Matrix3d::Zero();
strain_at_gp(n_n, mem, cache, a, u_e, Strain);
- gmm::copy(Strain, Strain_at_gp[a]);
+ Strain_at_gp[a] = Strain;
compute_stress(n_n, Strain, Stress_at_gp[a], H, D, T, thermal);
}
@@ -145,16 +126,10 @@ void strain_stress_x_y_z_at_node(Element *elem, std::vector<gmm::dense_matrix<do
elem->gp2node(values_at_gp, values_at_node);
for (size_t k = 0; k < n_n; ++k)
- {
- gmm::dense_matrix<double> Strain_tmp(3, 3);
- gmm::dense_matrix<double> Stress_tmp(3, 3);
for (size_t i = 0; i < 3; ++i)
for (size_t j = 0; j < 3; ++j)
{
- Strain_tmp(i, j) = values_at_node(k, i * 3 + j);
- Stress_tmp(i, j) = values_at_node(k, 9 + i * 3 + j);
+ Strain_at_node[k](i, j) = values_at_node(k, i * 3 + j);
+ Stress_at_node[k](i, j) = values_at_node(k, 9 + i * 3 + j);
}
- gmm::copy(Strain_tmp, Strain_at_node[k]);
- gmm::copy(Stress_tmp, Stress_at_node[k]);
- }
}
diff --git a/tbox/_src/tboxw.h b/tbox/_src/tboxw.h
index ce854492278e43ec67e41a4acaaefd42a8bb379f..0aed2c43f7871adc58e5d45631ecc42f51bc16ec 100644
--- a/tbox/_src/tboxw.h
+++ b/tbox/_src/tboxw.h
@@ -15,7 +15,7 @@
*/
#include "tbox.h"
-#include "gmm_extra.h"
+#include "eigen_extra.h"
#include "std_extra.h"
#include "wCache.h"
diff --git a/tbox/src/CMakeLists.txt b/tbox/src/CMakeLists.txt
index bc7e64d036c4c56e83f6bcd74c0914d37b82b955..682a8a677d081de3910001ce544a796f4e0d1cf1 100644
--- a/tbox/src/CMakeLists.txt
+++ b/tbox/src/CMakeLists.txt
@@ -68,12 +68,22 @@ ENDIF()
# -- Eigen --
FIND_PACKAGE(EIGEN REQUIRED)
TARGET_INCLUDE_DIRECTORIES(tbox PUBLIC ${EIGEN_INCLUDE_DIRS})
+TARGET_COMPILE_DEFINITIONS(tbox PUBLIC EIGEN_DONT_PARALLELIZE)
+IF(MKL_INCLUDE_DIRS AND MKL_LIBRARIES)
+ MESSAGE(STATUS "Linking Eigen with MKL")
+ TARGET_COMPILE_DEFINITIONS(tbox PUBLIC EIGEN_USE_MKL_ALL)
+ELSE()
+ MESSAGE(STATUS "Linking Eigen with ${BLA_VENDOR}")
+ TARGET_COMPILE_DEFINITIONS(tbox PUBLIC EIGEN_USE_BLAS)
+ #TARGET_COMPILE_DEFINITIONS(tbox PUBLIC EIGEN_USE_LAPACKE) # why this line fails?
+ENDIF()
+TARGET_LINK_LIBRARIES(tbox ${LAPACK_LIBRARIES})
# -- GMM --
FIND_PACKAGE(GMM REQUIRED)
TARGET_INCLUDE_DIRECTORIES(tbox PUBLIC ${GMM_INCLUDE_DIRS})
-TARGET_COMPILE_DEFINITIONS(tbox PUBLIC GMM_USES_LAPACK ${GMM_DEFINITIONS})
-TARGET_LINK_LIBRARIES(tbox ${LAPACK_LIBRARIES})
+#TARGET_COMPILE_DEFINITIONS(tbox PUBLIC GMM_USES_LAPACK ${GMM_DEFINITIONS})
+#TARGET_LINK_LIBRARIES(tbox ${LAPACK_LIBRARIES})
# -- MUMPS --
FIND_PACKAGE(MUMPS REQUIRED)
diff --git a/tbox/src/gmm_extra.cpp b/tbox/src/eigen_extra.cpp
similarity index 71%
rename from tbox/src/gmm_extra.cpp
rename to tbox/src/eigen_extra.cpp
index c16d1fc2cdf9b1d8103b4c822a3ee2e9f4a5fc6a..980a4fa419e1249587bba31db8efed99ea220af6 100644
--- a/tbox/src/gmm_extra.cpp
+++ b/tbox/src/eigen_extra.cpp
@@ -14,7 +14,7 @@
* limitations under the License.
*/
-#include "gmm_extra.h"
+#include "eigen_extra.h"
#include <fstream>
using namespace tbox;
@@ -26,17 +26,13 @@ using namespace tbox;
*/
TBOX_API void
-tomatlab(std::string const &fname, gmm::csr_matrix<double> const &A)
+tomatlab(std::string const &fname, Eigen::SparseMatrix<double> const &A)
{
std::ofstream f(fname.c_str());
- for (size_t i = 0; i < mat_nrows(A); ++i)
- for (size_t j = 0; j < mat_ncols(A); ++j)
- {
- double v = A(i, j);
- if (v != 0.0)
- f << i + 1 << ' ' << j + 1 << ' ' << v << '\n';
- }
+ for (size_t k = 0; k < A.outerSize(); ++k)
+ for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it)
+ f << it.row() + 1 << ' ' << it.col() + 1 << ' ' << it.value() << '\n';
f.close();
}
diff --git a/tbox/src/gmm_extra.h b/tbox/src/eigen_extra.h
similarity index 80%
rename from tbox/src/gmm_extra.h
rename to tbox/src/eigen_extra.h
index 9796bf7616daa52375098837b776098b99d1fc4c..09559a7352d5decefff4b6f9904e977d9d5cc0ad 100644
--- a/tbox/src/gmm_extra.h
+++ b/tbox/src/eigen_extra.h
@@ -14,18 +14,18 @@
* limitations under the License.
*/
-#ifndef GMM_EXTRA_H
-#define GMM_EXTRA_H
+#ifndef EIGEN_EXTRA_H
+#define EIGEN_EXTRA_H
#include "tbox.h"
-#include "gmm_gmm.h"
#include <string>
+#include <Eigen/Sparse>
namespace tbox
{
TBOX_API void
-tomatlab(std::string const &fname, gmm::csr_matrix<double> const &A);
+tomatlab(std::string const &fname, Eigen::SparseMatrix<double> const &A);
}
-#endif //GMM_EXTRA_H
+#endif //EIGEN_EXTRA_H
diff --git a/tbox/src/wMshDeform.cpp b/tbox/src/wMshDeform.cpp
index a54c9fa6ab87a4317ee9e2be0fdd76db05c68bf0..2b57c30a6efb31632acbcfa42441480c3a66e074 100644
--- a/tbox/src/wMshDeform.cpp
+++ b/tbox/src/wMshDeform.cpp
@@ -23,9 +23,10 @@
#include "wElement.h"
#include "wPt.h"
#include "wMem.h"
+#include "wMumpsInterface.h"
+
+#include <Eigen/Sparse>
-#include "gmm_gmm.h"
-#include "gmm_extra.h"
//#include <tbb/task_scheduler_init.h>
//#include <tbb/tbb.h>
//#include <tbb/parallel_do.h>
@@ -236,11 +237,12 @@ void MshDeform::deform()
// tbb::spin_mutex mutex;
// Stiffness matrix
- gmm::csr_matrix<double> K;
- // Temporary matrix
- gmm::row_matrix<gmm::wsvector<double>> Kt(nDim * mshSize, nDim * mshSize);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> K(nDim * mshSize, nDim * mshSize);
+ // List of triplets to build stifness matrix
+ std::vector<Eigen::Triplet<double>> T;
+ T.reserve(nDim * mshSize);
// force (RHS) and displacement (UKN) vectors
- std::vector<double> f(nDim * mshSize, 0.), q(nDim * mshSize, 0.);
+ Eigen::VectorXd f = Eigen::VectorXd::Zero(nDim * mshSize), q= Eigen::VectorXd::Zero(nDim * mshSize);
// Build stiffness matrix
std::cout << "--- Deforming mesh using linear elasticiy equations ---\n"
@@ -253,7 +255,7 @@ void MshDeform::deform()
Eigen::MatrixXd H;
if (nDim == 2)
{
- H = Eigen::Matrix3d::Zero();
+ H = Eigen::Matrix3d::Zero();
double nu = 0.; // Poisson ratio [-1, 0.5]
double E = 1 / e->computeV(); // Young modulus
H(0, 0) = 1 - nu;
@@ -265,7 +267,7 @@ void MshDeform::deform()
}
else
{
- H = Eigen::MatrixXd::Zero(6, 6);
+ H = Eigen::MatrixXd::Zero(6, 6);
double nu = 0.; // Poisson ratio [-1, 0.5]
double E = 1 / e->computeV(); // Young modulus
H(0, 0) = 1 - nu;
@@ -284,13 +286,8 @@ void MshDeform::deform()
}
// Elementary stiffness matrix
- Eigen::MatrixXd KKe;
- e->buildK_linElastic(KKe, H);
- gmm::dense_matrix<double> Ke(KKe.rows(), KKe.cols());
- for (size_t ii = 0; ii < KKe.rows(); ++ii)
- for (size_t jj = 0; jj < KKe.cols(); ++jj)
- Ke(ii, jj) = KKe(ii, jj);
-
+ Eigen::MatrixXd Ke;
+ e->buildK_linElastic(Ke, H);
// Assembly
// tbb::spin_mutex::scoped_lock lock(mutex);
@@ -304,13 +301,15 @@ void MshDeform::deform()
for (int n = 0; n < nDim; ++n)
{
size_t rowj = nDim * (e->nodes[j]->row) + n;
- Kt(rowi, rowj) += Ke(nDim * i + m, nDim * j + n);
+ T.push_back(Eigen::Triplet<double>(rowi, rowj, Ke(nDim * i + m, nDim * j + n)));
}
}
}
}
}
// });
+ // Build matrix without BCs
+ K.setFromTriplets(T.begin(), T.end());
// Apply "periodic condition" on internal boundaries
// -> for each node pair (a,b)
for (auto p : intBndPair)
@@ -319,17 +318,20 @@ void MshDeform::deform()
{
size_t row0 = nDim * (p.first->row) + m;
size_t row1 = nDim * (p.second->row) + m;
- // -> Add the lines a and b on b, and reset a line
- for (size_t j = 0; j < Kt.ncols(); j++)
+ // first grab "a" row values...
+ std::vector<Eigen::Triplet<double>> tUp;
+ tUp.reserve(nDim * mshSize / 10);
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K, row0); it; ++it)
{
- Kt(row1, j) = Kt(row1, j) + Kt(row0, j);
- Kt(row0, j) = 0.;
+ tUp.push_back(Eigen::Triplet<double>(it.row(), it.col(), it.value()));
+ it.valueRef() = 0.;
}
- f[row1] += f[row0];
- // -> Write the equation on a line (posA = posB)
- Kt(row0, row0) = 1.0;
- Kt(row0, row1) = -1.0;
- f[row0] = 0.;
+ // ...then add them to "b" row values
+ for (auto t : tUp)
+ K.coeffRef(row1, t.col()) += t.value();
+ // write equation on "a" line (pos"a" = pos"b")
+ K.coeffRef(row0, row0) = 1.0;
+ K.coeffRef(row0, row1) = -1.0;
}
}
// Apply Dirichlet condition on fixed boundaries
@@ -338,18 +340,26 @@ void MshDeform::deform()
for (int m = 0; m < nDim; ++m)
{
size_t rowi = nDim * (fxdBndNodes[i]->row) + m;
- Kt.row(rowi) = gmm::wsvector<double>(nDim * mshSize); // clear row
- Kt(rowi, rowi) = 1.;
- f[rowi] = 0.;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K, rowi); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
}
}
// Apply Dirichlet condition on symmetry boundaries
for (size_t i = 0; i < symBndNodes.size(); i++)
{
size_t rowi = nDim * (symBndNodes[i]->row) + blck;
- Kt.row(rowi) = gmm::wsvector<double>(nDim * mshSize); // clear row
- Kt(rowi, rowi) = 1.;
- f[rowi] = 0.;
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K, rowi); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
}
// Apply Dirichlet condition on moving boundaries
for (size_t i = 0; i < movBndNodes.size(); i++)
@@ -357,20 +367,25 @@ void MshDeform::deform()
for (int m = 0; m < nDim; ++m)
{
size_t rowi = nDim * (movBndNodes[i]->row) + m;
- Kt.row(rowi) = gmm::wsvector<double>(nDim * mshSize); // clear row
- Kt(rowi, rowi) = 1.;
- f[rowi] = movBndNodes[i]->pos.x[m] - initPos[nDim * i + m];
+ for (Eigen::SparseMatrix<double, Eigen::RowMajor>::InnerIterator it(K, rowi); it; ++it)
+ {
+ if (it.row() == it.col())
+ it.valueRef() = 1.;
+ else
+ it.valueRef() = 0.;
+ }
+ f(rowi) = movBndNodes[i]->pos.x[m] - initPos[nDim * i + m];
}
}
+ // Clean matrix and turn to compressed row format
+ K.prune(0.);
+ K.makeCompressed();
std::cout << "done" << std::endl;
-
// Solve the SoE
std ::cout << "Solving equations... " << std::flush;
- gmm::copy(Kt, K);
- gmm::iteration iter(1e-8);
- int restart = 50;
- gmm::ilu_precond<gmm::csr_matrix<double>> PR(K);
- gmm::gmres(K, q, f, PR, restart, iter);
+ Eigen::BiCGSTAB<Eigen::SparseMatrix<double, Eigen::RowMajor>, Eigen::IncompleteLUT<double>> solver; // something more efficient might exist
+ solver.compute(K);
+ q = solver.solve(f);
std::cout << "done" << std::endl;
// Update position
@@ -378,13 +393,13 @@ void MshDeform::deform()
for (auto n : msh->nodes)
{
for (int m = 0; m < nDim; m++)
- n->pos.x[m] += q[nDim * (n->row) + m];
+ n->pos.x[m] += q(nDim * (n->row) + m);
}
// Reset moving boundary position
for (auto n : movBndNodes)
{
for (int m = 0; m < nDim; m++)
- n->pos.x[m] -= q[nDim * (n->row) + m];
+ n->pos.x[m] -= q(nDim * (n->row) + m);
}
// Update element memory
for (auto e : fldElems)
diff --git a/tbox/src/wMumpsInterface.cpp b/tbox/src/wMumpsInterface.cpp
index e2488a7576b022aa75369d98e1ed1aac2673a667..74c1b45a25412c1e12fa4f3280c5c9425fd2b9f3 100644
--- a/tbox/src/wMumpsInterface.cpp
+++ b/tbox/src/wMumpsInterface.cpp
@@ -34,45 +34,38 @@ MumpsInterface::MumpsInterface(std::map<std::string, int> const &opts)
/**
* @brief Setup the system of equations
- *
- * Works only for row-major matrices
* @authors Adrien Crovato
*/
-void MumpsInterface::setSoE(gmm::csr_matrix<double> const &A, std::vector<double> const &b)
+void MumpsInterface::setSoE(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::VectorXd> const &b)
{
// Check that matrix A is square
- if (gmm::mat_nrows(A) != gmm::mat_ncols(A))
+ if (A.rows() != A.cols())
throw std::runtime_error("tbox::MumpsInterface::setMat: Matrix A is not square!\n");
// Build data structure
- // [AC] the matrix size should not change at each iteration, but can the number of nonzero element change?
- // if so, we need to clear and resize at each iteration, otherwise we can split this method
- size_t nz = gmm::nnz(A);
- size_t nr = gmm::mat_nrows(A);
- size_t nc = gmm::mat_ncols(A);
+ // the matrix size, structure and number of nonzero entries can change at each iteration!
+ size_t nz = A.nonZeros();
rows.clear();
cols.clear();
matx.clear();
- vect.clear();
+ vect.setZero();
rows.reserve(nz);
cols.reserve(nz);
matx.reserve(nz);
- vect.resize(b.size());
- gmm::copy(b, vect);
- // assume order of A to be row-major
- for (size_t i = 0; i < nr; ++i)
- for (size_t j = 0; j < nc; ++j)
- if (A(i, j) != 0.)
- {
- rows.push_back(static_cast<MUMPS_INT>(i + 1));
- cols.push_back(static_cast<MUMPS_INT>(j + 1));
- matx.push_back(A(i, j));
- }
+ vect = b;
+ // Get the nonzero indices and values
+ for (size_t k = 0; k < A.outerSize(); ++k)
+ for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it)
+ {
+ rows.push_back(static_cast<MUMPS_INT>(it.row() + 1));
+ cols.push_back(static_cast<MUMPS_INT>(it.col() + 1));
+ matx.push_back(it.value());
+ }
// Set sizes
- id.n = int(nr);
+ id.n = int(A.rows());
id.nz = int(nz);
// Set matrices
id.a = &matx[0];
- id.rhs = &vect[0];
+ id.rhs = vect.data();
// Set id
id.irn = &rows[0];
id.jcn = &cols[0];
@@ -100,12 +93,23 @@ void MumpsInterface::factorize()
* @brief Solve the SoE
* @authors Adrien Crovato, Romain Boman
*/
-void MumpsInterface::solve(std::vector<double> &x)
+void MumpsInterface::solve(Eigen::Map<Eigen::VectorXd> &x)
{
id.job = 3;
dmumps_c(&id);
check();
- gmm::copy(vect, x);
+ x = vect;
+}
+/**
+ * @brief Analyze the structure of the matrix, factorize it and solve the SoE
+ * @authors Adrien Crovato, Romain Boman
+ */
+void MumpsInterface::compute(Eigen::Map<Eigen::VectorXd> &x)
+{
+ id.job = 6;
+ dmumps_c(&id);
+ check();
+ x = vect;
}
/**
diff --git a/tbox/src/wMumpsInterface.h b/tbox/src/wMumpsInterface.h
index 494fb1ff9e06eba2fbc813138455687bf692be1d..9035bdd49b3d5786b3330f9fc2af493fb6455fd4 100644
--- a/tbox/src/wMumpsInterface.h
+++ b/tbox/src/wMumpsInterface.h
@@ -20,7 +20,7 @@
#include "tbox.h"
#include "wObject.h"
#include <dmumps_c.h>
-#include "gmm_gmm.h"
+#include <Eigen/Sparse>
#ifndef SWIG
@@ -34,10 +34,10 @@ namespace tbox
*/
/**
- * @brief MUMPS interface for tbox (gmm)
+ * @brief MUMPS interface for Eigen
*
* Reference: https://github.com/rboman/math0471/blob/master/mumps/main.cpp
- * Reference gmm_MUMPS_interface.h
+ * Reference gmm_MUMPS_interface.h (old)
* List of supported options (refer to setOptions() for description):
* - verbosity
* - ordering
@@ -52,7 +52,7 @@ private:
std::vector<MUMPS_INT> rows; ///< rows id
std::vector<MUMPS_INT> cols; ///< columns id
std::vector<double> matx; ///< matrix nonzero entries
- std::vector<double> vect; ///< rhs vector
+ Eigen::VectorXd vect; ///< rhs vector
void setOptions(std::map<std::string, int> const &opts);
void check();
@@ -60,10 +60,11 @@ private:
public:
MumpsInterface(std::map<std::string, int> const &opts);
- void setSoE(gmm::csr_matrix<double> const &A, std::vector<double> const &b);
+ void setSoE(Eigen::SparseMatrix<double> const &A, Eigen::Map<Eigen::VectorXd> const &b);
void analyze();
void factorize();
- void solve(std::vector<double> &x);
+ void solve(Eigen::Map<Eigen::VectorXd> &x);
+ void compute(Eigen::Map<Eigen::VectorXd> &x);
void finalize();
virtual void write(std::ostream &out) const override;
diff --git a/waves/src/wForwardEuler.cpp b/waves/src/wForwardEuler.cpp
index 77b9ac252648eb447fd86a1139db2b5a2126113c..db4527259962c58112960aa55486177d6ebf33a0 100644
--- a/waves/src/wForwardEuler.cpp
+++ b/waves/src/wForwardEuler.cpp
@@ -25,7 +25,6 @@
#include "wSource.h"
#include "wResults.h"
#include "wMshExport.h"
-#include "gmm_gmm.h"
using namespace tbox;
using namespace waves;
@@ -41,6 +40,7 @@ void ForwardEuler::start(std::shared_ptr<MshExport> mshWriter)
auto msh = pbl->msh;
u1.resize(msh->nodes.size());
v1.resize(msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> u1_(u1.data(), u1.size()), v1_(v1.data(), v1.size());
// assign a row for each node
int i = 0;
@@ -49,11 +49,12 @@ void ForwardEuler::start(std::shared_ptr<MshExport> mshWriter)
//std::for_each(msh->nodes.begin(), msh->nodes.end(), [&](Node *n) { n->row = i++; });
// matrix assembly
- gmm::csr_matrix<double> K;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> K(msh->nodes.size(), msh->nodes.size());
std::vector<double> Md(msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> Md_(Md.data(), Md.size());
buildKM(K, Md, u1);
- gmm::csr_matrix<double> S;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> S(msh->nodes.size(), msh->nodes.size());
buildS(S);
// setup results
@@ -66,7 +67,6 @@ void ForwardEuler::start(std::shared_ptr<MshExport> mshWriter)
std::cout << "\n\nTime integration...\n";
std::cout << *this << '\n';
- std::vector<double> tmp1(msh->nodes.size()), tmp2(msh->nodes.size());
double t = 0.0;
int nt = 1;
results.nt = nt;
@@ -86,16 +86,13 @@ void ForwardEuler::start(std::shared_ptr<MshExport> mshWriter)
nt++;
// equation #1
- gmm::add(u1, gmm::scaled(v1, dt), u1); // u1 = u0 + dt*v0
+ u1_ += v1_ * dt; // u1 = u0 + dt*v0
for (auto s : pbl->srcs)
s->apply(t, u1);
// equation #2
- gmm::mult(K, u1, tmp1);
- gmm::mult(S, v1, tmp2);
- for (size_t i = 0; i < gmm::vect_size(Md); ++i)
- v1[i] = v1[i] - dt / Md[i] * (tmp1[i] + tmp2[i]);
+ v1_.array() -= dt * (K * u1_ + S * v1_).array() / Md_.array();
if (nt % savefreq == 0)
{
diff --git a/waves/src/wNewmark.cpp b/waves/src/wNewmark.cpp
index f09f1d2d8c4ac1ee21af2c1e0c9fb0303b710112..0ae1088f6afd6f80540f9afdae636a95846db392 100644
--- a/waves/src/wNewmark.cpp
+++ b/waves/src/wNewmark.cpp
@@ -25,7 +25,6 @@
#include "wSource.h"
#include "wResults.h"
#include "wMshExport.h"
-#include "gmm_gmm.h"
using namespace tbox;
using namespace waves;
@@ -43,6 +42,7 @@ void Newmark::start(std::shared_ptr<MshExport> mshWriter)
auto msh = pbl->msh;
u1.resize(msh->nodes.size());
v1.resize(msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> u1_(u1.data(), u1.size()), v1_(v1.data(), v1.size());
// assign a row for each node
int i = 0;
@@ -50,11 +50,12 @@ void Newmark::start(std::shared_ptr<MshExport> mshWriter)
n->row = i++;
// matrix assembly
- gmm::csr_matrix<double> K;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> K(msh->nodes.size(), msh->nodes.size());
std::vector<double> Md(msh->nodes.size());
+ Eigen::Map<Eigen::VectorXd> Md_(Md.data(), Md.size());
buildKM(K, Md, u1);
- gmm::csr_matrix<double> S;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> S(msh->nodes.size(), msh->nodes.size());
buildS(S);
// setup results
@@ -67,7 +68,6 @@ void Newmark::start(std::shared_ptr<MshExport> mshWriter)
std::cout << "\n\nTime integration...\n";
std::cout << *this << '\n';
- std::vector<double> tmp1(msh->nodes.size()), tmp2(msh->nodes.size());
double t = 0.0;
int nt = 1;
results.nt = nt;
@@ -85,16 +85,13 @@ void Newmark::start(std::shared_ptr<MshExport> mshWriter)
nt++;
// equation #1
- gmm::add(u1, gmm::scaled(v1, dt), u1); // u1 = u0 + dt*v0
+ u1_ += v1_ * dt; // u1 = u0 + dt*v0
for (auto s : pbl->srcs)
s->apply(t, u1);
// equation #2
- gmm::mult(K, u1, tmp1);
- gmm::mult(S, v1, tmp2);
- for (size_t i = 0; i < gmm::vect_size(Md); ++i)
- v1[i] = v1[i] - dt / Md[i] * (tmp1[i] + tmp2[i]);
+ v1_.array() -= dt * (K * u1_ + S * v1_).array() / Md_.array();
if (nt % savefreq == 0)
{
diff --git a/waves/src/wRungeKutta.cpp b/waves/src/wRungeKutta.cpp
index 015c0f95eb72f39a1d8a2b3532abbb8867acee7c..30b2a4adfd2d44cf7d065e0fdd41429c62fe0aad 100644
--- a/waves/src/wRungeKutta.cpp
+++ b/waves/src/wRungeKutta.cpp
@@ -25,7 +25,6 @@
#include "wSource.h"
#include "wResults.h"
#include "wMshExport.h"
-#include "gmm_gmm.h"
using namespace tbox;
using namespace waves;
@@ -48,11 +47,11 @@ void RungeKutta::start(std::shared_ptr<MshExport> mshWriter)
n->row = i++;
// matrix assembly
- gmm::csr_matrix<double> K;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> K(msh->nodes.size(), msh->nodes.size());
std::vector<double> Md(msh->nodes.size());
buildKM(K, Md, u1);
- gmm::csr_matrix<double> S;
+ Eigen::SparseMatrix<double, Eigen::RowMajor> S(msh->nodes.size(), msh->nodes.size());
buildS(S);
// setup results
@@ -94,14 +93,14 @@ void RungeKutta::start(std::shared_ptr<MshExport> mshWriter)
t += dt;
nt++;
- predictor(up1, vp1, u1, v1, u1, v1, tmp1, tmp2, dt, Md, K, S, 1.0, t);
- predictor(up2, vp2, up1, vp1, u1, v1, tmp1, tmp2, dt, Md, K, S, 1.0, t);
- predictor(up3, vp3, up2, vp2, u1, v1, tmp1, tmp2, dt, Md, K, S, 1.0, t);
+ predictor(up1, vp1, u1, v1, u1, v1, dt, Md, K, S, 1.0, t);
+ predictor(up2, vp2, up1, vp1, u1, v1, dt, Md, K, S, 1.0, t);
+ predictor(up3, vp3, up2, vp2, u1, v1, dt, Md, K, S, 1.0, t);
- predictor(u2, v2, u1, v1, u1, v1, tmp1, tmp2, dt, Md, K, S, 1. / 2., t);
- predictor(u1, v1, up1, vp1, u2, v2, tmp1, tmp2, dt, Md, K, S, 1. / 3., t);
- predictor(u2, v2, up2, vp2, u1, v1, tmp1, tmp2, dt, Md, K, S, 1. / 8., t);
- predictor(u1, v1, up3, vp3, u2, v2, tmp1, tmp2, dt, Md, K, S, 1. / 24., t);
+ predictor(u2, v2, u1, v1, u1, v1, dt, Md, K, S, 1. / 2., t);
+ predictor(u1, v1, up1, vp1, u2, v2, dt, Md, K, S, 1. / 3., t);
+ predictor(u2, v2, up2, vp2, u1, v1, dt, Md, K, S, 1. / 8., t);
+ predictor(u1, v1, up3, vp3, u2, v2, dt, Md, K, S, 1. / 24., t);
//u1[snod->row] = sin(2*M_PI*20*t); //BC bis
//if(source) source->apply(t, u1);
@@ -128,13 +127,15 @@ void RungeKutta::start(std::shared_ptr<MshExport> mshWriter)
void RungeKutta::predictor(std::vector<double> &up2, std::vector<double> &vp2, // premier membre
std::vector<double> &up1, std::vector<double> &vp1, // predicteur
std::vector<double> &u0, std::vector<double> &v0, // sol au temps precedent
- std::vector<double> &tmp1, std::vector<double> &tmp2, double dt,
+ double dt,
std::vector<double> &Md,
- gmm::csr_matrix<double> &K,
- gmm::csr_matrix<double> &S, double ak, double t)
+ Eigen::SparseMatrix<double, Eigen::RowMajor> &K,
+ Eigen::SparseMatrix<double, Eigen::RowMajor> &S, double ak, double t)
{
+ // Define maps
+ Eigen::Map<Eigen::VectorXd> u0_(u0.data(), u0.size()), v0_(v0.data(), v0.size()), up1_(up1.data(), up1.size()), vp1_(vp1.data(), vp1.size()), up2_(up2.data(), up2.size()), vp2_(vp2.data(), vp2.size()), Md_(Md.data(), Md.size());
// equation #1
- gmm::add(u0, gmm::scaled(vp1, ak * dt), up2);
+ up2_ = u0_ + vp1_ * ak * dt;
//up2[snod->row] = sin(2*M_PI*20*t); //BC bis
//if(source) source->apply(t, up2);
@@ -142,8 +143,5 @@ void RungeKutta::predictor(std::vector<double> &up2, std::vector<double> &vp2, /
std::for_each(pbl->srcs.begin(), pbl->srcs.end(), [&](std::shared_ptr<Source> s) { s->apply(t, up2); });
// equation #2
- gmm::mult(K, up1, tmp1);
- gmm::mult(S, vp1, tmp2);
- for (size_t i = 0; i < gmm::vect_size(Md); ++i)
- vp2[i] = v0[i] - ak * dt / Md[i] * (tmp1[i] + tmp2[i]);
+ vp2_.array() = v0_.array() - ak * dt * (K * up1_ + S * vp1_).array() / Md_.array();
}
diff --git a/waves/src/wRungeKutta.h b/waves/src/wRungeKutta.h
index bdc622cdd67ff7eae3d46cf29ae224bfcd0222f9..f8ee1eeda3ec7242e227c5702d8d314628fd88d7 100644
--- a/waves/src/wRungeKutta.h
+++ b/waves/src/wRungeKutta.h
@@ -19,7 +19,7 @@
#include "waves.h"
#include "wTimeIntegration.h"
-#include "gmm_gmm.h"
+#include <Eigen/Sparse>
#include <vector>
#include <iostream>
@@ -43,10 +43,10 @@ private:
void predictor(std::vector<double> &up2, std::vector<double> &vp2, // premier membre
std::vector<double> &up1, std::vector<double> &vp1, // predicteur
std::vector<double> &u0, std::vector<double> &v0, // sol au temps precedent
- std::vector<double> &tmp1, std::vector<double> &tmp2, double dt,
+ double dt,
std::vector<double> &Md,
- gmm::csr_matrix<double> &K,
- gmm::csr_matrix<double> &S, double ak, double t);
+ Eigen::SparseMatrix<double, Eigen::RowMajor> &K,
+ Eigen::SparseMatrix<double, Eigen::RowMajor> &S, double ak, double t);
};
} // namespace waves
diff --git a/waves/src/wTimeIntegration.cpp b/waves/src/wTimeIntegration.cpp
index 91d62a73880f8e1844ddfae763cbc8a04e8890fd..8985c38c7f8280255f62926c5227161fe9508f13 100644
--- a/waves/src/wTimeIntegration.cpp
+++ b/waves/src/wTimeIntegration.cpp
@@ -26,7 +26,6 @@
#include "wBoundary.h"
#include "wFct0.h"
#include "wMshExport.h"
-#include "gmm_gmm.h"
#include <typeinfo>
#include <tbb/task_scheduler_init.h>
@@ -45,7 +44,7 @@ TimeIntegration::TimeIntegration(std::shared_ptr<Problem> _pbl) : pbl(_pbl)
ttot = 10.0;
dt = 1.0e-1 / 5;
savefreq = 5;
- nthreads = 4;
+ nthreads = 1;
wgmsh = Wgmsh::WG_NONE;
verbose = false;
@@ -88,13 +87,14 @@ void TimeIntegration::dummyIC()
}
}
-void TimeIntegration::buildS(gmm::csr_matrix<double> &S2)
+void TimeIntegration::buildS(Eigen::SparseMatrix<double, Eigen::RowMajor> &S2)
{
tbb::spin_mutex mutex;
tbb::task_scheduler_init init(nthreads);
auto msh = pbl->msh;
- gmm::row_matrix<gmm::wsvector<double>> S(msh->nodes.size(), msh->nodes.size());
+ std::vector<Eigen::Triplet<double>> T; // list of triplets to build S matrix
+ T.reserve(msh->nodes.size());
std::cout << "building S (TBB/lambda) using " << nthreads << " threads...\n";
@@ -108,13 +108,9 @@ void TimeIntegration::buildS(gmm::csr_matrix<double> &S2)
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::MatrixXd SSe;
+ Eigen::MatrixXd Se;
// ** Se matrix => S vector
- e->buildS(SSe);
- gmm::dense_matrix<double> Se(SSe.rows(), SSe.cols());
- for (size_t ii = 0; ii < SSe.rows(); ++ii)
- for (size_t jj = 0; jj < SSe.cols(); ++jj)
- Se(ii, jj) = SSe(ii, jj);
+ e->buildS(Se);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -125,21 +121,22 @@ void TimeIntegration::buildS(gmm::csr_matrix<double> &S2)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- S(nodi->row, nodj->row) += bnd->c * Se(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, bnd->c * Se(ii, jj)));
}
}
});
}
+ // Build, clean and turn to compressed row format
+ S2.setFromTriplets(T.begin(), T.end());
+ S2.prune(0.);
+ S2.makeCompressed();
- // copy in a read-only CSR
- gmm::copy(S, S2);
-
- std::cout << "S (" << gmm::mat_nrows(S2) << "," << gmm::mat_ncols(S2) << ") nnz=" << gmm::nnz(S2) << "\n";
+ std::cout << "S (" << S2.rows() << "," << S2.cols() << ") nnz=" << S2.nonZeros() << "\n";
chrono1.read();
std::cout << "[cpu] " << chrono1 << '\n';
}
-void TimeIntegration::buildKM(gmm::csr_matrix<double> &K2,
+void TimeIntegration::buildKM(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2,
std::vector<double> &Md, std::vector<double> const &u)
{
Timer chrono;
@@ -149,14 +146,15 @@ void TimeIntegration::buildKM(gmm::csr_matrix<double> &K2,
std::cout << "[cpu:buildKM] " << chrono << '\n';
}
-void TimeIntegration::buildKM_tbb_lambda(gmm::csr_matrix<double> &K2,
+void TimeIntegration::buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2,
std::vector<double> &Md, std::vector<double> const &u)
{
tbb::spin_mutex mutex;
tbb::task_scheduler_init init(nthreads);
auto msh = pbl->msh;
- gmm::row_matrix<gmm::wsvector<double>> K(msh->nodes.size(), msh->nodes.size());
+ std::vector<Eigen::Triplet<double>> T; // list of triplets to build K matrix
+ T.reserve(msh->nodes.size());
std::cout << "building K/M (TBB/lambda) using " << nthreads << " threads...\n";
@@ -169,22 +167,14 @@ void TimeIntegration::buildKM_tbb_lambda(gmm::csr_matrix<double> &K2,
return;
//std::cout << "processing element #" << e->no << "\n";
- Eigen::MatrixXd MMe;
- Eigen::MatrixXd KKe;
+ Eigen::MatrixXd Me;
+ Eigen::MatrixXd Ke;
// ** Me matrix => Md vector
- e->buildM(MMe);
- gmm::dense_matrix<double> Me(MMe.rows(), MMe.cols());
- for (size_t ii = 0; ii < MMe.rows(); ++ii)
- for (size_t jj = 0; jj < MMe.cols(); ++jj)
- Me(ii, jj) = MMe(ii, jj);
+ e->buildM(Me);
// ** Ke matrix => K matrix
- e->buildK(KKe, u, Fct0C(1.0));
- gmm::dense_matrix<double> Ke(KKe.rows(), KKe.cols());
- for (size_t ii = 0; ii < KKe.rows(); ++ii)
- for (size_t jj = 0; jj < KKe.cols(); ++jj)
- Ke(ii, jj) = KKe(ii, jj);
+ e->buildK(Ke, u, Fct0C(1.0));
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -195,28 +185,30 @@ void TimeIntegration::buildKM_tbb_lambda(gmm::csr_matrix<double> &K2,
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- K(nodi->row, nodj->row) += (mat->c * mat->c) * Ke(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, (mat->c * mat->c) * Ke(ii, jj)));
Md[nodi->row] += Me(ii, jj);
}
}
});
}
+ // Build, clean and turn to compressed row format
+ K2.setFromTriplets(T.begin(), T.end());
+ K2.prune(0.);
+ K2.makeCompressed();
- // copy in a read-only CSR
- gmm::copy(K, K2);
-
- std::cout << "K (" << gmm::mat_nrows(K2) << "," << gmm::mat_ncols(K2) << ") nnz=" << gmm::nnz(K2) << "\n";
+ std::cout << "S (" << K2.rows() << "," << K2.cols() << ") nnz=" << K2.nonZeros() << "\n";
chrono1.read();
std::cout << "[cpu] " << chrono1 << '\n';
}
-void TimeIntegration::build(MATTYPE type, gmm::csr_matrix<double> &A2)
+void TimeIntegration::build(MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2)
{
tbb::spin_mutex mutex;
tbb::task_scheduler_init init(nthreads); // TODO mettre ça ailleurs...
auto msh = pbl->msh;
- gmm::row_matrix<gmm::wsvector<double>> A(msh->nodes.size(), msh->nodes.size());
+ std::vector<Eigen::Triplet<double>> T; // list of triplets to build K matrix
+ T.reserve(msh->nodes.size());
std::cout << "building " << type << " using " << nthreads << " threads...\n";
@@ -227,12 +219,8 @@ void TimeIntegration::build(MATTYPE type, gmm::csr_matrix<double> &A2)
[&](Element *e) {
if (e->type() != ELTYPE::HEX8)
return;
- Eigen::MatrixXd AAe;
- e->build(type, AAe);
- gmm::dense_matrix<double> Ae(AAe.rows(), AAe.cols());
- for (size_t ii = 0; ii < AAe.rows(); ++ii)
- for (size_t jj = 0; jj < AAe.cols(); ++jj)
- Ae(ii, jj) = AAe(ii, jj);
+ Eigen::MatrixXd Ae;
+ e->build(type, Ae);
// assembly
tbb::spin_mutex::scoped_lock lock(mutex);
@@ -242,16 +230,17 @@ void TimeIntegration::build(MATTYPE type, gmm::csr_matrix<double> &A2)
for (size_t jj = 0; jj < e->nodes.size(); ++jj)
{
Node *nodj = e->nodes[jj];
- A(nodi->row, nodj->row) += Ae(ii, jj);
+ T.push_back(Eigen::Triplet<double>(nodi->row, nodj->row, Ae(ii, jj)));
}
}
});
}
+ // Build, clean and turn to compressed row format
+ A2.setFromTriplets(T.begin(), T.end());
+ A2.prune(0.);
+ A2.makeCompressed();
- // copy in a read-only CSR
- gmm::copy(A, A2);
-
- std::cout << type << " (" << gmm::mat_nrows(A2) << "," << gmm::mat_ncols(A2) << ") nnz=" << gmm::nnz(A2) << "\n";
+ std::cout << "S (" << A2.rows() << "," << A2.cols() << ") nnz=" << A2.nonZeros() << "\n";
chrono1.read();
std::cout << "[cpu] " << chrono1 << '\n';
}
diff --git a/waves/src/wTimeIntegration.h b/waves/src/wTimeIntegration.h
index 65ac1593ec0a7c864c797b46b9a3135281ca634b..fd9c4b4ea2b5c1689f832d51b7aa909b23b97d6e 100644
--- a/waves/src/wTimeIntegration.h
+++ b/waves/src/wTimeIntegration.h
@@ -24,7 +24,7 @@
#include <iostream>
#include <vector>
#include <memory>
-#include "gmm_gmm.h"
+#include <Eigen/Sparse>
namespace waves
{
@@ -71,14 +71,14 @@ public:
void dummyIC();
- void buildS(gmm::csr_matrix<double> &S2);
+ void buildS(Eigen::SparseMatrix<double, Eigen::RowMajor> &S2);
- void buildKM(gmm::csr_matrix<double> &K2, std::vector<double> &Md, std::vector<double> const &u);
- void buildKM_tbb_lambda(gmm::csr_matrix<double> &K2, std::vector<double> &Md, std::vector<double> const &u);
+ void buildKM(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, std::vector<double> &Md, std::vector<double> const &u);
+ void buildKM_tbb_lambda(Eigen::SparseMatrix<double, Eigen::RowMajor> &K2, std::vector<double> &Md, std::vector<double> const &u);
void setGUI(DisplayHook &hook) { dhook = &hook; }
- void build(tbox::MATTYPE type, gmm::csr_matrix<double> &A2);
+ void build(tbox::MATTYPE type, Eigen::SparseMatrix<double, Eigen::RowMajor> &A2);
void stop() { stopit = true; }