Geometric Multigrid Solver for Unstructured Finite Elements

Recently I've spent good chunk of time implementing and studying the efficacy of utilizing geometric Multigrid (GMG) for solving partial differential equations (PDEs) on unstructured grids distributed over arbitrary processor counts.

Geometric Multigrid is regarded as a particularly effective solution technique; it's algorithmically optimal ($O(n)$) in both memory requirements as well as computational cost and it naturally complements iterative solution techniques from the realm of linear algebra as it can be understood as a preconditioning strategy. Furthermore, if one employs load balanced domain decomposition techniques GMG retains its optimal properties under parallelization. The overall methodology operates by smoothing the solution errors (or residuals) on a collection of coarser sub-problems whose information can then be used to repeatedly solve the originally desired fine grid problem in a highly efficient fashion. The method is algorithmically scaleable, and it converges to the solution in a fixed number of iterations that is independent of the mesh resolution. A key aspect and difficulty of such a technique is the development of an implementation in an efficient, parallel manner which is independent of the input parameters to the PDE which is to be solved. As such, much of the solver infrastructure related to this research is implemented in the open source libMesh finite element library, while examples of the optimality of the technique for more complex examples were demonstrated in the GRINS multiphysics application built on top of libMesh.

In order to better understand the method we can tesselate our PDE domain and define a sequence of grids denoted $\{G_i\}$, and denote by $\#G_i$ the number of mesh elements within a given grid $G_i$. For simplicity we assume that this grid sequence is ordered based on the count of the finite elements so that we have: $\#G_f >... > \#G_i > ... > \#G_c$ where we adopt the notation so that $G_f$ and $G_c$ correspond respectively to the finest and coarsest grids. This grid sequences then serves as a fundamental component to the Multigrid technique.

As such problems are often too large to be solved directly and on one computer, the libMesh library is then used to partition the domain (i.e. each grid) between processors, refine the mesh, and then manage local-to-global degree of freedom information as well as ghosted element info which must be periodically communicated between processors. The following diagram attempts to depict these steps in an abstracted manner.

Upon establishing the grid hierarchy the Multigrid technique utilizes the degree of freedom map managed by libMesh to construct projection and restriction operators which become fundamental pieces allowing the algorithm to transfer information between grids.

As a simple example of the solver in action we can consider Poisson's equation, the prototypical elliptic PDE which models diffusive processes. Poisson's equation has a strong formulation given by:

$$\nabla^2 \varphi = f$$

where $f(x,y,z)$ is some specified function and $\varphi$ is sought. For a simple scenario we can let $\varphi$ represent temperature, and select a two dimensional square domain $\Gamma$ with insulated boundaries alongside an internal forcing with a heat source given by:

$$f(x,y,z) = 2\pi^2 \sin(\pi x)\sin(\pi y)$$

This gives rise to a radially symmetric solution as depicted in the following image:

Utilizing such a basic configuration allows us to inspect the convergence rates by comparing our computed solution to an exact solution via the method of manufactured solutions. We know from a priori finite element error estimates that we can expect the error in a Poisson problem to decay like $O(h^2)$, where $h$ denotes the grid spacing. On the other hand Multigrid convergence theory dictates that when utilizing the full Multigrid method (FMG) one can expect an $O(n)$ convergence rate. With this information we can construct an experiment where we monitor the error decay while solving the problem on successively smaller grid resolutions (via uniform refinement of the $N$ finite elements) each time utilizing exactly one application of a Multigrid iteration; we compare the suboptimal V-cycle convergence to the optimal FMG technique in the following plot:

As expected, one application of FMG produces a solution on the order of the discretization error while the V-cycle convergence deteriorates after increasing the problem size to $\sim 10^4$ elements.

Cleverly constructing the solver infrastructure in an element-agnostic and parallelizable manner allows us to leverage Multigrid in order to study more exotic physics at a high resolution. For example we can consider thermally driven fluid flow which incorporates a Boussinesq buoyancy approximation that couples the energy balance equation to that of the incompressible Navier Stokes equations.

$$ \begin{aligned} \frac{\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \cdot \nabla)\boldsymbol{u} &- \mu \nabla^2 \boldsymbol{u} + \frac{1}{\rho}\nabla p = - \beta \textbf{g} T \\ \nabla \cdot \boldsymbol{u} &= 0 \\ \frac{\partial T}{\partial t} + \boldsymbol{u}\cdot \nabla T &- \nabla \cdot (k\nabla T) = \frac{J}{\rho c_p} \end{aligned} $$

