[Solver] Semi-implicit Navier-Stokes solver

CMakeLists.txt

@@ -97,7 +97,6 @@ if (MKL_FOUND)
   endif()
 endif()
 
-
 # KaHIP
 if (MPI_FOUND)
   if (NO_KAHIP)
@@ -138,10 +137,12 @@ include_directories(BEFORE
   ${mpm_SOURCE_DIR}/include/functions/
   ${mpm_SOURCE_DIR}/include/generators/
   ${mpm_SOURCE_DIR}/include/io/
+  ${mpm_SOURCE_DIR}/include/linear_solvers/
   ${mpm_SOURCE_DIR}/include/loads_bcs/
   ${mpm_SOURCE_DIR}/include/materials/
   ${mpm_SOURCE_DIR}/include/particles/
   ${mpm_SOURCE_DIR}/include/solvers/
+  ${mpm_SOURCE_DIR}/include/utilities/
   ${mpm_SOURCE_DIR}/external/
   ${mpm_SOURCE_DIR}/tests/include/
 )
@@ -161,6 +162,7 @@ SET(mpm_src
   ${mpm_SOURCE_DIR}/src/io/logger.cc
   ${mpm_SOURCE_DIR}/src/io/partio_writer.cc
   ${mpm_SOURCE_DIR}/src/io/vtk_writer.cc
+  ${mpm_SOURCE_DIR}/src/linear_solver.cc
   ${mpm_SOURCE_DIR}/src/material.cc
   ${mpm_SOURCE_DIR}/src/mpm.cc
   ${mpm_SOURCE_DIR}/src/nodal_properties.cc

include/cell.h

@@ -214,6 +214,66 @@ class Cell {
   //! Return rank
   unsigned rank() const;
 
+  //! Assign free surface
+  //! \param[in] free_surface boolean indicating free surface cell
+  void assign_free_surface(bool free_surface) { free_surface_ = free_surface; };
+
+  //! Return free surface bool
+  //! \retval free_surface_ indicating free surface cell
+  bool free_surface() { return free_surface_; };
+
+  //! Assign volume traction to node
+  //! \param[in] volume_fraction cell volume fraction
+  void assign_volume_fraction(double volume_fraction) {
+    volume_fraction_ = volume_fraction;
+  };
+
+  //! Return cell volume fraction
+  //! \retval volume_fraction_ cell volume fraction
+  double volume_fraction() { return volume_fraction_; };
+
+  //! Map cell volume to the nodes
+  //! \param[in] phase to map volume
+  bool map_cell_volume_to_nodes(unsigned phase);
+
+  //! Initialize local elemental matrices
+  bool initialise_element_matrix();
+
+  //! Return local node indices
+  Eigen::VectorXi local_node_indices();
+
+  //! Return local laplacian
+  const Eigen::MatrixXd& laplacian_matrix() { return laplacian_matrix_; };
+
+  //! Compute local laplacian matrix (Used in poisson equation)
+  //! \param[in] grad_shapefn shape function gradient
+  //! \param[in] pvolume volume weight
+  //! \param[in] multiplier multiplier
+  void compute_local_laplacian(const Eigen::MatrixXd& grad_shapefn,
+                               double pvolume,
+                               double multiplier = 1.0) noexcept;
+
+  //! Return local laplacian RHS matrix
+  const Eigen::MatrixXd& poisson_right_matrix() {
+    return poisson_right_matrix_;
+  };
+
+  //! Compute local poisson RHS matrix (Used in poisson equation)
+  //! \param[in] shapefn shape function
+  //! \param[in] grad_shapefn shape function gradient
+  //! \param[in] pvolume volume weight
+  void compute_local_poisson_right(const Eigen::VectorXd& shapefn,
+                                   const Eigen::MatrixXd& grad_shapefn,
+                                   double pvolume,
+                                   double multiplier = 1.0) noexcept;
+
+  //! Return local correction matrix
+  const Eigen::MatrixXd& correction_matrix() { return correction_matrix_; };
+
+  //! Compute local correction matrix (Used to correct velocity)
+  void compute_local_correction_matrix(const Eigen::VectorXd& shapefn,
+                                       const Eigen::MatrixXd& grad_shapefn,
+                                       double pvolume) noexcept;
   //! Return previous mpi rank
   unsigned previous_mpirank() const;
 
@@ -264,6 +324,16 @@ class Cell {
   //! Normal of face
   //! first-> face_id, second->vector of the normal
   std::map<unsigned, Eigen::VectorXd> face_normals_;
+  //! Free surface bool
+  bool free_surface_{false};
+  //! Volume fraction
+  double volume_fraction_{0.0};
+  //! Local laplacian matrix
+  Eigen::MatrixXd laplacian_matrix_;
+  //! Local poisson RHS matrix
+  Eigen::MatrixXd poisson_right_matrix_;
+  //! Local correction RHS matrix
+  Eigen::MatrixXd correction_matrix_;
   //! Logger
   std::unique_ptr<spdlog::logger> console_;
 };  // Cell class

include/cell.tcc

@@ -412,7 +412,7 @@ inline Eigen::Matrix<double, 2, 1> mpm::Cell<2>::local_coordinates_point(
     if (indices.size() == 3) {
       //   2 0
       //     |\
-      //     | \  
+      //     | \
       //   c |  \ b
       //     |   \
       //     |    \
@@ -844,3 +844,107 @@ template <unsigned Tdim>
 inline unsigned mpm::Cell<Tdim>::previous_mpirank() const {
   return this->previous_mpirank_;
 }
+
+//! Map cell volume to nodes
+template <unsigned Tdim>
+bool mpm::Cell<Tdim>::map_cell_volume_to_nodes(unsigned phase) {
+  bool status = true;
+  try {
+    // Check if cell volume is set
+    if (volume_ == std::numeric_limits<double>::lowest())
+      this->compute_volume();
+
+    for (unsigned i = 0; i < nodes_.size(); ++i) {
+      nodes_[i]->update_volume(true, phase, volume_ / nnodes_);
+    }
+
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+    status = false;
+  }
+  return status;
+}
+
+//! Return local node indices
+template <unsigned Tdim>
+Eigen::VectorXi mpm::Cell<Tdim>::local_node_indices() {
+  Eigen::VectorXi indices;
+  try {
+    indices.resize(nodes_.size());
+    indices.setZero();
+    unsigned node_idx = 0;
+    for (auto node_itr = nodes_.cbegin(); node_itr != nodes_.cend();
+         ++node_itr) {
+      indices(node_idx) = (*node_itr)->active_id();
+      node_idx++;
+    }
+  } catch (std::exception& exception) {
+    console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+  }
+  return indices;
+}
+
+//! Initialise element matrix
+template <unsigned Tdim>
+bool mpm::Cell<Tdim>::initialise_element_matrix() {
+  bool status = true;
+  if (this->status()) {
+    try {
+      // Initialse Laplacian matrix (NxN)
+      laplacian_matrix_.resize(nnodes_, nnodes_);
+      laplacian_matrix_.setZero();
+
+      // Initialse poisson RHS matrix (Nx(N*Tdim))
+      poisson_right_matrix_.resize(nnodes_, nnodes_ * Tdim);
+      poisson_right_matrix_.setZero();
+
+      // Initialse correction RHS matrix (NxTdim)
+      correction_matrix_.resize(nnodes_, nnodes_ * Tdim);
+      correction_matrix_.setZero();
+
+    } catch (std::exception& exception) {
+      console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
+      status = false;
+    }
+  }
+  return status;
+}
+
+//! Compute local matrix of laplacian
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::compute_local_laplacian(
+    const Eigen::MatrixXd& grad_shapefn, double pvolume,
+    double multiplier) noexcept {
+
+  std::lock_guard<std::mutex> guard(cell_mutex_);
+  laplacian_matrix_ +=
+      grad_shapefn * grad_shapefn.transpose() * multiplier * pvolume;
+}
+
+//! Compute local poisson RHS matrix
+//! Used in poisson equation RHS for Navier Stokes solver
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::compute_local_poisson_right(
+    const Eigen::VectorXd& shapefn, const Eigen::MatrixXd& grad_shapefn,
+    double pvolume, double multiplier) noexcept {
+
+  std::lock_guard<std::mutex> guard(cell_mutex_);
+  for (unsigned i = 0; i < Tdim; i++) {
+    poisson_right_matrix_.block(0, i * nnodes_, nnodes_, nnodes_) +=
+        shapefn * grad_shapefn.col(i).transpose() * multiplier * pvolume;
+  }
+}
+
+//! Compute local correction matrix
+//! Used to compute corrector of nodal velocity for Navier Stokes solver
+template <unsigned Tdim>
+void mpm::Cell<Tdim>::compute_local_correction_matrix(
+    const Eigen::VectorXd& shapefn, const Eigen::MatrixXd& grad_shapefn,
+    double pvolume) noexcept {
+
+  std::lock_guard<std::mutex> guard(cell_mutex_);
+  for (unsigned i = 0; i < Tdim; i++) {
+    correction_matrix_.block(0, i * nnodes_, nnodes_, nnodes_) +=
+        shapefn * grad_shapefn.col(i).transpose() * pvolume;
+  }
+}

include/io/io_mesh.h

@@ -98,6 +98,11 @@ class IOMesh {
       read_friction_constraints(
           const std::string& friction_constraints_file) = 0;
 
+  //! Read pressure constraints file
+  //! \param[in] pressure_constraints_files file name with pressure constraints
+  virtual std::vector<std::tuple<mpm::Index, double>> read_pressure_constraints(
+      const std::string& _pressure_constraints_file) = 0;
+
   //! Read forces file
   //! \param[in] forces_files file name with nodal concentrated force
   virtual std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(

include/io/io_mesh_ascii.h

@@ -91,6 +91,12 @@ class IOMeshAscii : public IOMesh<Tdim> {
   std::vector<std::tuple<mpm::Index, unsigned, double>> read_forces(
       const std::string& forces_file) override;
 
+  //! Read pressure constraints file
+  //! \param[in] pressure_constraints_files file name with pressure
+  //! constraints
+  std::vector<std::tuple<mpm::Index, double>> read_pressure_constraints(
+      const std::string& pressure_constraints_file) override;
+
  private:
   //! Logger
   std::shared_ptr<spdlog::logger> console_;

include/io/io_mesh_ascii.tcc

@@ -486,6 +486,49 @@ std::vector<std::tuple<mpm::Index, unsigned, int, double>>
   return constraints;
 }
 
+//! Read pressure constraints file
+template <unsigned Tdim>
+std::vector<std::tuple<mpm::Index, double>>
+    mpm::IOMeshAscii<Tdim>::read_pressure_constraints(
+        const std::string& pressure_constraints_file) {
+  // Particle pressure constraints
+  std::vector<std::tuple<mpm::Index, double>> constraints;
+  constraints.clear();
+
+  // input file stream
+  std::fstream file;
+  file.open(pressure_constraints_file.c_str(), std::ios::in);
+
+  try {
+    if (file.is_open() && file.good()) {
+      // Line
+      std::string line;
+      while (std::getline(file, line)) {
+        boost::algorithm::trim(line);
+        std::istringstream istream(line);
+        // ignore comment lines (# or !) or blank lines
+        if ((line.find('#') == std::string::npos) &&
+            (line.find('!') == std::string::npos) && (line != "")) {
+          while (istream.good()) {
+            // ID
+            mpm::Index id;
+            // Pressure
+            double pressure;
+            // Read stream
+            istream >> id >> pressure;
+            constraints.emplace_back(std::make_tuple(id, pressure));
+          }
+        }
+      }
+    }
+    file.close();
+  } catch (std::exception& exception) {
+    console_->error("Read pressure constraints: {}", exception.what());
+    file.close();
+  }
+  return constraints;
+}
+
 //! Return particles force
 template <unsigned Tdim>
 std::vector<std::tuple<mpm::Index, unsigned, double>>

include/io/logger.h

@@ -41,6 +41,10 @@ struct Logger {
 
   // Create a logger for MPM Explicit USL
   static const std::shared_ptr<spdlog::logger> mpm_explicit_usl_logger;
+
+  // Create a logger for MPM Semi-implicit Navier Stokes
+  static const std::shared_ptr<spdlog::logger>
+      mpm_semi_implicit_navier_stokes_logger;
 };
 
 }  // namespace mpm