Meshfree and Particle Methods

Particle Methods

Particle methods approximate flow equations on points that move with the flow, and thus solve advection exactly. This is a great advantage over Eulerian approaches, for which an accurate treatment of advection is still subject to improvement. In addition, particle methods have many more advantages, for instance can diffusion in high dimensional problems be resembled efficiently. The main focus in my research is their Lagrangian property, i.e. the exact solution of advection. Here advection can refer to the material derivative ∂t+(u⋅∇) in fluid flow, or the derivative of the flux function f'(u) in a conservation law u_t+(f(u))x = 0.

The advantages of Lagrangian particles come at a price:

The aspect of particle management becomes particularly interesting for scalar conservation laws in one space dimension u_t+(f(u))x = 0. The method particleclaw approximates the solution using only characteristic particle movement, and particle merges to represent shocks.

Another interesting application is the numerical simulation of macroscopic traffic models. Particles moving with the local flow velocity resemble the behavior of cars, although the total number of particles need not equal the total number of cars. This connection can be used to obtain insight into the macroscopic limit of microscopic traffic models and vice versa.

Meshfree Finite Difference Methods

Meshfree methods approximate a partial differential equations on an unstructured set of points without any specified connections between points (a point cloud). Meshfree methods exist in form of finite element, finite volume, and finite difference approaches. The former two are commonly defined via basis functions that are constructed by a partition of unity on the point cloud. The latter can be defined by Taylor expansion. In my research, I focus of meshfree finite difference approaches. On regular grids, finite difference methods are very successful and very flexible. The goal is to carry their advantages over to the unstructured case.

Meshfree methods operate on fully unstructured data sets (as opposed to immersed interface, ghost fluid, or level set methods). Compared to unstructured mesh-based approaches, they avoid the task of construction and update of a mesh. Therefore, they are particularly suited for Lagrangian approaches of flow equations, that move the data structures with the flow velocity (such as particle methods). The advantages of meshfree methods come at a price. Two fundamental difficulties are:

An important problem is the approximation of the Laplace operator, for instance for the treatment of diffusion and in a pressure correction step for the computation of incompressible flows. For a given point cloud, at each interior point, the Laplace operator shall be approximated by a linear combination of function values of nearby points. The coefficients are gathered in a finite difference stencil. By Taylor expansion, one can derive algebraic conditions that a consistent meshfree approximation of the Laplacian has to satisfy. These are {3,6,10} conditions in {1D,2D,3D}. Therefore, for each point, at least {2,5,9} neighboring points have to be selected. If more neighbors are considered, many stencils exist, and a unique one can to be defined by a minimization problem. A desirable property for a Laplace stencil is positivity: The central coefficient is negative, and all neighboring coefficients are positive. This property guarantees flow of heat from warm to cold, and yields M-matrices in Poisson problems.

Generally, a meshfree finite difference approximation for a given point consists of two steps:

  1. Select neighbors. This is commonly done by geometric criteria. The figures below show examples of used constructions.
  2. Compute a stencil, defined by a minimization problem. Commonly used is a weighted least squares formulation. One fundamental problem with this approach is that it may yield stencils that are not positive. In the figures below, red points are selected, and blue markers indicate negative coefficients obtained by a least squares approach.

Circular neighborhood Voronoi neighborhood Quadrant neighborhood MPS neighborhood
Circular neighborhood. All points inside a given radius are considered.
Problems: For locally unbalanced point clouds, neighbors may not lie around the central point. A large number of neighbors may be selected. No guarantee for positive stencil. 		Voronoi neighborhood. Neighbors are defined by Delaunay edges.
Problems: Only four neighbors may be selected, but five are required. The computation of the Delaunay neighborhood can be costly. No guarantee for positive stencil. 		Quadrant neighborhood. In each quadrant (octant in 3D), the two closest points are selected.
Problems: Artificial dependence on local coordinates. No guarantee for positive stencil. 		MPS neighborhood. Comes out of a linear minimization approach. Stencil is minimal and positive by construction.

A new approach, introduced and analyzed in the references below, is the minimal positive stencil (MPS) method. For each interior point, all points inside a circular radius are considered. For these, a linear minimization problem with sign constraints (enforcing positive stencils) is formulated. Let s denote the desired stencil vector (without the central entry), then the MPS stencil is the solution of the linear program

min ∑ si/wi s.t. V⋅s = b, s≥ 0 .
CPU times for LP solvers

Here V⋅s = b is the linear system of constraints, and s≥ 0 indicates that each component is non-negative. The solution is positive due to the sign constraints, and minimal due to the fundamental theorem of linear programming. The approach falls into the context of L1 minimization, that in recent years has received significant attention in statistics. Efficient methods for the linear programs are still subject to research. Tests with CPLEX indicate that simplex methods generally are the fastest.

The MPS approach approximates −Δ by an M-matrix, which is an advantageous structure for many linear solvers. In particular, algebraic multigrid (AMG) methods perform very well on the arising linear systems. Numerical tests on a 3D geometry appearing in an industrial application were performed. While the linear minimization approach is slightly slower than classical least squares approaches, the optimal sparsity reduces the overall solution time significantly. In addition, less memory is required.

The computational performance of the new approach is measured in comparison with a classical least squares approach. CPU times and memory consumption are measured for a 3D Poisson problem in a complex geometry, for a sequence of resolutions. In the figures below, the blue curves represent a least squares approach, while the red curves show the results for the linear minimization approach. The MPS approach selects a minimal stencil out of a larger candidate set of neighbors. The dashed red curves indicate solution times with the whole candidate set stored in the sparse matrix. This approach unnecessarily stores many zero entries, and thus shows that the speed-up in the MPS method is due to the optimal sparsity.

CPU times matrix setup Memory consumption CPU times solve by BiCGstab CPU times solve by AMG
CPU times for setup of system matrix Memory consumption CPU times for solution with BiCGstab CPU times for solution with AMG

Related Publications

B. Seibold, Performance of algebraic multigrid methods for non-symmetric matrices arising in particle methods, Numer. Linear Algebra Appl., Vol. 17, 2010, pp. 433-451.
B. Seibold, Minimal positive stencils in meshfree finite difference methods for the Poisson equation, Comput. Methods Appl. Mech. Engrg., Vol. 198, 2008, pp. 592-601.
B. Seibold, M-Matrices in meshless finite difference methods, PhD thesis, Department of Mathematics, University of Kaiserslautern, 2006.