In this formulation we seek to solve for the velocity of the fluid $\boldsymbol{u}$, the pressure $p$, and temperature $T$ for all time $t \in (0,\tau)$. Here $\mu$ represents the kinematic viscosity, $\textbf{g}$ gravity forces,$\rho$ the density, while $\beta$ corresponds to the thermal expansion coefficient driving the buoyancy, $c_p$ the specific heat at constant pressure, $k$ is the thermal diffusivity constant, and $J$ is the rate per unit volume of internal heat production, which for our purposes will be assumed to be zero. For the sake of continuity from previous examples we retain the domain of the unit square, but now impose no slip boundary conditions for the velocity, adiabatic boundary conditions for the temperature on a pair of opposing walls, and assign Dirichlet temperature boundary conditions on the remaining temperature boundaries so that we have:

$$ \begin{aligned} \boldsymbol{u} &=0 \text{ on } \Gamma \times (0,\tau) \\ T &= T_D \text{ on } \Gamma_D \times (0,\tau) \\ \nabla T \cdot \textbf{n} &= 0 \text{ on } \Gamma_N \times (0,\tau) \end{aligned} $$

We solve this nonlinear system of equations using Newtons Method while discretizing in the usual finite element manner. Proceeding as such leads to the large linear system shown below, a system which must often be solved in parallel due to high dimensional nature of the of the linear problem.

$$ \begin{equation} \begin{bmatrix} F & B & M_1 \\ B^{\top} & 0 & 0 \\ M_2 & 0 & K \end{bmatrix} \begin{bmatrix} \boldsymbol{u} \\ p \\ T \end{bmatrix} = \begin{bmatrix} f_1 \\ f_2 \\ f_3 \end{bmatrix} \end{equation} $$

This system admits a block structure where the $F,B$ and $B^\top$ matrixes are the Jacobian components for a pure Navier-Stokes formulation, while $M_1$, $M_2$ and $K$ blocks arise from the temperature contribution Jacobian terms. With proper grouping of terms and variables one can notice that both the $(0,0)$ and $(2,2)$ block appear elliptic for moderate parameter values. Thus these blocks become natural candidates for the application of Multigrid.

Using Multigrid as a preconditioner to the appropriate components of the linear system in conjunction with a Krylov accelerator allows for optimal time-to-solutions and is achieved in a manner which is independent of input configurations and parameters. In the figure below we present the solution to the steady state case from above where can observe the formation of thermal vortexes, a phenomena which is particularly difficult for standard black box solvers to resolve.

As the solver infrastructure is quite modular we can leave the realm of flow based physics and also apply the Multigrid technique to elastodynamics problems. For example consider the system below:

$$ \begin{aligned} \rho_0 \ddot{\boldsymbol{\phi}} = \nabla^\textbf{X} \cdot [\textbf{P(F)}] &+ \rho_0 \textbf{b}_m \quad \text{ in } B\times (0,\tau) \\ \boldsymbol{\phi} &= \textbf{g} \quad \text{ in } \Gamma_d \times (0,\tau) \\ \textbf{P(F)}\textbf{N} &= \textbf{h} \quad \text{ in } \Gamma_n \times (0,\tau) \\ \boldsymbol{\phi} (\cdot,0) &= \textbf{X} \quad \text{ in } B \\ \dot{\boldsymbol{\phi}} (\cdot,0) &= \textbf{V}_0 \quad \text{ in } B \end{aligned} $$

Here we seek to solve for the motion $\phi$ where in the above formulation $B$ represents the reference configuration while $\textbf{X}$ corresponds to the Lagrangian coordinate system. Furthermore, $\textbf{P(F)}$ is the first Piola-Kirchoff frame indifferent stress response function having $\textbf{F} = \nabla^\textbf{X}\boldsymbol{\phi}$, $\rho_0$ is the reference mass density, $\textbf{b}_m$ is the material description of a specified spatial body force field per unit mass, $\textbf{N}$ is the outward unit normal field on $\partial B$, and $\textbf{g,h,V}_0$ are assigned fields.

If we select an initially near flat circular mesh and solve the above system using a pseudo-transient continuation method with pressure acting normal to the surface we can understand the physics as modeling the inflation of a balloon or bubble. The initial and final configurations of such a simulation are depicted below.

By now we have hopefully come to see that the implemented solver infrastructure is quite flexible and serves as a powerful and robust computational technique. Further possibilities abound for applications of such an infrastructure. Through the introduction of different domains, boundary conditions, and couplings between various physics one can explore a large swath of problems stemming the realm of PDEs. Having the ability to explore such scenarios from the command line in a fast, parallel, modular and reusable manner provides a particularly effective method for modeling and investigating the behavior and dynamics of various problems of interest.