Merge pull request #553 from vressegu/devMattPullRequest

Dev matt pull request
ITHACA-FV · May 18, 2024 · 535c149 · 535c149
2 parents 4c8ec5f + d3a76d5
commit 535c149
Show file tree

Hide file tree

Showing 34 changed files with 329,238 additions and 27 deletions.
diff --git a/src/ITHACA_CORE/Foam2Eigen/Foam2Eigen.C b/src/ITHACA_CORE/Foam2Eigen/Foam2Eigen.C
@@ -422,6 +422,61 @@ template volTensorField Foam2Eigen::Eigen2field(
     volTensorField& field_in, Eigen::VectorXd& eigen_vector, bool correctBC);
 
 template <template <class> class PatchField, class GeoMesh>
+GeometricField<vector, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
+    GeometricField<vector, PatchField, GeoMesh>& field_in,
+    Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary)
+{
+    GeometricField<vector, PatchField, GeoMesh> field_out(field_in);
+
+    for (auto i = 0; i < field_out.size(); i++)
+    {
+        field_out.ref()[i][0] = eigen_vector(i);
+        field_out.ref()[i][1] = eigen_vector(i + field_out.size());
+        field_out.ref()[i][2] = eigen_vector(i + field_out.size() * 2);
+    }
+
+    for(unsigned int id = 0; id < field_out.boundaryField().size(); id++)
+    {
+        unsigned int idBSize = field_out.boundaryField()[id].size();
+        for (unsigned int ith_bcell = 0; ith_bcell < idBSize; ith_bcell++)
+        {
+            ITHACAutilities::assignBC(field_out, id, eigen_vector_boundary[id]);
+        }
+    }
+
+    return field_out;
+}
+
+template volVectorField Foam2Eigen::Eigen2field(
+    volVectorField& field_in, Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary);
+
+template<template<class> class PatchField, class GeoMesh>
+GeometricField<scalar, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
+    GeometricField<scalar, PatchField, GeoMesh>& field_in,
+    Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary)
+{
+    GeometricField<scalar, PatchField, GeoMesh> field_out(field_in);
+
+    for (auto i = 0; i < field_out.size(); i++)
+    {
+        field_out.ref()[i] = eigen_vector(i);
+    }
+
+    for(unsigned int id = 0; id < field_out.boundaryField().size(); id++)
+    {
+        for (unsigned int ith_bcell = 0; ith_bcell < field_out.boundaryField()[id].size(); ith_bcell++)
+        {
+            ITHACAutilities::assignBC(field_out, id, eigen_vector_boundary[id]);
+        }
+    }
+
+    return field_out;
+}
+
+template volScalarField Foam2Eigen::Eigen2field(
+    volScalarField& field_in, Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary);
+
+template<template<class> class PatchField, class GeoMesh>
 GeometricField<vector, PatchField, GeoMesh> Foam2Eigen::Eigen2field(
     GeometricField<vector, PatchField, GeoMesh>& field_in,
     Eigen::VectorXd& eigen_vector, bool correctBC)

diff --git a/src/ITHACA_CORE/Foam2Eigen/Foam2Eigen.H b/src/ITHACA_CORE/Foam2Eigen/Foam2Eigen.H
@@ -40,6 +40,7 @@
 #include "fvCFD.H"
 #include "IOmanip.H"
 #include "ITHACAassert.H"
+#include "ITHACAutilities.H"
 #include <tuple>
 #include <sys/stat.h>
 #pragma GCC diagnostic push
@@ -356,6 +357,43 @@ class Foam2Eigen
         static GeometricField<vector, PatchField, GeoMesh> Eigen2field(
             GeometricField<vector, PatchField, GeoMesh>& field,
             Eigen::VectorXd& eigen_vector, bool correctBC = true);
+
+        //----------------------------------------------------------------------
+        /// @brief         Convert a vector in Eigen format into an OpenFOAM
+        ///                scalar GeometricField
+        ///
+        /// @param[in/out] field         OpenFOAM GeometricField
+        /// @param[in]     eigen_vector  Vector in Eigen format
+        /// @param[in]     eigen_vector_boundary  Vector in Eigen format
+        ///
+        /// @tparam        PatchField    fvPatchField or fvsPatchField.
+        /// @tparam        GeoMesh       volMesh or surfaceMesh.
+        ///
+        /// @return        OpenFOAM GeometricField
+        ///
+        template<template<class> class PatchField, class GeoMesh>
+        static GeometricField<scalar, PatchField, GeoMesh> Eigen2field(
+            GeometricField<scalar, PatchField, GeoMesh>& field,
+            Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary);
+
+        //----------------------------------------------------------------------
+        /// @brief         Convert a vector in Eigen format into an OpenFOAM
+        ///                vector GeometricField
+        ///
+        /// @param[in/out] field         OpenFOAM GeometricField
+        /// @param[in]     eigen_vector  Vector in Eigen format
+        /// @param[in]     eigen_vector_boundary  Vector in Eigen format
+        ///
+        /// @tparam        PatchField    fvPatchField or fvsPatchField.
+        /// @tparam        GeoMesh       volMesh or surfaceMesh.
+        ///
+        /// @return        OpenFOAM GeometricField
+        ///
+        template<template<class> class PatchField, class GeoMesh>
+        static GeometricField<vector, PatchField, GeoMesh> Eigen2field(
+            GeometricField<vector, PatchField, GeoMesh>& field,
+            Eigen::VectorXd& eigen_vector, List<Eigen::VectorXd>& eigen_vector_boundary);
+
 
         //----------------------------------------------------------------------
         /// @brief         Convert a vector in Eigen format into an OpenFOAM

