diff --git a/CMakeLists.txt b/CMakeLists.txt index 5fc7cfa1d..d196fd1a2 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/include/cell.h b/include/cell.h index 2380b7d3d..489affc33 100644 --- a/include/cell.h +++ b/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 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 console_; }; // Cell class diff --git a/include/cell.tcc b/include/cell.tcc index 05f839dab..9cd3f6401 100644 --- a/include/cell.tcc +++ b/include/cell.tcc @@ -412,7 +412,7 @@ inline Eigen::Matrix mpm::Cell<2>::local_coordinates_point( if (indices.size() == 3) { // 2 0 // |\ - // | \ + // | \ // c | \ b // | \ // | \ @@ -844,3 +844,107 @@ template inline unsigned mpm::Cell::previous_mpirank() const { return this->previous_mpirank_; } + +//! Map cell volume to nodes +template +bool mpm::Cell::map_cell_volume_to_nodes(unsigned phase) { + bool status = true; + try { + // Check if cell volume is set + if (volume_ == std::numeric_limits::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 +Eigen::VectorXi mpm::Cell::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 +bool mpm::Cell::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 +void mpm::Cell::compute_local_laplacian( + const Eigen::MatrixXd& grad_shapefn, double pvolume, + double multiplier) noexcept { + + std::lock_guard 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 +void mpm::Cell::compute_local_poisson_right( + const Eigen::VectorXd& shapefn, const Eigen::MatrixXd& grad_shapefn, + double pvolume, double multiplier) noexcept { + + std::lock_guard 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 +void mpm::Cell::compute_local_correction_matrix( + const Eigen::VectorXd& shapefn, const Eigen::MatrixXd& grad_shapefn, + double pvolume) noexcept { + + std::lock_guard 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; + } +} diff --git a/include/io/io_mesh.h b/include/io/io_mesh.h index 4150af5b0..f1c05af6d 100644 --- a/include/io/io_mesh.h +++ b/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> 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> read_forces( diff --git a/include/io/io_mesh_ascii.h b/include/io/io_mesh_ascii.h index 433661546..9a9cbf4e1 100644 --- a/include/io/io_mesh_ascii.h +++ b/include/io/io_mesh_ascii.h @@ -91,6 +91,12 @@ class IOMeshAscii : public IOMesh { std::vector> read_forces( const std::string& forces_file) override; + //! Read pressure constraints file + //! \param[in] pressure_constraints_files file name with pressure + //! constraints + std::vector> read_pressure_constraints( + const std::string& pressure_constraints_file) override; + private: //! Logger std::shared_ptr console_; diff --git a/include/io/io_mesh_ascii.tcc b/include/io/io_mesh_ascii.tcc index 0f09fdc1a..1059c6542 100644 --- a/include/io/io_mesh_ascii.tcc +++ b/include/io/io_mesh_ascii.tcc @@ -486,6 +486,49 @@ std::vector> return constraints; } +//! Read pressure constraints file +template +std::vector> + mpm::IOMeshAscii::read_pressure_constraints( + const std::string& pressure_constraints_file) { + // Particle pressure constraints + std::vector> 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 std::vector> diff --git a/include/io/logger.h b/include/io/logger.h index 082a0e874..1ecbce67f 100644 --- a/include/io/logger.h +++ b/include/io/logger.h @@ -41,6 +41,10 @@ struct Logger { // Create a logger for MPM Explicit USL static const std::shared_ptr mpm_explicit_usl_logger; + + // Create a logger for MPM Semi-implicit Navier Stokes + static const std::shared_ptr + mpm_semi_implicit_navier_stokes_logger; }; } // namespace mpm diff --git a/include/linear_solvers/assembler_base.h b/include/linear_solvers/assembler_base.h new file mode 100644 index 000000000..859acbfd5 --- /dev/null +++ b/include/linear_solvers/assembler_base.h @@ -0,0 +1,94 @@ +#ifndef MPM_ASSEMBLER_BASE_H_ +#define MPM_ASSEMBLER_BASE_H_ + +#include +#include +#include +#include +#include + +#include + +#include "mesh.h" +#include "node_base.h" + +namespace mpm { + +// Linear solver assembler base class +//! \brief Perform matrix assembly procedures +//! \details Build global matrices considering local element matrices +//! \tparam Tdim Dimension +template +class AssemblerBase { + public: + AssemblerBase() { + //! Global degrees of freedom + active_dof_ = 0; + } + + // Virtual destructor + virtual ~AssemblerBase() = default; + + //! Copy constructor + AssemblerBase(const AssemblerBase&) = default; + + //! Assignment operator + AssemblerBase& operator=(const AssemblerBase&) = default; + + //! Move constructor + AssemblerBase(AssemblerBase&&) = default; + + //! Assign mesh pointer + //! \param[in] mesh mesh pointer + void assign_mesh_pointer(const std::shared_ptr>& mesh) { + mesh_ = mesh; + } + + //! Create a pair between nodes and index in Matrix / Vector + virtual bool assign_global_node_indices(unsigned active_dof) = 0; + + //! Assemble laplacian matrix + virtual Eigen::SparseMatrix& laplacian_matrix() = 0; + + //! Assemble laplacian matrix + virtual bool assemble_laplacian_matrix(double dt) = 0; + + //! Assemble poisson RHS vector + virtual Eigen::VectorXd& poisson_rhs_vector() = 0; + + //! Assemble poisson RHS vector + virtual bool assemble_poisson_right(double dt) = 0; + + //! Assign free surface node id + virtual void assign_free_surface( + const std::set& free_surface_id) = 0; + + //! Assign pressure constraints + virtual bool assign_pressure_constraints(double beta, + double current_time) = 0; + + //! Apply pressure constraints to poisson equation + virtual void apply_pressure_constraints() = 0; + + //! Return pressure increment + virtual Eigen::VectorXd& pressure_increment() = 0; + + //! Assign pressure increment + virtual void assign_pressure_increment( + const Eigen::VectorXd& pressure_increment) = 0; + + //! Return correction matrix + virtual Eigen::SparseMatrix& correction_matrix() = 0; + + //! Assemble corrector RHS + virtual bool assemble_corrector_right(double dt) = 0; + + protected: + //! Number of total active_dof + unsigned active_dof_; + //! Mesh object + std::shared_ptr> mesh_; +}; +} // namespace mpm + +#endif // MPM_ASSEMBLER_BASE_H_ \ No newline at end of file diff --git a/include/linear_solvers/assembler_eigen_semi_implicit_navierstokes.h b/include/linear_solvers/assembler_eigen_semi_implicit_navierstokes.h new file mode 100644 index 000000000..2679b7da5 --- /dev/null +++ b/include/linear_solvers/assembler_eigen_semi_implicit_navierstokes.h @@ -0,0 +1,94 @@ +#ifndef MPM_ASSEMBLER_EIGEN_SEMI_IMPLICIT_NAVIERSTOKES_H_ +#define MPM_ASSEMBLER_EIGEN_SEMI_IMPLICIT_NAVIERSTOKES_H_ + +#include +#include + +// Speed log +#include "assembler_base.h" +#include "spdlog/spdlog.h" + +#include "cg_eigen.h" +#include "mesh.h" + +namespace mpm { +template +class AssemblerEigenSemiImplicitNavierStokes : public AssemblerBase { + public: + //! Constructor + AssemblerEigenSemiImplicitNavierStokes(); + + //! Create a pair between nodes and index in Matrix / Vector + bool assign_global_node_indices(unsigned active_dof) override; + + //! Return laplacian matrix + Eigen::SparseMatrix& laplacian_matrix() override { + return laplacian_matrix_; + } + + //! Assemble laplacian matrix + bool assemble_laplacian_matrix(double dt) override; + + //! Return poisson RHS vector + Eigen::VectorXd& poisson_rhs_vector() override { return poisson_rhs_vector_; } + + //! Assemble poisson RHS vector + bool assemble_poisson_right(double dt) override; + + //! Assign free surface node id + void assign_free_surface( + const std::set& free_surface_id) override { + free_surface_ = free_surface_id; + } + + //! Assign pressure constraints + bool assign_pressure_constraints(double beta, double current_time) override; + + //! Apply pressure constraints to poisson equation + void apply_pressure_constraints(); + + //! Return pressure increment + Eigen::VectorXd& pressure_increment() override { return pressure_increment_; } + + //! Assign pressure increment + void assign_pressure_increment( + const Eigen::VectorXd& pressure_increment) override { + pressure_increment_ = pressure_increment; + } + + //! Return correction matrix + Eigen::SparseMatrix& correction_matrix() override { + return correction_matrix_; + } + + //! Assemble corrector RHS + bool assemble_corrector_right(double dt) override; + + protected: + //! Logger + std::shared_ptr console_; + + private: + //! number of nodes + using AssemblerBase::active_dof_; + //! Mesh object + using AssemblerBase::mesh_; + //! Laplacian matrix + Eigen::SparseMatrix laplacian_matrix_; + //! Poisson RHS vector + Eigen::VectorXd poisson_rhs_vector_; + //! Free surface + std::set free_surface_; + //! Pressure constraints + Eigen::SparseVector pressure_constraints_; + //! \delta p^(t+1) = p^(t+1) - beta * p^(t) + Eigen::VectorXd pressure_increment_; + //! correction_matrix + Eigen::SparseMatrix correction_matrix_; + //! Global node indices + std::vector global_node_indices_; +}; +} // namespace mpm + +#include "assembler_eigen_semi_implicit_navierstokes.tcc" +#endif // MPM_ASSEMBLER_EIGEN_SEMI_IMPLICIT_NAVIERSTOKES_H_ diff --git a/include/linear_solvers/assembler_eigen_semi_implicit_navierstokes.tcc b/include/linear_solvers/assembler_eigen_semi_implicit_navierstokes.tcc new file mode 100644 index 000000000..f913109bd --- /dev/null +++ b/include/linear_solvers/assembler_eigen_semi_implicit_navierstokes.tcc @@ -0,0 +1,301 @@ +//! Construct a semi-implicit eigen matrix assembler +template +mpm::AssemblerEigenSemiImplicitNavierStokes< + Tdim>::AssemblerEigenSemiImplicitNavierStokes() + : mpm::AssemblerBase() { + //! Logger + console_ = spdlog::stdout_color_mt("AssemblerEigenSemiImplicitNavierStokes"); +} + +//! Assign global node indices +template +bool mpm::AssemblerEigenSemiImplicitNavierStokes< + Tdim>::assign_global_node_indices(unsigned nactive_node) { + bool status = true; + try { + // Total number of active node + active_dof_ = nactive_node; + + global_node_indices_ = mesh_->global_node_indices(); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Assemble Laplacian matrix +template +bool mpm::AssemblerEigenSemiImplicitNavierStokes< + Tdim>::assemble_laplacian_matrix(double dt) { + bool status = true; + try { + // Initialise Laplacian matrix + laplacian_matrix_.resize(active_dof_, active_dof_); + laplacian_matrix_.setZero(); + + // Reserve storage for sparse matrix + switch (Tdim) { + // For 2d: 10 entries /column + case (2): { + laplacian_matrix_.reserve(Eigen::VectorXi::Constant(active_dof_, 10)); + break; + } + // For 3d: 30 entries /column + case (3): { + laplacian_matrix_.reserve(Eigen::VectorXi::Constant(active_dof_, 30)); + break; + } + } + + // Cell pointer + const auto& cells = mesh_->cells(); + + // Iterate over cells + mpm::Index cid = 0; + for (auto cell_itr = cells.cbegin(); cell_itr != cells.cend(); ++cell_itr) { + if ((*cell_itr)->status()) { + // Node ids in each cell + const auto nids = global_node_indices_.at(cid); + + // Laplacian element of cell + const auto cell_laplacian = (*cell_itr)->laplacian_matrix(); + + // Assemble global laplacian matrix + for (unsigned i = 0; i < nids.size(); ++i) { + for (unsigned j = 0; j < nids.size(); ++j) { + laplacian_matrix_.coeffRef(global_node_indices_.at(cid)(i), + global_node_indices_.at(cid)(j)) += + cell_laplacian(i, j); + } + } + + ++cid; + } + } + + laplacian_matrix_ *= dt; + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +template +bool mpm::AssemblerEigenSemiImplicitNavierStokes::assemble_poisson_right( + double dt) { + bool status = true; + try { + // Initialise Poisson RHS matrix + Eigen::SparseMatrix poisson_right_matrix; + poisson_right_matrix.resize(active_dof_, active_dof_ * Tdim); + poisson_right_matrix.setZero(); + + // Reserve storage for sparse matrix + switch (Tdim) { + // For 2d: 10 entries /column + case (2): { + poisson_right_matrix.reserve( + Eigen::VectorXi::Constant(active_dof_ * Tdim, 10)); + break; + } + // For 3d: 30 entries /column + case (3): { + poisson_right_matrix.reserve( + Eigen::VectorXi::Constant(active_dof_ * Tdim, 30)); + break; + } + } + + // Cell pointer + const auto& cells = mesh_->cells(); + + // Iterate over cells + mpm::Index cid = 0; + for (auto cell_itr = cells.cbegin(); cell_itr != cells.cend(); ++cell_itr) { + if ((*cell_itr)->status()) { + // Node ids in each cell + const auto nids = global_node_indices_.at(cid); + + // Local Poisson RHS matrix + auto cell_poisson_right = (*cell_itr)->poisson_right_matrix(); + + // Assemble global poisson RHS matrix + for (unsigned i = 0; i < nids.size(); ++i) { + for (unsigned j = 0; j < nids.size(); ++j) { + for (unsigned k = 0; k < Tdim; ++k) { + poisson_right_matrix.coeffRef( + global_node_indices_.at(cid)(i), + global_node_indices_.at(cid)(j) + k * active_dof_) += + cell_poisson_right(i, j + k * nids.size()); + } + } + } + cid++; + } + } + + // Resize poisson right vector + poisson_rhs_vector_.resize(active_dof_); + poisson_rhs_vector_.setZero(); + + // Compute velocity + Eigen::MatrixXd fluid_velocity; + fluid_velocity.resize(active_dof_, Tdim); + fluid_velocity.setZero(); + + // Active nodes + const auto& active_nodes = mesh_->active_nodes(); + const unsigned fluid = mpm::ParticlePhase::SinglePhase; + unsigned node_index = 0; + + for (auto node_itr = active_nodes.cbegin(); node_itr != active_nodes.cend(); + ++node_itr) { + // Compute nodal intermediate force + fluid_velocity.row(node_index) = (*node_itr)->velocity(fluid).transpose(); + node_index++; + } + + // Resize fluid velocity + fluid_velocity.resize(active_dof_ * Tdim, 1); + + // Compute poisson RHS vector + poisson_rhs_vector_ = -poisson_right_matrix * fluid_velocity; + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Assemble corrector right matrix +template +bool mpm::AssemblerEigenSemiImplicitNavierStokes< + Tdim>::assemble_corrector_right(double dt) { + bool status = true; + try { + // Resize correction matrix + correction_matrix_.resize(active_dof_, active_dof_ * Tdim); + correction_matrix_.setZero(); + + // Reserve storage for sparse matrix + switch (Tdim) { + // For 2d: 10 entries /column + case (2): { + correction_matrix_.reserve( + Eigen::VectorXi::Constant(active_dof_ * Tdim, 10)); + break; + } + // For 3d: 30 entries /column + case (3): { + correction_matrix_.reserve( + Eigen::VectorXi::Constant(active_dof_ * Tdim, 30)); + break; + } + } + + // Cell pointer + unsigned nnodes_per_cell = global_node_indices_.at(0).size(); + const auto& cells = mesh_->cells(); + + // Iterate over cells + unsigned cid = 0; + for (auto cell_itr = cells.cbegin(); cell_itr != cells.cend(); ++cell_itr) { + if ((*cell_itr)->status()) { + auto cell_correction_matrix = (*cell_itr)->correction_matrix(); + for (unsigned k = 0; k < Tdim; k++) { + for (unsigned i = 0; i < nnodes_per_cell; i++) { + for (unsigned j = 0; j < nnodes_per_cell; j++) { + // Fluid + correction_matrix_.coeffRef( + global_node_indices_.at(cid)(i), + k * active_dof_ + global_node_indices_.at(cid)(j)) += + cell_correction_matrix(i, j + k * nnodes_per_cell); + } + } + } + cid++; + } + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Apply pressure constraints vector +template +void mpm::AssemblerEigenSemiImplicitNavierStokes< + Tdim>::apply_pressure_constraints() { + try { + // Modify poisson_rhs_vector_ + poisson_rhs_vector_ -= laplacian_matrix_ * pressure_constraints_; + + // Apply free surface + if (!free_surface_.empty()) { + const auto& nodes = mesh_->nodes(); + for (const auto& free_node : free_surface_) { + const auto column_index = nodes[free_node]->active_id(); + // Modify poisson_rhs_vector + poisson_rhs_vector_(column_index) = 0; + // Modify laplacian_matrix + laplacian_matrix_.row(column_index) *= 0; + laplacian_matrix_.col(column_index) *= 0; + laplacian_matrix_.coeffRef(column_index, column_index) = 1; + } + // Clear the vector + free_surface_.clear(); + } + + // Apply pressure constraints + for (Eigen::SparseVector::InnerIterator it(pressure_constraints_); + it; ++it) { + // Modify poisson_rhs_vector + poisson_rhs_vector_(it.index()) = it.value(); + // Modify laplacian_matrix + laplacian_matrix_.row(it.index()) *= 0; + laplacian_matrix_.col(it.index()) *= 0; + laplacian_matrix_.coeffRef(it.index(), it.index()) = 1; + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } +} + +//! Assign pressure constraints +template +bool mpm::AssemblerEigenSemiImplicitNavierStokes< + Tdim>::assign_pressure_constraints(double beta, double current_time) { + bool status = false; + try { + // Resize pressure constraints vector + pressure_constraints_.resize(active_dof_); + pressure_constraints_.reserve(int(0.5 * active_dof_)); + + // Nodes container + const auto& nodes = mesh_->active_nodes(); + // Iterate over nodes to get pressure constraints + for (auto node = nodes.cbegin(); node != nodes.cend(); ++node) { + // Assign total pressure constraint + const double pressure_constraint = + (*node)->pressure_constraint(0, current_time); + + // Check if there is a pressure constraint + if (pressure_constraint != std::numeric_limits::max()) { + // Insert the pressure constraints + pressure_constraints_.insert((*node)->active_id()) = + pressure_constraint; + } + } + status = true; + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return status; +} diff --git a/include/linear_solvers/cg_eigen.h b/include/linear_solvers/cg_eigen.h new file mode 100644 index 000000000..359df75b7 --- /dev/null +++ b/include/linear_solvers/cg_eigen.h @@ -0,0 +1,48 @@ +#ifndef MPM_CG_EIGEN_H_ +#define MPM_CG_EIGEN_H_ + +#include + +#include "factory.h" +#include "solver_base.h" +#include +#include +#include + +namespace mpm { + +//! MPM Eigen CG class +//! \brief Conjugate Gradient solver class using Eigen +template +class CGEigen : public SolverBase { + public: + //! Constructor + //! \param[in] max_iter Maximum number of iterations + //! \param[in] tolerance Tolerance for solver to achieve convergence + CGEigen(unsigned max_iter, double tolerance) + : mpm::SolverBase(max_iter, tolerance) { + //! Logger + console_ = spdlog::stdout_color_mt("EigenSolver"); + }; + + //! Matrix solver with default initial guess + Eigen::VectorXd solve(const Eigen::SparseMatrix& A, + const Eigen::VectorXd& b, + std::string solver_type) override; + + //! Return the type of solver + std::string solver_type() const { return "Eigen"; } + + protected: + //! Maximum number of iterations + using SolverBase::max_iter_; + //! Tolerance + using SolverBase::tolerance_; + //! Logger + using SolverBase::console_; +}; +} // namespace mpm + +#include "cg_eigen.tcc" + +#endif // MPM_CG_EIGEN_H_ diff --git a/include/linear_solvers/cg_eigen.tcc b/include/linear_solvers/cg_eigen.tcc new file mode 100644 index 000000000..783b247fb --- /dev/null +++ b/include/linear_solvers/cg_eigen.tcc @@ -0,0 +1,43 @@ +//! Conjugate Gradient with default initial guess +template +Eigen::VectorXd mpm::CGEigen::solve( + const Eigen::SparseMatrix& A, const Eigen::VectorXd& b, + std::string solver_type) { + Eigen::VectorXd x; + try { + + if (solver_type == "cg") { + Eigen::ConjugateGradient> solver; + + solver.setMaxIterations(max_iter_); + solver.setTolerance(tolerance_); + solver.compute(A); + + x = solver.solve(b); + + if (solver.info() != Eigen::Success) { + throw std::runtime_error("Fail to solve linear systems!\n"); + } + + } else if (solver_type == "lscg") { + + // Another option is LDLT, but not accurate as lscg + // Eigen::SimplicialLDLT> solver; + Eigen::LeastSquaresConjugateGradient> solver; + + solver.setMaxIterations(max_iter_); + solver.setTolerance(tolerance_); + solver.compute(A); + + x = solver.solve(b); + + if (solver.info() != Eigen::Success) { + throw std::runtime_error("Fail to solve linear systems!\n"); + } + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return x; +} \ No newline at end of file diff --git a/include/linear_solvers/solver_base.h b/include/linear_solvers/solver_base.h new file mode 100644 index 000000000..111b2cb08 --- /dev/null +++ b/include/linear_solvers/solver_base.h @@ -0,0 +1,32 @@ +#ifndef MPM_SOLVER_BASE_H_ +#define MPM_SOLVER_BASE_H_ + +#include "spdlog/spdlog.h" +namespace mpm { +template +class SolverBase { + public: + //! Constructor with min and max iterations and tolerance + SolverBase(unsigned max_iter, double tolerance) { + //! Logger + max_iter_ = max_iter; + tolerance_ = tolerance; + console_ = spdlog::stdout_color_mt("SolverBase"); + }; + + //! Matrix solver with default initial guess + virtual Eigen::VectorXd solve(const Eigen::SparseMatrix& A, + const Eigen::VectorXd& b, + std::string solver_type) = 0; + + protected: + //! Maximum number of iterations + unsigned max_iter_; + //! Tolerance + double tolerance_; + //! Logger + std::shared_ptr console_; +}; +} // namespace mpm + +#endif \ No newline at end of file diff --git a/include/loads_bcs/constraints.h b/include/loads_bcs/constraints.h index 0efdef6f6..df5f62ac9 100644 --- a/include/loads_bcs/constraints.h +++ b/include/loads_bcs/constraints.h @@ -6,6 +6,7 @@ #include "friction_constraint.h" #include "logger.h" #include "mesh.h" +#include "pressure_constraint.h" #include "velocity_constraint.h" namespace mpm { @@ -48,6 +49,21 @@ class Constraints { const std::vector>& friction_constraints); + //! Assign nodal pressure constraints + //! \param[in] mfunction Math function if defined + //! \param[in] setid Node set id + //! \param[in] phase Index corresponding to the phase + //! \param[in] pconstraint Pressure constraint at node + bool assign_nodal_pressure_constraint( + const std::shared_ptr& mfunction, + const std::shared_ptr& pconstraint); + + //! Assign nodal pressure constraints to nodes + //! \param[in] pressure_constraints Constraint at node, pressure + bool assign_nodal_pressure_constraints( + const unsigned phase, + const std::vector>& pressure_constraints); + private: //! Mesh object std::shared_ptr> mesh_; diff --git a/include/loads_bcs/constraints.tcc b/include/loads_bcs/constraints.tcc index 5813e609f..7126d6f92 100644 --- a/include/loads_bcs/constraints.tcc +++ b/include/loads_bcs/constraints.tcc @@ -105,3 +105,59 @@ bool mpm::Constraints::assign_nodal_friction_constraints( } return status; } + +//! Assign nodal pressure constraints +template +bool mpm::Constraints::assign_nodal_pressure_constraint( + const std::shared_ptr& mfunction, + const std::shared_ptr& pconstraint) { + bool status = true; + try { + int set_id = pconstraint->setid(); + auto nset = mesh_->nodes(set_id); + if (nset.size() != 0) { + unsigned phase = pconstraint->phase(); + double pressure = pconstraint->pressure(); + + for (auto nitr = nset.cbegin(); nitr != nset.cend(); ++nitr) { + if (!(*nitr)->assign_pressure_constraint(phase, pressure, mfunction)) + throw std::runtime_error( + "Failed to initialise pressure constraint at node"); + } + } else + throw std::runtime_error( + "Nodal assignment of pressure constraint failed"); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Assign nodal pressure constraints to nodes +template +bool mpm::Constraints::assign_nodal_pressure_constraints( + const unsigned phase, + const std::vector>& pressure_constraints) { + bool status = false; + try { + for (const auto& pressure_constraint : pressure_constraints) { + // Node id + mpm::Index nid = std::get<0>(pressure_constraint); + // Pressure + double pressure = std::get<1>(pressure_constraint); + + // Apply constraint + status = mesh_->node(nid)->assign_pressure_constraint(phase, pressure, + nullptr); + + if (!status) + throw std::runtime_error("Node or pressure constraint is invalid"); + } + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} diff --git a/include/loads_bcs/pressure_constraint.h b/include/loads_bcs/pressure_constraint.h new file mode 100644 index 000000000..4ea68bb06 --- /dev/null +++ b/include/loads_bcs/pressure_constraint.h @@ -0,0 +1,35 @@ +#ifndef MPM_PRESSURE_CONSTRAINT_H_ +#define MPM_PRESSURE_CONSTRAINT_H_ + +namespace mpm { + +//! PressureConstraint class to store pressure constraint on a set +//! \brief PressureConstraint class to store a constraint on a set +//! \details PressureConstraint stores the constraint as a static value +class PressureConstraint { + public: + // Constructor + //! \param[in] setid set id + //! \param[in] pressure Constraint pressure + PressureConstraint(int setid, unsigned phase, double pressure) + : setid_{setid}, phase_{phase}, pressure_{pressure} {}; + + // Set id + int setid() const { return setid_; } + + // Return phase + unsigned phase() const { return phase_; } + + // Return pressure + double pressure() const { return pressure_; } + + private: + // ID + int setid_; + // Phase + unsigned phase_; + // Velocity + double pressure_; +}; +} // namespace mpm +#endif // MPM_PRESSURE_CONSTRAINT_H_ diff --git a/include/mesh.h b/include/mesh.h index 2de4f49e9..6c16e3dc4 100644 --- a/include/mesh.h +++ b/include/mesh.h @@ -10,6 +10,7 @@ // Eigen #include "Eigen/Dense" +#include // MPI #ifdef USE_MPI #include "mpi.h" @@ -40,6 +41,7 @@ using Json = nlohmann::json; #include "node.h" #include "particle.h" #include "particle_base.h" +#include "radial_basis_function.h" #include "traction.h" #include "vector.h" #include "velocity_constraint.h" @@ -101,6 +103,9 @@ class Mesh { //! Return the number of nodes mpm::Index nnodes() const { return nodes_.size(); } + //! Return container of nodes + mpm::Vector> nodes() { return nodes_; } + //! Return the number of nodes in rank mpm::Index nnodes_rank(); @@ -249,6 +254,12 @@ class Mesh { template void iterate_over_particles(Toper oper); + //! Iterate over particles with predicate + //! \tparam Toper Callable object typically a baseclass functor + //! \tparam Tpred Predicate + template + void iterate_over_particles_predicate(Toper oper, Tpred pred); + //! Iterate over particle set //! \tparam Toper Callable object typically a baseclass functor //! \param[in] set_id particle set id @@ -441,6 +452,58 @@ class Mesh { //! Inject particles void inject_particles(double current_time); + //! Compute free surface + //! \param[in] method Type of method to use + //! \param[in] volume_tolerance for volume_fraction approach + //! \retval status Status of compute_free_surface + bool compute_free_surface( + std::string method, + double volume_tolerance = std::numeric_limits::epsilon()); + + //! Compute free surface by density method + //! \details Using simple approach of volume fraction approach as (Kularathna + //! & Soga, 2017) and density ratio comparison (Hamad, 2015). This method is + //! fast, but less accurate. + //! \param[in] volume_tolerance for volume_fraction approach + //! \retval status Status of compute_free_surface + bool compute_free_surface_by_density( + double volume_tolerance = std::numeric_limits::epsilon()); + + //! Compute free surface by geometry method + //! \details Using a more expensive approach using neighbouring particles and + //! current geometry. This method combine multiple checks in order to simplify + //! and fasten the process: (1) Volume fraction approach as (Kularathna & Soga + //! 2017), (2) Density comparison approach as (Hamad, 2015), and (3) Geometry + //! based approach as (Marrone et al. 2010) + //! \param[in] volume_tolerance for volume_fraction approach + //! \retval status Status of compute_free_surface + bool compute_free_surface_by_geometry( + double volume_tolerance = std::numeric_limits::epsilon()); + + //! Get free surface node set + std::set free_surface_nodes(); + + //! Get free surface cell set + std::set free_surface_cells(); + + //! Get free surface particle set + std::set free_surface_particles(); + + //! Create a list of active nodes in mesh and assign active node id + //! (rank-wise) + unsigned assign_active_nodes_id(); + + //! Return container of active nodes + mpm::Vector> active_nodes() { return active_nodes_; } + + //! Return global node indices + std::vector global_node_indices() const; + + //! Compute correction force in the node + bool compute_nodal_correction_force( + const Eigen::SparseMatrix& correction_matrix, + const Eigen::VectorXd& pressure_increment, double dt); + // Create the nodal properties' map void create_nodal_properties(); diff --git a/include/mesh.tcc b/include/mesh.tcc index c603a2fea..f090d0477 100644 --- a/include/mesh.tcc +++ b/include/mesh.tcc @@ -93,6 +93,45 @@ void mpm::Mesh::find_active_nodes() { if ((*nitr)->status()) this->active_nodes_.add(*nitr); } +//! Create a list of active nodes in mesh and assign active node id (rank-wise) +template +unsigned mpm::Mesh::assign_active_nodes_id() { + // Clear existing list of active nodes + this->active_nodes_.clear(); + Index active_id = 0; + + for (auto nitr = nodes_.cbegin(); nitr != nodes_.cend(); ++nitr) { + if ((*nitr)->status()) { + this->active_nodes_.add(*nitr); + (*nitr)->assign_active_id(active_id); + active_id++; + } else { + (*nitr)->assign_active_id(std::numeric_limits::max()); + } + } + + return active_id; +} + +//! Return global node indices +template +std::vector mpm::Mesh::global_node_indices() const { + // Vector of node_pairs + std::vector node_indices; + try { + // Iterate over cells + for (auto citr = cells_.cbegin(); citr != cells_.cend(); ++citr) { + if ((*citr)->status()) { + node_indices.emplace_back((*citr)->local_node_indices()); + } + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return node_indices; +} + //! Iterate over active nodes template template @@ -1022,6 +1061,16 @@ void mpm::Mesh::iterate_over_particles(Toper oper) { oper(*pitr); } +//! Iterate over particles +template +template +void mpm::Mesh::iterate_over_particles_predicate(Toper oper, Tpred pred) { +#pragma omp parallel for schedule(runtime) + for (auto pitr = particles_.cbegin(); pitr != particles_.cend(); ++pitr) { + if (pred(*pitr)) oper(*pitr); + } +} + //! Iterate over particle set template template @@ -1885,6 +1934,461 @@ void mpm::Mesh::create_nodal_properties() { } } +//! Compute nodal correction force +template +bool mpm::Mesh::compute_nodal_correction_force( + const Eigen::SparseMatrix& correction_matrix, + const Eigen::VectorXd& pressure_increment, double dt) { + bool status = true; + try { + //! Active node size + const auto nactive_node = active_nodes_.size(); + + // Part of nodal correction force of one direction + Eigen::MatrixXd correction_force; + correction_force.resize(nactive_node, Tdim); + + // Iterate over each direction + for (unsigned i = 0; i < Tdim; ++i) { + correction_force.block(0, i, nactive_node, 1) = + -correction_matrix.block(0, nactive_node * i, nactive_node, + nactive_node) * + pressure_increment; + } + + // Iterate over each active node + VectorDim nodal_correction_force; + for (auto nitr = active_nodes_.cbegin(); nitr != active_nodes_.cend(); + ++nitr) { + unsigned active_id = (*nitr)->active_id(); + nodal_correction_force = (correction_force.row(active_id)).transpose(); + + // Compute correction force for each node + map_nodes_[(*nitr)->id()]->compute_nodal_correction_force( + nodal_correction_force); + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return status; +}; + +//! Compute free surface cells, nodes, and particles +template +bool mpm::Mesh::compute_free_surface(std::string method, + double volume_tolerance) { + if (method == "density") { + return this->compute_free_surface_by_density(volume_tolerance); + } else if (method == "geometry") { + return this->compute_free_surface_by_geometry(volume_tolerance); + } else { + console_->info( + "The selected free-surface computation method: {}\n is not available. " + "Using density approach as default method.", + method); + return this->compute_free_surface_by_density(volume_tolerance); + } +} + +//! Compute free surface cells, nodes, and particles by density and geometry +template +bool mpm::Mesh::compute_free_surface_by_geometry( + double volume_tolerance) { + bool status = true; + try { + // Reset free surface cell and particles + this->iterate_over_cells(std::bind(&mpm::Cell::assign_free_surface, + std::placeholders::_1, false)); + + VectorDim temp_normal; + temp_normal.setZero(); + this->iterate_over_particles_predicate( + std::bind(&mpm::ParticleBase::assign_normal, + std::placeholders::_1, temp_normal), + std::bind(&mpm::ParticleBase::free_surface, + std::placeholders::_1)); + + this->iterate_over_particles( + std::bind(&mpm::ParticleBase::assign_free_surface, + std::placeholders::_1, false)); + + // Reset volume fraction + this->iterate_over_cells(std::bind(&mpm::Cell::assign_volume_fraction, + std::placeholders::_1, 0.0)); + + // Compute and assign volume fraction to each cell + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + if ((*citr)->status()) { + // Compute volume fraction + double cell_volume_fraction = 0.0; + for (const auto p_id : (*citr)->particles()) + cell_volume_fraction += map_particles_[p_id]->volume(); + + cell_volume_fraction = cell_volume_fraction / (*citr)->volume(); + (*citr)->assign_volume_fraction(cell_volume_fraction); + } + } + + // First, we detect the cell with possible free surfaces + // Compute boundary cells and nodes based on geometry + std::set free_surface_candidate_cells; + std::set free_surface_candidate_nodes; + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + // Cell contains particles + if ((*citr)->status()) { + bool candidate_cell = false; + const auto& node_id = (*citr)->nodes_id(); + if ((*citr)->volume_fraction() < volume_tolerance) { + candidate_cell = true; + for (const auto id : node_id) { + map_nodes_[id]->assign_free_surface(true); + } + } else { + // Loop over neighbouring cells + for (const auto n_id : (*citr)->neighbours()) { + if (!map_cells_[n_id]->status()) { + candidate_cell = true; + const auto& n_node_id = map_cells_[n_id]->nodes_id(); + + // Detect common node id + std::set common_node_id; + std::set_intersection( + node_id.begin(), node_id.end(), n_node_id.begin(), + n_node_id.end(), + std::inserter(common_node_id, common_node_id.begin())); + + // Assign free surface nodes + if (!common_node_id.empty()) { + for (const auto common_id : common_node_id) { + map_nodes_[common_id]->assign_free_surface(true); + } + } + } + } + } + + // Assign free surface cell + if (candidate_cell) { + (*citr)->assign_free_surface(true); + free_surface_candidate_cells.insert((*citr)->id()); + free_surface_candidate_nodes.insert(node_id.begin(), node_id.end()); + } + } + } + + // Compute particle neighbours for particles at candidate cells + std::vector free_surface_candidate_particles; + for (const auto cell_id : free_surface_candidate_cells) { + this->find_particle_neighbours(map_cells_[cell_id]); + const auto& particle_ids = map_cells_[cell_id]->particles(); + free_surface_candidate_particles.insert( + free_surface_candidate_particles.end(), particle_ids.begin(), + particle_ids.end()); + } + + // Compute boundary particles based on density function + // Lump cell volume to nodes + this->iterate_over_cells(std::bind( + &mpm::Cell::map_cell_volume_to_nodes, std::placeholders::_1, 0)); + + // Compute nodal value of mass density + this->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::compute_density, std::placeholders::_1), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + std::set free_surface_candidate_particles_second; + for (const auto p_id : free_surface_candidate_particles) { + const auto& particle = map_particles_[p_id]; + bool status = particle->compute_free_surface_by_density(); + if (status) free_surface_candidate_particles_second.insert(p_id); + } + + // Find free surface particles through geometry + std::set free_surface_particles; + for (const auto p_id : free_surface_candidate_particles_second) { + // Initialize renormalization matrix + Eigen::Matrix renormalization_matrix_inv; + renormalization_matrix_inv.setZero(); + + // Loop over neighbours + const auto& particle = map_particles_[p_id]; + const auto& p_coord = particle->coordinates(); + const auto& neighbour_particles = particle->neighbours(); + const double smoothing_length = 1.33 * particle->diameter(); + for (const auto n_id : neighbour_particles) { + const auto& n_coord = map_particles_[n_id]->coordinates(); + const VectorDim rel_coord = n_coord - p_coord; + + // Compute kernel gradient + const VectorDim kernel_gradient = + mpm::RadialBasisFunction::gradient(smoothing_length, + -rel_coord, "gaussian"); + + // Inverse of renormalization matrix B + renormalization_matrix_inv += + (particle->mass() / particle->mass_density()) * kernel_gradient * + rel_coord.transpose(); + } + + // Compute lambda: minimum eigenvalue of B_inverse + Eigen::SelfAdjointEigenSolver es( + renormalization_matrix_inv); + double lambda = es.eigenvalues().minCoeff(); + + // Categorize particle based on lambda + bool free_surface = false; + bool secondary_check = false; + bool interior = false; + if (lambda <= 0.2) + free_surface = true; + else if (lambda > 0.2 && lambda <= 0.75) + secondary_check = true; + else + interior = true; + + // Compute numerical normal vector + VectorDim normal; + normal.setZero(); + if (!interior) { + VectorDim temporary_vec; + temporary_vec.setZero(); + for (const auto n_id : neighbour_particles) { + const auto& n_coord = map_particles_[n_id]->coordinates(); + const VectorDim rel_coord = n_coord - p_coord; + + // Compute kernel gradient + const VectorDim kernel_gradient = + mpm::RadialBasisFunction::gradient(smoothing_length, + -rel_coord, "gaussian"); + + // Sum of kernel by volume + temporary_vec += + (particle->mass() / particle->mass_density()) * kernel_gradient; + } + normal = -renormalization_matrix_inv.inverse() * temporary_vec; + if (normal.norm() > std::numeric_limits::epsilon()) + normal.normalize(); + else + normal.setZero(); + } + + // If secondary check is needed + if (secondary_check) { + // Construct scanning region + // TODO: spacing distance should be a function of porosity + const double spacing_distance = smoothing_length; + VectorDim t_coord = p_coord + spacing_distance * normal; + + // Check all neighbours + for (const auto n_id : neighbour_particles) { + const auto& n_coord = map_particles_[n_id]->coordinates(); + const VectorDim rel_coord_np = n_coord - p_coord; + const double distance_np = rel_coord_np.norm(); + const VectorDim rel_coord_nt = n_coord - t_coord; + const double distance_nt = rel_coord_nt.norm(); + + free_surface = true; + if (distance_np < std::sqrt(2) * spacing_distance) { + if (std::acos(normal.dot(rel_coord_np) / distance_np) < M_PI / 4) { + free_surface = false; + break; + } + } else { + if (distance_nt < spacing_distance) { + free_surface = false; + break; + } + } + } + } + + // Assign normal only to validated free surface + if (free_surface) { + particle->assign_free_surface(true); + particle->assign_normal(normal); + free_surface_particles.insert(p_id); + } + } + + // Compute node sign distance function + for (const auto node_id : free_surface_candidate_nodes) { + const auto& node_coord = map_nodes_[node_id]->coordinates(); + double closest_distance = std::numeric_limits::max(); + double signed_distance = std::numeric_limits::max(); + for (const auto fs_id : free_surface_particles) { + const auto& fs_particle = map_particles_[fs_id]; + const VectorDim fs_coord = + fs_particle->coordinates() + + 0.5 * fs_particle->diameter() * fs_particle->normal(); + const VectorDim rel_coord = fs_coord - node_coord; + const double distance = rel_coord.norm(); + if (distance < closest_distance) { + closest_distance = distance; + signed_distance = (rel_coord).dot(fs_particle->normal()); + } + } + + // Assign signed distance to node + map_nodes_[node_id]->assign_signed_distance(signed_distance); + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return status; +} + +//! Compute free surface cells, nodes, and particles by density +template +bool mpm::Mesh::compute_free_surface_by_density(double volume_tolerance) { + bool status = true; + try { + // Reset free surface cell + this->iterate_over_cells(std::bind(&mpm::Cell::assign_free_surface, + std::placeholders::_1, false)); + + // Reset volume fraction + this->iterate_over_cells(std::bind(&mpm::Cell::assign_volume_fraction, + std::placeholders::_1, 0.0)); + + // Reset free surface particle + this->iterate_over_particles( + std::bind(&mpm::ParticleBase::assign_free_surface, + std::placeholders::_1, false)); + + // Compute and assign volume fraction to each cell + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + if ((*citr)->status()) { + // Compute volume fraction + double cell_volume_fraction = 0.0; + for (const auto p_id : (*citr)->particles()) + cell_volume_fraction += map_particles_[p_id]->volume(); + + cell_volume_fraction = cell_volume_fraction / (*citr)->volume(); + (*citr)->assign_volume_fraction(cell_volume_fraction); + } + } + + // Compute boundary cells and nodes based on geometry + std::set boundary_cells; + std::set boundary_nodes; + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); + ++citr) { + + if ((*citr)->status()) { + bool cell_at_interface = false; + const auto& node_id = (*citr)->nodes_id(); + bool internal = true; + + //! Check internal cell + for (const auto c_id : (*citr)->neighbours()) { + if (!map_cells_[c_id]->status()) { + internal = false; + break; + } + } + + //! Check volume fraction only for boundary cell + if (!internal) { + if ((*citr)->volume_fraction() < volume_tolerance) { + cell_at_interface = true; + for (const auto id : node_id) { + map_nodes_[id]->assign_free_surface(cell_at_interface); + boundary_nodes.insert(id); + } + } else { + for (const auto n_id : (*citr)->neighbours()) { + if (map_cells_[n_id]->volume_fraction() < volume_tolerance) { + cell_at_interface = true; + const auto& n_node_id = map_cells_[n_id]->nodes_id(); + + // Detect common node id + std::set common_node_id; + std::set_intersection( + node_id.begin(), node_id.end(), n_node_id.begin(), + n_node_id.end(), + std::inserter(common_node_id, common_node_id.begin())); + + // Assign free surface nodes + if (!common_node_id.empty()) { + for (const auto common_id : common_node_id) { + map_nodes_[common_id]->assign_free_surface( + cell_at_interface); + boundary_nodes.insert(common_id); + } + } + } + } + } + + // Assign free surface cell + if (cell_at_interface) { + (*citr)->assign_free_surface(cell_at_interface); + boundary_cells.insert((*citr)->id()); + } + } + } + } + + // Compute boundary particles based on density function + // Lump cell volume to nodes + this->iterate_over_cells(std::bind( + &mpm::Cell::map_cell_volume_to_nodes, std::placeholders::_1, 0)); + + // Compute nodal value of mass density + this->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::compute_density, std::placeholders::_1), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + // Evaluate free surface particles + std::set boundary_particles; + for (auto pitr = this->particles_.cbegin(); pitr != this->particles_.cend(); + ++pitr) { + bool status = (*pitr)->compute_free_surface_by_density(); + if (status) { + (*pitr)->assign_free_surface(status); + boundary_particles.insert((*pitr)->id()); + } + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return status; +} + +//! Get free surface node set +template +std::set mpm::Mesh::free_surface_nodes() { + std::set id_set; + for (auto nitr = this->nodes_.cbegin(); nitr != this->nodes_.cend(); ++nitr) + if ((*nitr)->free_surface()) id_set.insert((*nitr)->id()); + return id_set; +} + +//! Get free surface cell set +template +std::set mpm::Mesh::free_surface_cells() { + std::set id_set; + for (auto citr = this->cells_.cbegin(); citr != this->cells_.cend(); ++citr) + if ((*citr)->free_surface()) id_set.insert((*citr)->id()); + return id_set; +} + +//! Get free surface particle set +template +std::set mpm::Mesh::free_surface_particles() { + std::set id_set; + for (auto pitr = this->particles_.cbegin(); pitr != this->particles_.cend(); + ++pitr) + if ((*pitr)->free_surface()) id_set.insert((*pitr)->id()); + return id_set; +} + // Initialise the nodal properties' map template void mpm::Mesh::initialise_nodal_properties() { diff --git a/include/node.h b/include/node.h index 31a6b9a03..58a8083e3 100644 --- a/include/node.h +++ b/include/node.h @@ -245,6 +245,78 @@ class Node : public NodeBase { //! Set ghost id void ghost_id(Index gid) override { ghost_id_ = gid; } + //! Return real density at a given node for a given phase + //! \param[in] phase Index corresponding to the phase + double density(unsigned phase) override { return density_(phase); } + + //! Compute nodal density + void compute_density() override; + + //! Compute nodal correction force term + bool compute_nodal_correction_force( + const VectorDim& correction_force) override; + + //! Assign free surface + void assign_free_surface(bool free_surface) override { + free_surface_ = free_surface; + } + + //! Return free surface bool + bool free_surface() override { return free_surface_; } + + //! Assign signed distance + void assign_signed_distance(double signed_distance) override { + signed_distance_ = signed_distance; + } + + //! Return signed distance + double signed_distance() override { return signed_distance_; } + + //! Assign active id + void assign_active_id(Index id) override { active_id_ = id; } + + //! Return active id + mpm::Index active_id() override { return active_id_; } + + //! Return nodal pressure constraint + //! \param[in] phase Index corresponding to the phase + //! \param[in] current_time current time of the analysis + //! \retval pressure constraint at proper time for given phase + double pressure_constraint(const unsigned phase, + const double current_time) override { + if (pressure_constraints_.find(phase) != pressure_constraints_.end()) { + const double scalar = + (pressure_function_.find(phase) != pressure_function_.end()) + ? pressure_function_[phase]->value(current_time) + : 1.0; + + return scalar * pressure_constraints_[phase]; + } else + return std::numeric_limits::max(); + } + + //! Assign pressure constraint + //! \param[in] phase Index corresponding to the phase + //! \param[in] pressure Applied pressure constraint + //! \param[in] function math function + bool assign_pressure_constraint( + const unsigned phase, const double pressure, + const std::shared_ptr& function) override; + + //! Update pressure increment at the node + void update_pressure_increment(const Eigen::VectorXd& pressure_increment, + unsigned phase, + double current_time = 0.) override; + + //! Return nodal pressure increment + double pressure_increment() const override { return pressure_increment_; } + + //! Compute navier-stokes semi-implicit acceleration and velocity + //! \param[in] phase Index corresponding to the phase + //! \param[in] dt Timestep in analysis + //! \retval status Computation status + bool compute_acceleration_velocity_navierstokes_semi_implicit( + unsigned phase, double dt) override; //! Update nodal property at the nodes from particle //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] property Property name @@ -299,6 +371,8 @@ class Node : public NodeBase { Eigen::Matrix acceleration_; //! Velocity constraints std::map velocity_constraints_; + //! Pressure constraint + std::map pressure_constraints_; //! Rotation matrix for general velocity constraints Eigen::Matrix rotation_matrix_; //! Material ids whose information was passed to this node @@ -313,12 +387,26 @@ class Node : public NodeBase { Eigen::Matrix concentrated_force_; //! Mathematical function for force std::shared_ptr force_function_{nullptr}; + //! Mathematical function for pressure + std::map> pressure_function_; //! Nodal property pool std::shared_ptr property_handle_{nullptr}; //! Logger std::unique_ptr console_; //! MPI ranks std::set mpi_ranks_; + //! Global index for active node + Index active_id_{std::numeric_limits::max()}; + //! Interpolated density + Eigen::Matrix density_; + //! p^(t+1) - beta * p^(t) + double pressure_increment_; + //! Correction force + Eigen::Matrix correction_force_; + //! Free surface + bool free_surface_{false}; + //! Signed distance + double signed_distance_; }; // Node class } // namespace mpm diff --git a/include/node.tcc b/include/node.tcc index 7d3a9cfcb..478d4cca0 100644 --- a/include/node.tcc +++ b/include/node.tcc @@ -32,6 +32,9 @@ void mpm::Node::initialise() noexcept { velocity_.setZero(); momentum_.setZero(); acceleration_.setZero(); + free_surface_ = false; + pressure_increment_ = 0.; + correction_force_.setZero(); status_ = false; material_ids_.clear(); } @@ -537,6 +540,121 @@ bool mpm::Node::mpi_rank(unsigned rank) { return status.second; } +//! Compute mass density (Z. Wiezckowski, 2004) +//! density = mass / lumped volume +template +void mpm::Node::compute_density() { + const double tolerance = 1.E-16; // std::numeric_limits::lowest(); + + for (unsigned phase = 0; phase < Tnphases; ++phase) { + if (mass_(phase) > tolerance) { + if (volume_(phase) > tolerance) + density_(phase) = mass_(phase) / volume_(phase); + + // Check to see if value is below threshold + if (std::abs(density_(phase)) < 1.E-15) density_(phase) = 0.; + + } else + throw std::runtime_error("Nodal mass is zero or below threshold"); + } +} + +//! Assign pressure constraint +template +bool mpm::Node::assign_pressure_constraint( + const unsigned phase, const double pressure, + const std::shared_ptr& function) { + bool status = true; + try { + // Constrain directions can take values between 0 and Tnphases + if (phase < Tnphases) { + this->pressure_constraints_.insert(std::make_pair( + static_cast(phase), static_cast(pressure))); + // Assign pressure function + if (function != nullptr) + this->pressure_function_.insert( + std::make_pair>( + static_cast(phase), + static_cast>(function))); + } else + throw std::runtime_error("Pressure constraint phase is out of bounds"); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Update pressure increment at the node +template +void mpm::Node::update_pressure_increment( + const Eigen::VectorXd& pressure_increment, unsigned phase, + double current_time) { + this->pressure_increment_ = pressure_increment(active_id_); + + // If pressure boundary, increment is zero + if (pressure_constraints_.find(phase) != pressure_constraints_.end() || + this->free_surface()) + this->pressure_increment_ = 0; +} + +//! Compute nodal correction force +template +bool mpm::Node::compute_nodal_correction_force( + const VectorDim& correction_force) { + bool status = true; + + try { + // Compute correction force for water phase + correction_force_.col(0) = correction_force; + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Compute semi-implicit acceleration and velocity +template +bool mpm::Node:: + compute_acceleration_velocity_navierstokes_semi_implicit(unsigned phase, + double dt) { + bool status = true; + const double tolerance = std::numeric_limits::min(); + try { + + Eigen::Matrix acceleration_corrected = + correction_force_.col(phase) / mass_(phase); + + // Acceleration + this->acceleration_.col(phase) += acceleration_corrected; + + // Update velocity + velocity_.col(phase) += acceleration_corrected * dt; + + // Apply friction constraints + this->apply_friction_constraints(dt); + + // Apply velocity constraints, which also sets acceleration to 0, + // when velocity is set. + this->apply_velocity_constraints(); + + // Set a threshold + for (unsigned i = 0; i < Tdim; ++i) { + if (std::abs(velocity_.col(phase)(i)) < tolerance) + velocity_.col(phase)(i) = 0.; + if (std::abs(acceleration_.col(phase)(i)) < tolerance) + acceleration_.col(phase)(i) = 0.; + } + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + //! Update nodal property at the nodes from particle template void mpm::Node::update_property( diff --git a/include/node_base.h b/include/node_base.h index 17eddfdca..c52caf833 100644 --- a/include/node_base.h +++ b/include/node_base.h @@ -234,6 +234,61 @@ class NodeBase { //! Set ghost id virtual void ghost_id(Index gid) = 0; + //! Return real density at a given node for a given phase + //! \param[in] phase Index corresponding to the phase + virtual double density(unsigned phase) = 0; + + //! Compute nodal density + virtual void compute_density() = 0; + + //! Assign free surface + virtual void assign_free_surface(bool free_surface) = 0; + + //! Return free surface bool + virtual bool free_surface() = 0; + + //! Assign signed distance + virtual void assign_signed_distance(double signed_distance) = 0; + + //! Return signed distance + virtual double signed_distance() = 0; + + //! Assign active id + virtual void assign_active_id(Index id) = 0; + + //! Return active id + virtual mpm::Index active_id() = 0; + + //! Compute navier-stokes semi-implicit acceleration and velocity + //! \param[in] phase Index corresponding to the phase + //! \param[in] dt Timestep in analysis + //! \retval status Computation status + virtual bool compute_acceleration_velocity_navierstokes_semi_implicit( + unsigned phase, double dt) = 0; + + //! Return pressure constraint + virtual double pressure_constraint(const unsigned phase, + const double current_time) = 0; + + //! Assign pressure constraint + //! \param[in] phase Index corresponding to the phase + //! \param[in] pressure Applied pressure constraint + //! \param[in] function math function + virtual bool assign_pressure_constraint( + const unsigned phase, double pressure, + const std::shared_ptr& function) = 0; + + //! Update pressure increment at the node + virtual void update_pressure_increment( + const Eigen::VectorXd& pressure_increment, unsigned phase, + double current_time = 0.) = 0; + + //! Return nodal pressure increment + virtual double pressure_increment() const = 0; + + //! Compute nodal correction force term + virtual bool compute_nodal_correction_force( + const VectorDim& correction_force) = 0; //! Update nodal property at the nodes from particle //! \param[in] update A boolean to update (true) or assign (false) //! \param[in] property Property name diff --git a/include/particles/fluid_particle.h b/include/particles/fluid_particle.h new file mode 100644 index 000000000..17fdbd9ae --- /dev/null +++ b/include/particles/fluid_particle.h @@ -0,0 +1,104 @@ +#ifndef MPM_FLUID_PARTICLE_H_ +#define MPM_FLUID_PARTICLE_H_ + +#include +#include +#include +#include +#include + +#include "logger.h" +#include "particle.h" + +namespace mpm { + +//! Fluid Particle class +//! \brief Class with function specific to fluid particles +//! \tparam Tdim Dimension +template +class FluidParticle : public mpm::Particle { + public: + //! Define a vector of size dimension + using VectorDim = Eigen::Matrix; + + //! Construct a particle with id and coordinates + //! \param[in] id Particle id + //! \param[in] coord Coordinates of the particles + FluidParticle(Index id, const VectorDim& coord); + + //! Destructor + ~FluidParticle() override{}; + + //! Delete copy constructor + FluidParticle(const FluidParticle&) = delete; + + //! Delete assignment operator + FluidParticle& operator=(const FluidParticle&) = delete; + + //! Compute stress + void compute_stress() noexcept override; + + //! Map internal force + inline void map_internal_force() noexcept override; + + //! ---------------------------------------------------------------- + //! Semi-Implicit integration functions based on Chorin's Projection + //! ---------------------------------------------------------------- + + //! Assigning beta parameter to particle + //! \param[in] pressure parameter determining type of projection + void assign_projection_parameter(double parameter) override { + this->projection_param_ = parameter; + }; + + //! Map laplacian element matrix to cell (used in poisson equation LHS) + bool map_laplacian_to_cell() override; + + //! Map poisson rhs element matrix to cell (used in poisson equation RHS) + bool map_poisson_right_to_cell() override; + + //! Map correction matrix element matrix to cell (used to correct velocity) + bool map_correction_matrix_to_cell() override; + + //! Update pressure after solving poisson equation + bool compute_updated_pressure() override; + + private: + //! Compute turbulent stress + virtual Eigen::Matrix compute_turbulent_stress(); + + private: + //! Cell + using ParticleBase::cell_; + //! Nodes + using ParticleBase::nodes_; + //! Fluid material + using ParticleBase::material_; + //! State variables + using ParticleBase::state_variables_; + //! Shape functions + using Particle::shapefn_; + //! dN/dX + using Particle::dn_dx_; + //! Fluid strain rate + using Particle::strain_rate_; + //! Effective stress of soil skeleton + using Particle::stress_; + //! Particle mass density + using Particle::mass_density_; + //! Particle mass density + using Particle::mass_; + //! Particle total volume + using Particle::volume_; + //! Projection parameter for semi-implicit update + double projection_param_{0.0}; + //! Pressure constraint + double pressure_constraint_{std::numeric_limits::max()}; + //! Logger + std::unique_ptr console_; +}; // FluidParticle class +} // namespace mpm + +#include "fluid_particle.tcc" + +#endif // MPM_FLUID_PARTICLE_H__ diff --git a/include/particles/fluid_particle.tcc b/include/particles/fluid_particle.tcc new file mode 100644 index 000000000..65050f5ed --- /dev/null +++ b/include/particles/fluid_particle.tcc @@ -0,0 +1,206 @@ +//! Construct a two phase particle with id and coordinates +template +mpm::FluidParticle::FluidParticle(Index id, const VectorDim& coord) + : mpm::Particle(id, coord) { + + // Logger + std::string logger = + "FluidParticle" + std::to_string(Tdim) + "d::" + std::to_string(id); + console_ = std::make_unique(logger, mpm::stdout_sink); +} + +// Compute stress +template +void mpm::FluidParticle::compute_stress() noexcept { + // Run particle compute stress + mpm::Particle::compute_stress(); + + // Calculate fluid turbulent stress + this->stress_ += this->compute_turbulent_stress(); +} + +// Compute turbulent stress +template +Eigen::Matrix + mpm::FluidParticle::compute_turbulent_stress() { + // Compute turbulent stress depends on the model + Eigen::Matrix tstress; + tstress.setZero(); + + // Apply LES Smagorinsky closure + const double smagorinsky_constant = 0.2; + const double grid_spacing = std::pow(cell_->volume(), 1 / (double)Tdim); + const auto strain_rate = this->strain_rate(); + double local_strain_rate = 0.0; + local_strain_rate += + 2 * (strain_rate[0] * strain_rate[0] + strain_rate[1] * strain_rate[1] + + strain_rate[2] * strain_rate[2]) + + strain_rate[3] * strain_rate[3] + strain_rate[4] * strain_rate[4] + + strain_rate[5] * strain_rate[5]; + local_strain_rate = std::sqrt(local_strain_rate); + + // Turbulent Eddy Viscosity + const double turbulent_viscosity = + std::pow(smagorinsky_constant * grid_spacing, 2) * local_strain_rate; + + // Turbulent stress + tstress(0) = 2. * turbulent_viscosity / mass_density_ * strain_rate(0); + tstress(1) = 2. * turbulent_viscosity / mass_density_ * strain_rate(1); + tstress(2) = 2. * turbulent_viscosity / mass_density_ * strain_rate(2); + tstress(3) = turbulent_viscosity / mass_density_ * strain_rate(3); + tstress(4) = turbulent_viscosity / mass_density_ * strain_rate(4); + tstress(5) = turbulent_viscosity / mass_density_ * strain_rate(5); + + return tstress; +} + +//! Map internal force +template <> +inline void mpm::FluidParticle<1>::map_internal_force() noexcept { + // initialise a vector of total stress (deviatoric + turbulent - pressure) + Eigen::Matrix total_stress = this->stress_; + total_stress(0) -= + this->projection_param_ * + this->state_variables(mpm::ParticlePhase::SinglePhase)["pressure"]; + // Compute nodal internal forces + for (unsigned i = 0; i < nodes_.size(); ++i) { + // Compute force: -pstress * volume + Eigen::Matrix force; + force[0] = dn_dx_(i, 0) * total_stress[0]; + force *= -1 * this->volume_; + + nodes_[i]->update_internal_force(true, mpm::ParticlePhase::SinglePhase, + force); + } +} + +//! Map internal force +template <> +inline void mpm::FluidParticle<2>::map_internal_force() noexcept { + // initialise a vector of total stress (deviatoric + turbulent - pressure) + Eigen::Matrix total_stress = this->stress_; + total_stress(0) -= + this->projection_param_ * + this->state_variables(mpm::ParticlePhase::SinglePhase)["pressure"]; + total_stress(1) -= + this->projection_param_ * + this->state_variables(mpm::ParticlePhase::SinglePhase)["pressure"]; + + // Compute nodal internal forces + for (unsigned i = 0; i < nodes_.size(); ++i) { + // Compute force: -pstress * volume + Eigen::Matrix force; + force[0] = dn_dx_(i, 0) * total_stress[0] + dn_dx_(i, 1) * total_stress[3]; + force[1] = dn_dx_(i, 1) * total_stress[1] + dn_dx_(i, 0) * total_stress[3]; + + force *= -1. * this->volume_; + + nodes_[i]->update_internal_force(true, mpm::ParticlePhase::SinglePhase, + force); + } +} + +//! Map internal force +template <> +inline void mpm::FluidParticle<3>::map_internal_force() noexcept { + // initialise a vector of total stress (deviatoric + turbulent - pressure) + Eigen::Matrix total_stress = this->stress_; + total_stress(0) -= + this->projection_param_ * + this->state_variables(mpm::ParticlePhase::SinglePhase)["pressure"]; + total_stress(1) -= + this->projection_param_ * + this->state_variables(mpm::ParticlePhase::SinglePhase)["pressure"]; + total_stress(2) -= + this->projection_param_ * + this->state_variables(mpm::ParticlePhase::SinglePhase)["pressure"]; + + // Compute nodal internal forces + for (unsigned i = 0; i < nodes_.size(); ++i) { + // Compute force: -pstress * volume + Eigen::Matrix force; + force[0] = dn_dx_(i, 0) * total_stress[0] + dn_dx_(i, 1) * total_stress[3] + + dn_dx_(i, 2) * total_stress[5]; + + force[1] = dn_dx_(i, 1) * total_stress[1] + dn_dx_(i, 0) * total_stress[3] + + dn_dx_(i, 2) * total_stress[4]; + + force[2] = dn_dx_(i, 2) * total_stress[2] + dn_dx_(i, 1) * total_stress[4] + + dn_dx_(i, 0) * total_stress[5]; + + force *= -1. * this->volume_; + + nodes_[i]->update_internal_force(true, mpm::ParticlePhase::SinglePhase, + force); + } +} + +//! Map laplacian element matrix to cell (used in poisson equation LHS) +template +bool mpm::FluidParticle::map_laplacian_to_cell() { + bool status = true; + try { + // Compute local matrix of Laplacian + cell_->compute_local_laplacian(dn_dx_, volume_); + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Map poisson rhs element matrix to cell (used in poisson equation RHS) +template +bool mpm::FluidParticle::map_poisson_right_to_cell() { + bool status = true; + try { + // Compute local poisson rhs matrix + cell_->compute_local_poisson_right( + shapefn_, dn_dx_, volume_, + this->material(mpm::ParticlePhase::SinglePhase) + ->template property(std::string("density"))); + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +// Compute updated pressure of the particle based on nodal pressure +template +bool mpm::FluidParticle::compute_updated_pressure() { + bool status = true; + try { + double pressure_increment = 0; + for (unsigned i = 0; i < nodes_.size(); ++i) { + pressure_increment += shapefn_(i) * nodes_[i]->pressure_increment(); + } + + // Get interpolated nodal pressure + state_variables_[mpm::ParticlePhase::SinglePhase].at("pressure") = + state_variables_[mpm::ParticlePhase::SinglePhase].at("pressure") * + projection_param_ + + pressure_increment; + + // Overwrite pressure if free surface + if (this->free_surface()) + state_variables_[mpm::ParticlePhase::SinglePhase].at("pressure") = 0.0; + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Map correction matrix element matrix to cell (used to correct velocity) +template +bool mpm::FluidParticle::map_correction_matrix_to_cell() { + bool status = true; + try { + cell_->compute_local_correction_matrix(shapefn_, dn_dx_, volume_); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + } + return status; +} \ No newline at end of file diff --git a/include/particles/particle.h b/include/particles/particle.h index 849b72f56..e9f7843bb 100644 --- a/include/particles/particle.h +++ b/include/particles/particle.h @@ -59,16 +59,6 @@ class Particle : public ParticleBase { const HDF5Particle& particle, const std::shared_ptr>& material) override; - //! Assign material history variables - //! \param[in] state_vars State variables - //! \param[in] material Material associated with the particle - //! \param[in] phase Index to indicate material phase - //! \retval status Status of cloning HDF5 particle - bool assign_material_state_vars( - const mpm::dense_map& state_vars, - const std::shared_ptr>& material, - unsigned phase = mpm::ParticlePhase::Solid) override; - //! Retrun particle data as HDF5 //! \retval particle HDF5 data of the particle HDF5Particle hdf5() const override; @@ -121,6 +111,14 @@ class Particle : public ParticleBase { //! Return volume double volume() const override { return volume_; } + //! Return the approximate particle diameter + double diameter() const override { + double diameter = 0.; + if (Tdim == 2) diameter = 2.0 * std::sqrt(volume_ / M_PI); + if (Tdim == 3) diameter = 2.0 * std::pow(volume_ * 0.75 / M_PI, (1 / 3)); + return diameter; + } + //! Return size of particle in natural coordinates VectorDim natural_size() const override { return natural_size_; } @@ -131,7 +129,6 @@ class Particle : public ParticleBase { void update_volume() noexcept override; //! Return mass density - //! \param[in] phase Index corresponding to the phase double mass_density() const override { return mass_density_; } //! Compute mass as volume * density @@ -222,7 +219,6 @@ class Particle : public ParticleBase { bool assign_traction(unsigned direction, double traction) override; //! Return traction of the particle - //! \param[in] phase Index corresponding to the phase VectorDim traction() const override { return traction_; } //! Map traction force @@ -234,6 +230,27 @@ class Particle : public ParticleBase { void compute_updated_position(double dt, bool velocity_update = false) noexcept override; + //! Assign material history variables + //! \param[in] state_vars State variables + //! \param[in] material Material associated with the particle + //! \param[in] phase Index to indicate material phase + //! \retval status Status of cloning HDF5 particle + bool assign_material_state_vars( + const mpm::dense_map& state_vars, + const std::shared_ptr>& material, + unsigned phase = mpm::ParticlePhase::Solid) override; + + //! Assign a state variable + //! \param[in] var State variable + //! \param[in] value State variable to be assigned + //! \param[in] phase Index to indicate phase + void assign_state_variable( + const std::string& var, double value, + unsigned phase = mpm::ParticlePhase::Solid) override { + if (state_variables_[phase].find(var) != state_variables_[phase].end()) + state_variables_[phase].at(var) = value; + } + //! Return a state variable //! \param[in] var State variable //! \param[in] phase Index to indicate phase @@ -255,13 +272,18 @@ class Particle : public ParticleBase { bool compute_pressure_smoothing( unsigned phase = mpm::ParticlePhase::Solid) noexcept override; + //! Assign a state variable + //! \param[in] value Particle pressure to be assigned + //! \param[in] phase Index to indicate phase + void assign_pressure(double pressure, + unsigned phase = mpm::ParticlePhase::Solid) override { + this->assign_state_variable("pressure", pressure, phase); + } + //! Return pressure of the particles //! \param[in] phase Index to indicate phase double pressure(unsigned phase = mpm::ParticlePhase::Solid) const override { - return (state_variables_[phase].find("pressure") != - state_variables_[phase].end()) - ? state_variables_[phase].at("pressure") - : std::numeric_limits::quiet_NaN(); + return this->state_variable("pressure", phase); } //! Return tensor data of particles @@ -278,6 +300,26 @@ class Particle : public ParticleBase { //! Assign material id of this particle to nodes void append_material_id_to_nodes() const override; + //! Assign free surface + void assign_free_surface(bool free_surface) override { + free_surface_ = free_surface; + }; + + //! Return free surface bool + bool free_surface() override { return free_surface_; }; + + //! Compute free surface in particle level by density ratio comparison + //! \param[in] density_ratio_tolerance Tolerance of density ratio comparison + //! \retval status Status of compute_free_surface + bool compute_free_surface_by_density( + double density_ratio_tolerance = 0.70) override; + + //! Assign normal vector + void assign_normal(const VectorDim& normal) override { normal_ = normal; }; + + //! Return normal vector + VectorDim normal() override { return normal_; }; + //! Return the number of neighbour particles unsigned nneighbours() const override { return neighbours_.size(); }; @@ -306,7 +348,7 @@ class Particle : public ParticleBase { inline Eigen::Matrix compute_strain_rate( const Eigen::MatrixXd& dn_dx, unsigned phase) noexcept; - private: + protected: //! particle id using ParticleBase::id_; //! coordinates @@ -371,6 +413,10 @@ class Particle : public ParticleBase { std::unique_ptr console_; //! Map of vector properties std::map> properties_; + //! Free surface + bool free_surface_{false}; + //! Free surface + Eigen::Matrix normal_; }; // Particle class } // namespace mpm diff --git a/include/particles/particle.tcc b/include/particles/particle.tcc index c3d945b69..31b316eea 100644 --- a/include/particles/particle.tcc +++ b/include/particles/particle.tcc @@ -240,6 +240,7 @@ void mpm::Particle::initialise() { stress_.setZero(); traction_.setZero(); velocity_.setZero(); + normal_.setZero(); volume_ = std::numeric_limits::max(); volumetric_strain_centroid_ = 0.; @@ -248,6 +249,7 @@ void mpm::Particle::initialise() { this->properties_["strains"] = [&]() { return strain(); }; this->properties_["velocities"] = [&]() { return velocity(); }; this->properties_["displacements"] = [&]() { return displacement(); }; + this->properties_["normals"] = [&]() { return normal(); }; } //! Initialise particle material container @@ -839,6 +841,10 @@ bool mpm::Particle::compute_pressure_smoothing(unsigned phase) noexcept { pressure += shapefn_[i] * nodes_[i]->pressure(phase); state_variables_[phase]["pressure"] = pressure; + + // If free_surface particle, overwrite pressure to zero + if (free_surface_) state_variables_[phase]["pressure"] = 0.0; + status = true; } return status; @@ -865,6 +871,27 @@ void mpm::Particle::append_material_id_to_nodes() const { nodes_[i]->append_material_id(this->material_id()); } +//! Compute free surface in particle level by density ratio comparison +template +bool mpm::Particle::compute_free_surface_by_density( + double density_ratio_tolerance) { + bool status = false; + // Check if particle has a valid cell ptr + if (cell_ != nullptr) { + // Simple approach of density comparison (Hamad, 2015) + // Get interpolated nodal density + double nodal_mass_density = 0; + for (unsigned i = 0; i < nodes_.size(); ++i) + nodal_mass_density += + shapefn_[i] * nodes_[i]->density(mpm::ParticlePhase::Solid); + + // Compare smoothen density to actual particle mass density + if ((nodal_mass_density / mass_density_) <= density_ratio_tolerance) + status = true; + } + return status; +}; + //! Assign neighbour particles template void mpm::Particle::assign_neighbours( diff --git a/include/particles/particle_base.h b/include/particles/particle_base.h index da0d22a99..cfbffb953 100644 --- a/include/particles/particle_base.h +++ b/include/particles/particle_base.h @@ -19,7 +19,12 @@ template class Material; //! Particle phases -enum ParticlePhase : unsigned int { Solid = 0, Liquid = 1, Gas = 2 }; +enum ParticlePhase : unsigned int { + SinglePhase = 0, + Solid = 0, + Liquid = 1, + Gas = 2 +}; //! ParticleBase class //! \brief Base class that stores the information about particleBases @@ -64,16 +69,6 @@ class ParticleBase { const HDF5Particle& particle, const std::shared_ptr>& material) = 0; - //! Assign material history variables - //! \param[in] state_vars State variables - //! \param[in] material Material associated with the particle - //! \param[in] phase Index to indicate material phase - //! \retval status Status of cloning HDF5 particle - virtual bool assign_material_state_vars( - const mpm::dense_map& state_vars, - const std::shared_ptr>& material, - unsigned phase = mpm::ParticlePhase::Solid) = 0; - //! Retrun particle data as HDF5 //! \retval particle HDF5 data of the particle virtual HDF5Particle hdf5() const = 0; @@ -123,6 +118,9 @@ class ParticleBase { //! Return volume virtual double volume() const = 0; + //! Return the approximate particle diameter + virtual double diameter() const = 0; + //! Return size of particle in natural coordinates virtual VectorDim natural_size() const = 0; @@ -167,6 +165,12 @@ class ParticleBase { return material_id_[phase]; } + //! Assign material state variables + virtual bool assign_material_state_vars( + const mpm::dense_map& state_vars, + const std::shared_ptr>& material, + unsigned phase = mpm::ParticlePhase::Solid) = 0; + //! Return state variables //! \param[in] phase Index to indicate material phase mpm::dense_map state_variables( @@ -174,6 +178,16 @@ class ParticleBase { return state_variables_[phase]; } + //! Assign a state variable + virtual void assign_state_variable( + const std::string& var, double value, + unsigned phase = mpm::ParticlePhase::Solid) = 0; + + //! Return a state variable + virtual double state_variable( + const std::string& var, + unsigned phase = mpm::ParticlePhase::Solid) const = 0; + //! Assign status void assign_status(bool status) { status_ = status; } @@ -189,6 +203,10 @@ class ParticleBase { //! Return mass virtual double mass() const = 0; + //! Assign pressure + virtual void assign_pressure(double pressure, + unsigned phase = mpm::ParticlePhase::Solid) = 0; + //! Return pressure virtual double pressure(unsigned phase = mpm::ParticlePhase::Solid) const = 0; @@ -208,7 +226,7 @@ class ParticleBase { virtual double dvolumetric_strain() const = 0; //! Initial stress - virtual void initial_stress(const Eigen::Matrix&) = 0; + virtual void initial_stress(const Eigen::Matrix& stress) = 0; //! Compute stress virtual void compute_stress() noexcept = 0; @@ -252,11 +270,6 @@ class ParticleBase { virtual void compute_updated_position( double dt, bool velocity_update = false) noexcept = 0; - //! Return a state variable - virtual double state_variable( - const std::string& var, - unsigned phase = mpm::ParticlePhase::Solid) const = 0; - //! Return tensor data of particles //! \param[in] property Property string //! \retval vecdata Tensor data of particle property @@ -271,6 +284,22 @@ class ParticleBase { //! Assign material id of this particle to nodes virtual void append_material_id_to_nodes() const = 0; + //! Assign particle free surface + virtual void assign_free_surface(bool free_surface) = 0; + + //! Assign particle free surface + virtual bool free_surface() = 0; + + //! Compute free surface in particle level by density ratio comparison + virtual bool compute_free_surface_by_density( + double density_ratio_tolerance = 0.70) = 0; + + //! Assign normal vector + virtual void assign_normal(const VectorDim& normal) = 0; + + //! Return normal vector + virtual VectorDim normal() = 0; + //! Return the number of neighbour particles virtual unsigned nneighbours() const = 0; @@ -282,6 +311,49 @@ class ParticleBase { //! Return neighbour ids virtual std::vector neighbours() const = 0; + //! Assigning beta parameter to particle + //! \param[in] pressure parameter determining type of projection + virtual void assign_projection_parameter(double parameter) { + throw std::runtime_error( + "Calling the base class function (assign_projection_parameter) in " + "ParticleBase:: illegal operation!"); + }; + + //! Map laplacian element matrix to cell (used in poisson equation LHS) + virtual bool map_laplacian_to_cell() { + throw std::runtime_error( + "Calling the base class function (map_laplacian_to_cell) in " + "ParticleBase:: " + "illegal operation!"); + return 0; + }; + + //! Map poisson rhs element matrix to cell (used in poisson equation RHS) + virtual bool map_poisson_right_to_cell() { + throw std::runtime_error( + "Calling the base class function (map_poisson_right_to_cell) in " + "ParticleBase:: " + "illegal operation!"); + return 0; + }; + + //! Map correction matrix element matrix to cell (used to correct velocity) + virtual bool map_correction_matrix_to_cell() { + throw std::runtime_error( + "Calling the base class function (map_correction_matrix_to_cell) in " + "ParticleBase:: " + "illegal operation!"); + return 0; + }; + + //! Update pressure after solving poisson equation + virtual bool compute_updated_pressure() { + throw std::runtime_error( + "Calling the base class function (compute_updated_pressure) in " + "ParticleBase:: illegal operation!"); + return 0; + }; + protected: //! particleBase id Index id_{std::numeric_limits::max()}; diff --git a/include/solvers/mpm_base.h b/include/solvers/mpm_base.h index ca017f1fc..bdb327569 100644 --- a/include/solvers/mpm_base.h +++ b/include/solvers/mpm_base.h @@ -80,14 +80,18 @@ class MPMBase : public MPM { //! Write HDF5 files void write_hdf5(mpm::Index step, mpm::Index max_steps) override; - //! Domain decomposition - //! \param[in] initial_step Start of simulation or later steps - void mpi_domain_decompose(bool initial_step = false) override; - //! Pressure smoothing //! \param[in] phase Phase to smooth pressure void pressure_smoothing(unsigned phase); + protected: + //! Locate particle + void locate_particle(); + + //! Domain decomposition + //! \param[in] initial_step Start of simulation or later steps + void mpi_domain_decompose(bool initial_step = false) override; + private: //! Return if a mesh will be isoparametric or not //! \retval isoparametric Status of mesh type @@ -147,6 +151,12 @@ class MPMBase : public MPM { const Json& mesh_prop, const std::shared_ptr>& particle_io); + //! Nodal pressure constraints + //! \param[in] mesh_prop Mesh properties + //! \param[in] mesh_io Mesh IO handle + void nodal_pressure_constraints( + const Json& mesh_prop, const std::shared_ptr>& mesh_io); + //! Particle entity sets //! \param[in] mesh_prop Mesh properties //! \param[in] check Check duplicates diff --git a/include/solvers/mpm_base.tcc b/include/solvers/mpm_base.tcc index 8426ccf4f..4ddbb5b9e 100644 --- a/include/solvers/mpm_base.tcc +++ b/include/solvers/mpm_base.tcc @@ -184,6 +184,9 @@ bool mpm::MPMBase::initialise_mesh() { // Read and assign friction constraints this->nodal_frictional_constraints(mesh_props, mesh_io); + // Read and assign friction constraints + this->nodal_pressure_constraints(mesh_props, mesh_io); + // Initialise cell auto cells_begin = std::chrono::steady_clock::now(); // Shape function name @@ -904,6 +907,72 @@ void mpm::MPMBase::nodal_frictional_constraints( } } +// Nodal pressure constraints +template +void mpm::MPMBase::nodal_pressure_constraints( + const Json& mesh_props, const std::shared_ptr>& mesh_io) { + try { + + unsigned phase = 0; + + // Read and assign pressure constraints + if (mesh_props.find("boundary_conditions") != mesh_props.end() && + mesh_props["boundary_conditions"].find("pressure_constraints") != + mesh_props["boundary_conditions"].end()) { + + // Iterate over pressure constraints + for (const auto& constraints : + mesh_props["boundary_conditions"]["pressure_constraints"]) { + + // Check if it is pressure increment constraints + if (constraints.find("increment_boundary") != constraints.end() && + constraints["increment_boundary"]) + phase += 1; + + // Pressure constraints are specified in a file + if (constraints.find("file") != constraints.end()) { + + std::string pressure_constraints_file = + constraints.at("file").template get(); + bool ppressure_constraints = + constraints_->assign_nodal_pressure_constraints( + phase, mesh_io->read_pressure_constraints( + io_->file_name(pressure_constraints_file))); + if (!ppressure_constraints) + throw std::runtime_error( + "Pressure constraints are not properly assigned"); + } else { + + // Get the math function + std::shared_ptr pfunction = nullptr; + if (constraints.find("math_function_id") != constraints.end()) + pfunction = math_functions_.at( + constraints.at("math_function_id").template get()); + + // Set id + int nset_id = constraints.at("nset_id").template get(); + // Pressure + double pressure = constraints.at("pressure").template get(); + // Phase if available + if (constraints.contains("phase")) + phase = constraints.at("phase").template get(); + // Add friction constraint to mesh + auto pressure_constraint = std::make_shared( + nset_id, phase, pressure); + // Add pressure constraint to mesh + constraints_->assign_nodal_pressure_constraint(pfunction, + pressure_constraint); + } + } + } else + throw std::runtime_error("Pressure constraints JSON not found"); + + } catch (std::exception& exception) { + console_->warn("#{}: Pressure conditions are undefined {} ", __LINE__, + exception.what()); + } +} + //! Cell entity sets template void mpm::MPMBase::cell_entity_sets(const Json& mesh_props, @@ -1096,6 +1165,20 @@ bool mpm::MPMBase::initialise_damping(const Json& damping_props) { return status; } +//! Locate particle +template +void mpm::MPMBase::locate_particle() { + // Locate particle + auto unlocatable_particles = mesh_->locate_particles_mesh(); + + if (!unlocatable_particles.empty() && this->locate_particles_) + throw std::runtime_error("Particle outside the mesh domain"); + // If unable to locate particles remove particles + if (!unlocatable_particles.empty() && !this->locate_particles_) + for (const auto& remove_particle : unlocatable_particles) + mesh_->remove_particle(remove_particle); +} + //! Domain decomposition template void mpm::MPMBase::mpi_domain_decompose(bool initial_step) { diff --git a/include/solvers/mpm_explicit.tcc b/include/solvers/mpm_explicit.tcc index 67261d041..955f21b17 100644 --- a/include/solvers/mpm_explicit.tcc +++ b/include/solvers/mpm_explicit.tcc @@ -292,15 +292,8 @@ bool mpm::MPMExplicit::solve() { if (this->stress_update_ == mpm::StressUpdate::USL) this->compute_stress_strain(phase); - // Locate particles - auto unlocatable_particles = mesh_->locate_particles_mesh(); - - if (!unlocatable_particles.empty() && this->locate_particles_) - throw std::runtime_error("Particle outside the mesh domain"); - // If unable to locate particles remove particles - if (!unlocatable_particles.empty() && !this->locate_particles_) - for (const auto& remove_particle : unlocatable_particles) - mesh_->remove_particle(remove_particle); + // Locate particle + this->locate_particle(); #ifdef USE_MPI #ifdef USE_GRAPH_PARTITIONING diff --git a/include/solvers/mpm_semi_implicit_navierstokes.h b/include/solvers/mpm_semi_implicit_navierstokes.h new file mode 100644 index 000000000..5e0fbe6af --- /dev/null +++ b/include/solvers/mpm_semi_implicit_navierstokes.h @@ -0,0 +1,95 @@ +#ifndef MPM_MPM_SEMI_IMPLICIT_NAVIER_STOKES_H_ +#define MPM_MPM_SEMI_IMPLICIT_NAVIER_STOKES_H_ + +#ifdef USE_GRAPH_PARTITIONING +#include "graph.h" +#endif + +#include "mpm_base.h" + +#include "assembler_base.h" +#include "cg_eigen.h" +#include "solver_base.h" + +namespace mpm { + +//! MPMSemiImplicit Navier Stokes class +//! \brief A class that implements the fractional step navier-stokes mpm +//! \tparam Tdim Dimension +template +class MPMSemiImplicitNavierStokes : public MPMBase { + public: + //! Default constructor + MPMSemiImplicitNavierStokes(const std::shared_ptr& io); + + //! Return matrix assembler pointer + std::shared_ptr> matrix_assembler() { + return assembler_; + } + + //! Solve + bool solve() override; + + //! Class private functions + private: + //! Initialise matrix + virtual bool initialise_matrix(); + + //! Initialise matrix + virtual bool reinitialise_matrix(); + + //! Compute poisson equation + bool compute_poisson_equation(std::string solver_type = "cg"); + + //! Compute corrected velocity + bool compute_correction_force(); + + //! Class private variables + private: + // Generate a unique id for the analysis + using mpm::MPMBase::uuid_; + //! Time step size + using mpm::MPMBase::dt_; + //! Current step + using mpm::MPMBase::step_; + //! Number of steps + using mpm::MPMBase::nsteps_; + //! Output steps + using mpm::MPMBase::output_steps_; + //! A unique ptr to IO object + using mpm::MPMBase::io_; + //! JSON analysis object + using mpm::MPMBase::analysis_; + //! JSON post-process object + using mpm::MPMBase::post_process_; + //! Logger + using mpm::MPMBase::console_; + //! Stress update + using mpm::MPMBase::stress_update_; + //! velocity update + using mpm::MPMBase::velocity_update_; + //! Gravity + using mpm::MPMBase::gravity_; + //! Mesh object + using mpm::MPMBase::mesh_; + //! Materials + using mpm::MPMBase::materials_; + //! Pressure smoothing + bool pressure_smoothing_{false}; + // Projection method paramter (beta) + double beta_{1}; + //! Assembler object + std::shared_ptr> assembler_; + //! Linear solver object + std::shared_ptr>> linear_solver_; + //! Method to detect free surface detection + std::string free_surface_detection_; + //! Volume tolerance for free surface + double volume_tolerance_{0}; + +}; // MPMSemiImplicit class +} // namespace mpm + +#include "mpm_semi_implicit_navierstokes.tcc" + +#endif // MPM_MPM_SEMI_IMPLICIT_NAVIER_STOKES_H_ diff --git a/include/solvers/mpm_semi_implicit_navierstokes.tcc b/include/solvers/mpm_semi_implicit_navierstokes.tcc new file mode 100644 index 000000000..1383489f7 --- /dev/null +++ b/include/solvers/mpm_semi_implicit_navierstokes.tcc @@ -0,0 +1,383 @@ +//! Constructor +template +mpm::MPMSemiImplicitNavierStokes::MPMSemiImplicitNavierStokes( + const std::shared_ptr& io) + : mpm::MPMBase(io) { + //! Logger + console_ = spdlog::get("MPMSemiImplicitNavierStokes"); +} + +//! MPM semi-implicit two phase solver +template +bool mpm::MPMSemiImplicitNavierStokes::solve() { + bool status = true; + + console_->info("MPM analysis type {}", io_->analysis_type()); + + // Initialise MPI rank and size + int mpi_rank = 0; + int mpi_size = 1; + +#ifdef USE_MPI + // Get MPI rank + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); + // Get number of MPI ranks + MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); +#endif + + // This solver consider only fluid variables + // NOTE: Due to indexing purposes + const unsigned fluid = mpm::ParticlePhase::SinglePhase; + + // Test if checkpoint resume is needed + bool resume = false; + if (analysis_.find("resume") != analysis_.end()) + resume = analysis_["resume"]["resume"].template get(); + + // Pressure smoothing + if (analysis_.find("pressure_smoothing") != analysis_.end()) + pressure_smoothing_ = analysis_["pressure_smoothing"].template get(); + + // Projection method parameter (beta) + if (analysis_.find("semi_implicit") != analysis_.end()) + beta_ = analysis_["semi_implicit"]["beta"].template get(); + + // Initialise material + bool mat_status = this->initialise_materials(); + if (!mat_status) { + status = false; + throw std::runtime_error("Initialisation of materials failed"); + } + + // Initialise mesh + bool mesh_status = this->initialise_mesh(); + if (!mesh_status) { + status = false; + throw std::runtime_error("Initialisation of mesh failed"); + } + + // Initialise particles + bool particle_status = this->initialise_particles(); + if (!particle_status) { + status = false; + throw std::runtime_error("Initialisation of particles failed"); + } + + // Initialise loading conditions + bool loading_status = this->initialise_loads(); + if (!loading_status) { + status = false; + throw std::runtime_error("Initialisation of loads failed"); + } + + // Initialise matrix + bool matrix_status = this->initialise_matrix(); + if (!matrix_status) { + status = false; + throw std::runtime_error("Initialisation of matrix failed"); + } + + // Compute mass for each phase + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::compute_mass, std::placeholders::_1)); + + // Assign beta to each particle + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::assign_projection_parameter, + std::placeholders::_1, beta_)); + + // Check point resume + if (resume) this->checkpoint_resume(); + + auto solver_begin = std::chrono::steady_clock::now(); + // Main loop + for (; step_ < nsteps_; ++step_) { + if (mpi_rank == 0) console_->info("Step: {} of {}.\n", step_, nsteps_); + +#pragma omp parallel sections + { + // Spawn a task for initialising nodes and cells +#pragma omp section + { + // Initialise nodes + mesh_->iterate_over_nodes( + std::bind(&mpm::NodeBase::initialise, std::placeholders::_1)); + + mesh_->iterate_over_cells( + std::bind(&mpm::Cell::activate_nodes, std::placeholders::_1)); + } + // Spawn a task for particles +#pragma omp section + { + // Iterate over each particle to compute shapefn + mesh_->iterate_over_particles(std::bind( + &mpm::ParticleBase::compute_shapefn, std::placeholders::_1)); + } + } // Wait to complete + + // Assign mass and momentum to nodes + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::map_mass_momentum_to_nodes, + std::placeholders::_1)); + + // Compute free surface cells, nodes, and particles + mesh_->compute_free_surface(free_surface_detection_, volume_tolerance_); + + // Spawn a task for initializing pressure at free surface +#pragma omp parallel sections + { +#pragma omp section + { + // Assign initial pressure for all free-surface particle + mesh_->iterate_over_particles_predicate( + std::bind(&mpm::ParticleBase::assign_pressure, + std::placeholders::_1, 0.0, fluid), + std::bind(&mpm::ParticleBase::free_surface, + std::placeholders::_1)); + } + } // Wait to complete + + // Compute nodal velocity at the begining of time step + mesh_->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::compute_velocity, + std::placeholders::_1), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + // Iterate over each particle to compute strain rate + mesh_->iterate_over_particles(std::bind( + &mpm::ParticleBase::compute_strain, std::placeholders::_1, dt_)); + + // Iterate over each particle to compute shear (deviatoric) stress + mesh_->iterate_over_particles(std::bind( + &mpm::ParticleBase::compute_stress, std::placeholders::_1)); + + // Spawn a task for external force +#pragma omp parallel sections + { +#pragma omp section + { + // Iterate over particles to compute nodal body force + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::map_body_force, + std::placeholders::_1, this->gravity_)); + + // Apply particle traction and map to nodes + mesh_->apply_traction_on_particles(this->step_ * this->dt_); + } + +#pragma omp section + { + // Spawn a task for internal force + // Iterate over each particle to compute nodal internal force + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::map_internal_force, + std::placeholders::_1)); + } + } // Wait for tasks to finish + + // Compute intermediate velocity + mesh_->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::compute_acceleration_velocity, + std::placeholders::_1, fluid, this->dt_), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + // Reinitialise system matrix to perform PPE + bool matrix_reinitialization_status = this->reinitialise_matrix(); + if (!matrix_reinitialization_status) { + status = false; + throw std::runtime_error("Reinitialisation of matrix failed"); + } + + // Compute poisson equation + this->compute_poisson_equation(); + + // Assign pressure to nodes + mesh_->iterate_over_nodes_predicate( + std::bind(&mpm::NodeBase::update_pressure_increment, + std::placeholders::_1, assembler_->pressure_increment(), + fluid, this->step_ * this->dt_), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + // Use nodal pressure to update particle pressure + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::compute_updated_pressure, + std::placeholders::_1)); + + // Compute correction force + this->compute_correction_force(); + + // Compute corrected acceleration and velocity + mesh_->iterate_over_nodes_predicate( + std::bind( + &mpm::NodeBase< + Tdim>::compute_acceleration_velocity_navierstokes_semi_implicit, + std::placeholders::_1, fluid, this->dt_), + std::bind(&mpm::NodeBase::status, std::placeholders::_1)); + + // Update particle position and kinematics + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::compute_updated_position, + std::placeholders::_1, this->dt_, velocity_update_)); + + // Apply particle velocity constraints + mesh_->apply_particle_velocity_constraints(); + + // Pressure smoothing + if (pressure_smoothing_) this->pressure_smoothing(fluid); + + // Locate particle + this->locate_particle(); + + if (step_ % output_steps_ == 0) { + // HDF5 outputs + this->write_hdf5(this->step_, this->nsteps_); +#ifdef USE_VTK + // VTK outputs + this->write_vtk(this->step_, this->nsteps_); +#endif + } + } + auto solver_end = std::chrono::steady_clock::now(); + console_->info("Rank {}, SemiImplicit_NavierStokes solver duration: {} ms", + mpi_rank, + std::chrono::duration_cast( + solver_end - solver_begin) + .count()); + + return status; +} + +// Semi-implicit functions +// Initialise matrix +template +bool mpm::MPMSemiImplicitNavierStokes::initialise_matrix() { + bool status = true; + try { + // Max iteration steps + unsigned max_iter = + analysis_["linear_solver"]["max_iter"].template get(); + // Tolerance + double tolerance = + analysis_["linear_solver"]["tolerance"].template get(); + // Get matrix assembler type + std::string assembler_type = analysis_["linear_solver"]["assembler_type"] + .template get(); + // Get matrix solver type + std::string solver_type = + analysis_["linear_solver"]["solver_type"].template get(); + // Create matrix assembler + assembler_ = + Factory>::instance()->create(assembler_type); + // Create matrix solver + linear_solver_ = + Factory>, unsigned, + double>::instance() + ->create(solver_type, std::move(max_iter), std::move(tolerance)); + // Assign mesh pointer to assembler + assembler_->assign_mesh_pointer(mesh_); + + // Get method to detect free surface detection + free_surface_detection_ = "density"; + if (analysis_["free_surface_detection"].contains("type")) + free_surface_detection_ = analysis_["free_surface_detection"]["type"] + .template get(); + // Get volume tolerance for free surface + volume_tolerance_ = analysis_["free_surface_detection"]["volume_tolerance"] + .template get(); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +// Reinitialise and resize matrices at the beginning of every time step +template +bool mpm::MPMSemiImplicitNavierStokes::reinitialise_matrix() { + bool status = true; + try { + // Assigning matrix id + const auto nactive_node = mesh_->assign_active_nodes_id(); + + // Assign global node indice + assembler_->assign_global_node_indices(nactive_node); + + // Assign pressure constraints + assembler_->assign_pressure_constraints(this->beta_, + this->step_ * this->dt_); + + // Initialise element matrix + mesh_->iterate_over_cells(std::bind( + &mpm::Cell::initialise_element_matrix, std::placeholders::_1)); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +// Compute poisson equation +template +bool mpm::MPMSemiImplicitNavierStokes::compute_poisson_equation( + std::string solver_type) { + bool status = true; + try { + // Construct local cell laplacian matrix + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::map_laplacian_to_cell, + std::placeholders::_1)); + + // Assemble global laplacian matrix + assembler_->assemble_laplacian_matrix(dt_); + + // Map Poisson RHS matrix + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::map_poisson_right_to_cell, + std::placeholders::_1)); + + // Assemble poisson RHS vector + assembler_->assemble_poisson_right(dt_); + + // Assign free surface to assembler + assembler_->assign_free_surface(mesh_->free_surface_nodes()); + + // Apply constraints + assembler_->apply_pressure_constraints(); + + // Solve matrix equation and assign solution to assembler + assembler_->assign_pressure_increment( + linear_solver_->solve(assembler_->laplacian_matrix(), + assembler_->poisson_rhs_vector(), solver_type)); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} + +//! Compute correction force +template +bool mpm::MPMSemiImplicitNavierStokes::compute_correction_force() { + bool status = true; + try { + // Map correction matrix from particles to cell + mesh_->iterate_over_particles( + std::bind(&mpm::ParticleBase::map_correction_matrix_to_cell, + std::placeholders::_1)); + + // Assemble correction matrix + assembler_->assemble_corrector_right(dt_); + + // Assign correction force + mesh_->compute_nodal_correction_force( + assembler_->correction_matrix(), assembler_->pressure_increment(), dt_); + + } catch (std::exception& exception) { + console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what()); + status = false; + } + return status; +} diff --git a/include/utilities/radial_basis_function.h b/include/utilities/radial_basis_function.h new file mode 100644 index 000000000..d872a2097 --- /dev/null +++ b/include/utilities/radial_basis_function.h @@ -0,0 +1,308 @@ +#ifndef MPM_RADIAL_BASIS_FUNCTION_H_ +#define MPM_RADIAL_BASIS_FUNCTION_H_ + +#include + +#include "Eigen/Dense" + +#include "logger.h" + +namespace mpm { + +// Namespace for radial basis function handling +namespace RadialBasisFunction { + +//! Cubic Spline Radial Basis Function +//! Source: Monaghan, 1985; Monaghan, 1992 +template +double cubic_spline(const double smoothing_length, const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 15.0 / (7.0 * M_PI * std::pow(smoothing_length, 2)); + else if (Tdim == 3) + multiplier = 3.0 / (2.0 * M_PI * std::pow(smoothing_length, 3)); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function + double basis_function = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 1.0) + basis_function *= + (2.0 / 3.0 - std::pow(radius, 2) + 0.5 * std::pow(radius, 3)); + else if (radius >= 1.0 && radius < 2.0) + basis_function *= (1.0 / 6.0 * std::pow((2.0 - radius), 3)); + else + basis_function = 0.0; + + return basis_function; +} + +//! Cubic Spline Radial Basis Function derivative +//! Source: Monaghan, 1985; Monaghan, 1992 +template +double cubic_spline_derivative(const double smoothing_length, + const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 15.0 / (7.0 * M_PI * std::pow(smoothing_length, 2)); + else if (Tdim == 3) + multiplier = 3.0 / (2.0 * M_PI * std::pow(smoothing_length, 3)); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function derivative + double dw_dr = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 1.0) + dw_dr *= (-2.0 * radius + 1.5 * std::pow(radius, 2)); + else if (radius >= 1.0 && radius < 2.0) + dw_dr *= -0.5 * std::pow((2.0 - radius), 2); + else + dw_dr = 0.0; + + return dw_dr; +} + +//! Quintic Spline Radial Basis Function +//! Source: Liu, 2010 +template +double quintic_spline(const double smoothing_length, + const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 1.0 / (478.0 * M_PI * std::pow(smoothing_length, 2)); + else if (Tdim == 3) + multiplier = 3.0 / (359.0 * M_PI * std::pow(smoothing_length, 3)); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function + double basis_function = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 1.0) + basis_function *= + (std::pow(3.0 - radius, 5) - 6.0 * std::pow(2.0 - radius, 5) + + 15.0 * std::pow(1.0 - radius, 5)); + else if (radius >= 1.0 && radius < 2.0) + basis_function *= + (std::pow(3.0 - radius, 5) - 6.0 * std::pow(2.0 - radius, 5)); + else if (radius >= 2.0 && radius < 3.0) + basis_function *= (std::pow(3.0 - radius, 5)); + else + basis_function = 0.0; + + return basis_function; +} + +//! Quintic Spline Radial Basis Function derivative +//! Source: Liu, 2010 +template +double quintic_spline_derivative(const double smoothing_length, + const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 1.0 / (478.0 * M_PI * std::pow(smoothing_length, 2)); + else if (Tdim == 3) + multiplier = 3.0 / (359.0 * M_PI * std::pow(smoothing_length, 3)); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function + double dw_dr = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 1.0) + dw_dr *= + (-5.0 * std::pow(3.0 - radius, 4) + 30. * std::pow(2.0 - radius, 4) - + 75. * std::pow(1.0 - radius, 4)); + else if (radius >= 1.0 && radius < 2.0) + dw_dr *= + (-5.0 * std::pow(3.0 - radius, 4) + 30. * std::pow(2.0 - radius, 4)); + else if (radius >= 2.0 && radius < 3.0) + dw_dr *= (-5.0 * std::pow(3.0 - radius, 4)); + else + dw_dr = 0.0; + + return dw_dr; +} + +//! Gaussian Kernel +//! Source: Liu, 2010 +template +double gaussian(const double smoothing_length, const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 1.0 / (M_PI * std::pow(smoothing_length, 2)); + else if (Tdim == 3) + multiplier = 1.0 / std::pow((std::sqrt(M_PI) * smoothing_length), 3); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function + double basis_function = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 3.0) + basis_function *= std::exp(-std::pow(radius, 2)); + else + basis_function = 0.0; + + return basis_function; +} + +//! Gaussian Kernel derivative +//! Source: Liu, 2010 +template +double gaussian_derivative(const double smoothing_length, + const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 1.0 / (M_PI * std::pow(smoothing_length, 2)); + else if (Tdim == 3) + multiplier = 1.0 / std::pow((std::sqrt(M_PI) * smoothing_length), 3); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function + double dw_dr = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 3.0) + dw_dr *= -2.0 * radius * std::exp(-std::pow(radius, 2)); + else + dw_dr = 0.0; + + return dw_dr; +} + +//! Super Gaussian Kernel +//! Source: Monaghan, 1992 +template +double super_gaussian(const double smoothing_length, + const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 1.0 / std::pow((std::sqrt(M_PI) * smoothing_length), 2); + else if (Tdim == 3) + multiplier = 1.0 / std::pow((std::sqrt(M_PI) * smoothing_length), 3); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function + double basis_function = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 3.0) + basis_function *= radius * (2.0 * radius * radius - (double)Tdim - 4.0) * + std::exp(-std::pow(radius, 2)); + else + basis_function = 0.0; + + return basis_function; +} + +//! Super Gaussian Kernel derivative +//! Source: Monaghan, 1992 +template +double super_gaussian_derivative(const double smoothing_length, + const double norm_distance) { + + // Assign multiplier depends on dimension + double multiplier; + if (Tdim == 2) + multiplier = 1.0 / std::pow((std::sqrt(M_PI) * smoothing_length), 2); + else if (Tdim == 3) + multiplier = 1.0 / std::pow((std::sqrt(M_PI) * smoothing_length), 3); + else + throw std::runtime_error("Tdim is invalid"); + + // Compute basis function + double dw_dr = multiplier; + const double radius = norm_distance / smoothing_length; + if (radius >= 0.0 && radius < 3.0) + dw_dr *= std::exp(-std::pow(radius, 2)) * + ((double)Tdim / 2.0 + 1.0 - radius * radius); + else + dw_dr = 0.0; + + return dw_dr; +} + +//! General Radial Basis Function Kernel call +template +double kernel(const double smoothing_length, const double norm_distance, + const std::string type = "cubic_spline") { + if (type == "cubic_spline") { + return mpm::RadialBasisFunction::cubic_spline(smoothing_length, + norm_distance); + } else if (type == "quintic_spline") { + return mpm::RadialBasisFunction::quintic_spline(smoothing_length, + norm_distance); + } else if (type == "gaussian") { + return mpm::RadialBasisFunction::gaussian(smoothing_length, + norm_distance); + } else if (type == "super_gaussian") { + return mpm::RadialBasisFunction::super_gaussian(smoothing_length, + norm_distance); + } else { + throw std::runtime_error( + "RadialBasisFunction kernel type is invalid. Available types are: " + "\"cubic_spline\", \"quintic_spline\", \"gaussian\", and, " + "\"super_gaussian\"."); + } +} + +//! General Radial Basis Function Kernel call +template +Eigen::Matrix gradient( + const double smoothing_length, + const Eigen::Matrix& relative_distance, + const std::string type = "cubic_spline") { + + // Compute norm distance + const double norm_distance = relative_distance.norm(); + double dw_dr; + if (type == "cubic_spline") { + dw_dr = mpm::RadialBasisFunction::cubic_spline_derivative( + smoothing_length, norm_distance); + } else if (type == "quintic_spline") { + dw_dr = mpm::RadialBasisFunction::quintic_spline_derivative( + smoothing_length, norm_distance); + } else if (type == "gaussian") { + dw_dr = mpm::RadialBasisFunction::gaussian_derivative( + smoothing_length, norm_distance); + } else if (type == "super_gaussian") { + dw_dr = mpm::RadialBasisFunction::super_gaussian_derivative( + smoothing_length, norm_distance); + } else { + throw std::runtime_error( + "RadialBasisFunction gradient type is invalid. Available types are: " + "\"cubic_spline\", \"quintic_spline\", \"gaussian\", and, " + "\"super_gaussian\"."); + } + + // Gradient = dw_dr * r / ||r|| / h + Eigen::Matrix gradient = relative_distance; + if (norm_distance > 1.e-12) + gradient *= dw_dr / (norm_distance * smoothing_length); + else + gradient *= 0.0; + + return gradient; +} + +} // namespace RadialBasisFunction + +} // namespace mpm +#endif \ No newline at end of file diff --git a/src/io/logger.cc b/src/io/logger.cc index 82e2b262c..26b18caeb 100644 --- a/src/io/logger.cc +++ b/src/io/logger.cc @@ -35,3 +35,8 @@ const std::shared_ptr mpm::Logger::mpm_explicit_usf_logger = // Create a logger for MPM Explicit USL const std::shared_ptr mpm::Logger::mpm_explicit_usl_logger = spdlog::stdout_color_st("MPMExplicitUSL"); + +// Create a logger for MPM Semi Implicit Two Phase +const std::shared_ptr + mpm::Logger::mpm_semi_implicit_navier_stokes_logger = + spdlog::stdout_color_st("MPMSemiImplicitNavierStokes"); \ No newline at end of file diff --git a/src/linear_solver.cc b/src/linear_solver.cc new file mode 100644 index 000000000..4314689e2 --- /dev/null +++ b/src/linear_solver.cc @@ -0,0 +1,23 @@ +#include "assembler_base.h" +#include "assembler_eigen_semi_implicit_navierstokes.h" + +#include "cg_eigen.h" +#include "solver_base.h" + +// Assembler collections +// Asssembler 2D +static Register, + mpm::AssemblerEigenSemiImplicitNavierStokes<2>> + assembler_eigen_semi_implicit_navierstokes_2d( + "EigenSemiImplicitNavierStokes2D"); +// Asssembler 3D +static Register, + mpm::AssemblerEigenSemiImplicitNavierStokes<3>> + assembler_eigen_semi_implicit_navierstokes_3d( + "EigenSemiImplicitNavierStokes3D"); + +// Linear Solver collections +// Eigen Conjugate Gradient +static Register>, + mpm::CGEigen>, unsigned, double> + solver_eigen_cg("CGEigen"); diff --git a/src/mpm.cc b/src/mpm.cc index 19b3f5dc9..8327f959d 100644 --- a/src/mpm.cc +++ b/src/mpm.cc @@ -4,6 +4,7 @@ #include "io.h" #include "mpm.h" #include "mpm_explicit.h" +#include "mpm_semi_implicit_navierstokes.h" namespace mpm { // Stress update method @@ -20,3 +21,13 @@ static Register, const std::shared_ptr&> // 3D Explicit MPM static Register, const std::shared_ptr&> mpm_explicit_3d("MPMExplicit3D"); + +// 2D SemiImplicit Navier Stokes MPM +static Register, + const std::shared_ptr&> + mpm_semi_implicit_navierstokes_2d("MPMSemiImplicitNavierStokes2D"); + +// 3D SemiImplicit Navier Stokes MPM +static Register, + const std::shared_ptr&> + mpm_semi_implicit_navierstokes_3d("MPMSemiImplicitNavierStokes3D"); \ No newline at end of file diff --git a/src/particle.cc b/src/particle.cc index 7488c8023..3ae854cb5 100644 --- a/src/particle.cc +++ b/src/particle.cc @@ -1,5 +1,6 @@ #include "particle.h" #include "factory.h" +#include "fluid_particle.h" #include "particle_base.h" // Particle2D (2 Dim) @@ -11,3 +12,13 @@ static Register, mpm::Particle<2>, mpm::Index, static Register, mpm::Particle<3>, mpm::Index, const Eigen::Matrix&> particle3d("P3D"); + +// Single phase (fluid) particle2D (2 Dim) +static Register, mpm::FluidParticle<2>, mpm::Index, + const Eigen::Matrix&> + particle2dfluid("P2DFLUID"); + +// Single phase (fluid) particle3D (3 Dim) +static Register, mpm::FluidParticle<3>, mpm::Index, + const Eigen::Matrix&> + particle3dfluid("P3DFLUID"); \ No newline at end of file diff --git a/tests/particle_test.cc b/tests/particle_test.cc index 0d2656f54..1e80c0a7a 100644 --- a/tests/particle_test.cc +++ b/tests/particle_test.cc @@ -1137,9 +1137,15 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { SECTION("Assign state variables") { // Assign material properties REQUIRE(particle->assign_material(mc_material) == true); - // Assing state variables + // Assign state variables REQUIRE(particle->assign_material_state_vars(state_variables, mc_material) == true); + // Assign and read a state variable + REQUIRE_NOTHROW(particle->assign_state_variable("phi", 30.)); + REQUIRE(particle->state_variable("phi") == 30.); + // Assign and read pressure though MC does not contain pressure + REQUIRE_NOTHROW(particle->assign_pressure(1000)); + REQUIRE(std::isnan(particle->pressure()) == true); } SECTION("Assign state variables fail on state variables size") { @@ -1157,7 +1163,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // Assign material properties REQUIRE(particle->assign_material(newtonian_material) == true); - // Assing state variables + // Assign state variables REQUIRE(particle->assign_material_state_vars(state_variables, mc_material) == false); } @@ -1177,7 +1183,7 @@ TEST_CASE("Particle is checked for 2D case", "[particle][2D]") { // Assign material properties REQUIRE(particle->assign_material(newtonian_material) == true); - // Assing state variables + // Assign state variables REQUIRE(particle->assign_material_state_vars(state_variables, mc_material) == false); } @@ -2389,9 +2395,15 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { SECTION("Assign state variables") { // Assign material properties REQUIRE(particle->assign_material(mc_material) == true); - // Assing state variables + // Assign state variables REQUIRE(particle->assign_material_state_vars(state_variables, mc_material) == true); + // Assign and read a state variable + REQUIRE_NOTHROW(particle->assign_state_variable("phi", 30.)); + REQUIRE(particle->state_variable("phi") == 30.); + // Assign and read pressure though MC does not contain pressure + REQUIRE_NOTHROW(particle->assign_pressure(1000)); + REQUIRE(std::isnan(particle->pressure()) == true); } SECTION("Assign state variables fail on state variables size") { @@ -2409,7 +2421,7 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { // Assign material properties REQUIRE(particle->assign_material(newtonian_material) == true); - // Assing state variables + // Assign state variables REQUIRE(particle->assign_material_state_vars(state_variables, mc_material) == false); } @@ -2429,7 +2441,7 @@ TEST_CASE("Particle is checked for 3D case", "[particle][3D]") { // Assign material properties REQUIRE(particle->assign_material(newtonian_material) == true); - // Assing state variables + // Assign state variables REQUIRE(particle->assign_material_state_vars(state_variables, mc_material) == false); }