Semi-Implicit Navier-Stokes solver for incompressible flow

[bodhinandach]

Summary

This RFC is to propose the Navier-Stokes solver in branch solver/navier-stokes. The NS solver utilized semi-implicit Chorin's projection scheme which is mainly to be used for modeling incompressible fluids.

Motivation

To extend the capability of treating internal incompressibility constraint more accurately particularly for simulations of fluid flows. The proposed method has been commonly studied in CFD community and proven to be one of the simplest yet efficient.

Design Detail

The following classes and functionality were implemented:

1. MPMSemiImplicitNavierStokes solver class (this at the moment only works with point 2.)
2. A FluidParticle derived from Particle with some specialized class derivation. (Will be refactor following #637, #654, #673)
3. Linear solvers and assembler classes based on Eigen.

a. The assembler class is made to assemble the global system of equations. First, all nodes will be assigned an active_id to specify the index of each node in the global systems of equations; this id is changing every time step as the particles move over the background grids. Then, when the explicit step is completed, a laplacian and RHS matrix are constructed in all cells. We store this local matrix with size of (nnodes x nnodes) as Eigen matrices in Cell along with the correction_matrix as:

  //! Local laplacian matrix
  Eigen::MatrixXd laplacian_matrix_;
  //! Local poisson RHS matrix
  Eigen::MatrixXd poisson_right_matrix_;
  //! Local correction RHS matrix
  Eigen::MatrixXd correction_matrix_;

These local matrices are by default has zero sizes, and only initialized when the initialise_element_matrix function is called by the solver. Each time set, the assembler job is to organize all the cell's local stiffness matrix to a global matrix and RHS vector and to apply constraints to the assembled systems through row-column modification.

b. The linear solver solves the given assembled coefficient matrix and RHS vector, and it is independent of dimension. The solver_type should be specified in the input .json.

virtual Eigen::VectorXd solve(const Eigen::SparseMatrix<double>& A,
                                const Eigen::VectorXd& b,
                                std::string solver_type) = 0;

c. Parallel solver capability is also added to improve the solver scalability and efficiency as previously noted in RFC #635.

4. Free-surface detection algorithms (both cheap and expensive detection were implemented) in Mesh class.
5. Some boundary conditions were implemented, e.g. pressure BC.
6. Some work on level-set signed distance computation if accurate free-surface detection method is considered.
7. Few testings have been added and few more will be added.

Drawbacks

No drawbacks in performance at the moment. It will enrich the capability of the software in many ways. However, the implementation might cause the function derivatives to be a bit messy, particularly in Cell class where we need to store and construct some local element matrices. Further thoughts and suggestions are welcome.

Rationale and Alternatives

After the merge of the NS solver, we can also extend the solver to be used for any incompressible Solid with minimum modification. Some further thoughts on including purely Implicit solver class is also possible for future development. They utilize very much the same structure and the capability can be extended from the proposed work.

Prior Art

Referring a lot to the work of @srhgk2:

1. Kularathna, S., & Soga, K. (2017). Implicit formulation of material point method for analysis of incompressible materials. Computer Methods in Applied Mechanics and Engineering, 313, 673-686.

Changelog

1. Add details of linear solver and assembler class structure.
2. Add capability of parallel solver and other MPI processes in fluid particle.

[kks32]

Related to the discussion in #633 how do we handle HDF5 data, and would we have mixed types of particles in the mesh container?

[bodhinandach]

Handling HDF5 data is a bit difficult for two-phase. I think the one that we implemented in TwoPhaseParticle is not the best way. I am thinking that it's good to find a way to derive ParticleHDF5 structure, we didn't manage to do it last time though. Right @tianchiTJ?

For FluidParticle, however, we can use the normal one cause there are no extra parameters needed. Mostly all the variables are the same with the normal Particle class.

[bodhinandach]

@kks32 for the particle in mesh container, we don't have to worry, since it is derived from the same base class, i.e. ParticleBase, it should be the same as it is now.

[kks32]

The container will be fine, but you'll be iterating through all types of particles in the same container, which will be inefficient. For example,

  // Assign beta to each particle
  mesh_->iterate_over_particles(
      std::bind(&mpm::ParticleBase<Tdim>::assign_projection_parameter,
                std::placeholders::_1, beta_));

This has to be filtered to run only for fluid particles

[bodhinandach]

You are right. The function assign_projection_parameter should be run in fluid as the current semi-implicit code only works for FluidParticle.

[kks32]

Please add code snippets and class outlines to the RFC, it should give an overview and a lot of details of implementation especially on the linear solvers, free-surface detection, and boundary conditions.

[stale[bot]]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

[stale[bot]]

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.