diff --git a/src/ITHACA_CORE/ITHACAutilities/ITHACAutilities.C b/src/ITHACA_CORE/ITHACAutilities/ITHACAutilities.C
@@ -133,4 +133,81 @@ List<T> combineList(List<List<T>>& doubleList)
 template List<label> combineList(List<List<label>>& doubleList);
 
 
+// Using the Eigen library, using the SVD decomposition method to solve the
+// matrix pseudo-inverse, the default error er is 0
+Eigen::MatrixXd pinv_eigen_based(Eigen::MatrixXd &origin, const float er) {
+  // perform svd decomposition
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd_holder(origin, Eigen::ComputeThinU | Eigen::ComputeThinV);
+  // Build SVD decomposition results
+  Eigen::MatrixXd U = svd_holder.matrixU();
+  Eigen::MatrixXd V = svd_holder.matrixV();
+  Eigen::MatrixXd D = svd_holder.singularValues();
+
+  // Build the S matrix
+  Eigen::MatrixXd S(V.cols(), U.cols());
+  S.setZero();
+
+  for (unsigned int i = 0; i < D.size(); ++i) {
+
+    if (D(i, 0) > er) {
+      S(i, i) = 1 / D(i, 0);
+    } else {
+      S(i, i) = 0;
+    }
+  }
+  return V * S * U.transpose();
+}
+
+
+Eigen::MatrixXd invertMatrix(Eigen::MatrixXd& matrixToInvert, const word inversionMethod){
+
+  Info << "Inversion method : " << inversionMethod << endl;
+  if(inversionMethod == "pinv_eigen_based"){
+    return pinv_eigen_based(matrixToInvert);
+  }
+  else if(inversionMethod == "direct"){
+    return matrixToInvert.inverse();
+  }
+  else if(inversionMethod == "fullPivLu"){
+    return matrixToInvert.fullPivLu().inverse();
+  }
+  else if(inversionMethod == "partialPivLu"){
+    return matrixToInvert.partialPivLu().inverse();
+  }
+  else if(inversionMethod == "householderQr"){
+    return matrixToInvert.householderQr().solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
+  }
+  else if(inversionMethod == "colPivHouseholderQr"){
+    return matrixToInvert.colPivHouseholderQr().inverse();
+  }
+  else if(inversionMethod == "fullPivHouseholderQr"){
+    return matrixToInvert.fullPivHouseholderQr().inverse();
+  }
+  else if(inversionMethod == "completeOrthogonalDecomposition"){
+    return matrixToInvert.completeOrthogonalDecomposition().pseudoInverse();
+  }
+  else if(inversionMethod == "jacobiSvd")
+  {
+    return matrixToInvert.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
+  }
+  else if(inversionMethod == "llt")
+  {
+    return matrixToInvert.llt().solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
+  }
+  else if(inversionMethod == "ldlt")
+  {
+    Eigen::LLT<Eigen::MatrixXd> lltOfA(matrixToInvert);
+    return lltOfA.solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
+  }
+  else if(inversionMethod == "bdcSvd")
+  {
+    return matrixToInvert.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(Eigen::MatrixXd::Identity(matrixToInvert.rows(), matrixToInvert.cols()));
+  } 
+  else{
+    Info << "Unkwown inversion method, solving with : completeOrthogonalDecomposition" << endl;
+    return matrixToInvert.completeOrthogonalDecomposition().pseudoInverse();
+  }
 }
+
+
+}
diff --git a/src/ITHACA_CORE/ITHACAutilities/ITHACAutilities.H b/src/ITHACA_CORE/ITHACAutilities/ITHACAutilities.H
@@ -125,6 +125,27 @@ bool isTurbulent();
 template<typename T>
 List<T> combineList(List<List<T>>& doubleList);
 
+//------------------------------------------------------------------------------
+/// @brief      Using the Eigen library, using the SVD decomposition method to solve the
+///             matrix pseudo-inverse, the default error er is 0
+///
+/// @param      origin  Matrix to invert
+/// @param      er      Error
+///
+/// @return     Inversed matrix
+///
+Eigen::MatrixXd pinv_eigen_based(Eigen::MatrixXd &origin, const float er = 0);
+
+//------------------------------------------------------------------------------
+/// @brief      Invert a matrix given the method name in the ITHACAdict
+///
+/// @param[in]      matrixToInvert  Matrix to invert
+/// @param[in]      methodName      Method 
+///
+/// @return         Inversed matrix
+///
+Eigen::MatrixXd invertMatrix(Eigen::MatrixXd& matrixToInvert, const word inversionMethod);
+
 
 };
 #endif
diff --git a/src/ITHACA_FOMPROBLEMS/Burgers/Burgers.C b/src/ITHACA_FOMPROBLEMS/Burgers/Burgers.C
@@ -98,6 +98,7 @@ void Burgers::truthSolve(word folder)
 
     while (_simple().loop())
     {
+        auto start = std::chrono::system_clock::now();
         Info << "Time = " << _runTime().timeName() << nl << endl;
 
         while (simple.correctNonOrthogonal())
@@ -112,6 +113,10 @@ void Burgers::truthSolve(word folder)
 
         phi = linearInterpolate(U) & mesh.Sf();
 
+        auto end = std::chrono::system_clock::now();
+        auto elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
+        Info << "Elapsed: " << elapsed.count() << " microseconds\n";
+
         if (checkWrite(runTime))
         {
             ITHACAstream::exportSolution(U, name(counter), folder + name(folderN));