PetIBM: toolbox and applications of the immersed-boundary method on distributed-memory architectures

In the IBM framework, a collection of Lagrangian markers defines the immersed boundary (where boundary conditions are enforced) and the fluid equations are solved over the extended domain (including the body domain). The Eulerian mesh remains unmodified when computing the flow around multiple moving immersed bodies, which removes the need for remeshing at every time step. PetIBM discretizes the fluid equations using a second-order finite-difference scheme, various optional time-integrators, and a fully discrete projection method (Perot (1993)). It implements two immersed-boundary algorithms: the immersed-boundary projection method (Taira and Colonius (2007)) and its decoupled version (Li et al. (2016)).

In the IBM framework, a collection of Lagrangian markers defines the immersed boundary (where boundary conditions are enforced) and the fluid equations are solved over the extended domain (including the body domain).The Eulerian mesh remains unmodified when computing the flow around multiple moving immersed bodies, which removes the need for remeshing at every time step.PetIBM discretizes the fluid equations using a second-order finite-difference scheme, various optional time-integrators, and a fully discrete projection method (Perot (1993)).It implements two immersed-boundary algorithms: the immersed-boundary projection method (Taira and Colonius (2007)) and its decoupled version (Li et al. (2016)).
Other open-source software packages offer immersed-boundary solvers: for example, IBAMR (Griffith et al. (2007), Bhalla et al. (2013)) is a long-standing C++ library with MPI parallelization that also provides adaptive mesh refinement.It can handle deforming immersed bodies and has been used in a variety of scenarios, including cardiac fluid dynamics, swimming, insect flight, and others.PetIBM and IBAMR use different immersed-boundary schemes, however.We developed PetIBM to work with the immersed-boundary projection method, which is based on the fully discrete formulation of Perot on staggered grids and thus eliminates the need for pressure boundary conditions, which have caused many headaches for CFD practitioners (Gresho and Sani (1987), Sani and Gresho (1994)).PetIBM features an operator-based design, providing routines to create and manipulate discrete operators (e.g., gradient, divergence, Laplacian, convection, diffusion, etc.), so it can be used as a toolbox for researching new solution methods.It is also capable of using graphics processing unit (GPU) architectures, a feature missing from other software, as far as we know.A previous project implementing immersed-boundary methods on GPU architecture is cuIBM (Krishnan, Mesnard, and Barba (2017)), but it is limited to two-dimensional problems that fit on a single GPU device.
PetIBM is written in C++ and relies on the PETSc library (Balay et al. (1997), Balay et al. (2017)) for data structures and parallel routines to run on memory-distributed architectures.PetIBM can solve one or several linear systems on multiple distributed CUDAcapable GPU devices with the NVIDIA AmgX library and AmgXWrapper (Chuang and Barba ( 2017)).The software package includes extended documentation as well as many examples to guide users.
PetIBM has already been used to generate results published in Mesnard and Barba (2017), a full replication of a study on the aerodynamics of a gliding snake species (Krishnan et al. (2014)).PetIBM is currently used to compute the three-dimensional flow of a glidingsnake model on the cloud platform Microsoft Azure.

Appendix: mathematical formulation
PetIBM solves the Navier-Stokes equations on an extended discretization grid that includes the interior of the immersed boundary.To model the presence of the boundary, a forcing term is added to the momentum equation and an additional equation for the noslip condition completes the system.Variants of the immersed-boundary method (IBM) depend on how one models the forcing.In PetIBM, we use regularized-delta functions to transfer data between the Eulerian grid and the Lagrangian boundary points.The system of equations is: where u is the velocity field, p is the fluid pressure, and Re is the Reynolds number.
Currently, PetIBM provides two application codes implementing different versions of the IBM: (1) an immersed-boundary projection method (IBPM) based on the work of Taira and Colonius ( 2007) and ( 2) a decoupled version of the IBPM proposed by Li et al. (2016).
Those two methods fit into the framework of the projection approach of Perot (1993).The equations are fully discretized (space and time) to form an algebraic system to be solved for the velocity u n+1 , the pressure field ϕ, and the Lagrangian forces f .The discretized system is: where D, G, and A are the divergence, gradient, and implicit operators, respectively.E and H are the interpolation and spreading operators, respectively, used to transfer the data between the Eulerian grid and the Lagrangian boundary points.On the righthand side, r n gathers all the explicit terms and u n+1 B is the known (prescribed) boundary velocity; bc 1 and bc 2 contain the boundary terms that arise from the discretization of momentum and continuity equations, respectively.
In the IBPM, we solve a modified Poisson system for the pressure field and Lagrangian forces, coupled together.This way, the divergence-free condition and no-slip constraint are simultaneously enforced on the velocity field at the end of the time step.The fully discretized system can be cast into the following: with In practice, we never form the full system.Instead, we apply a block-LU decomposition as follows: Thus, we retrieve the sequence of operations of the traditional projection method.We solve a system for an intermediate velocity field that is corrected, after solving a modified Poisson system for the variable λ, to enforce the divergence-free condition and the no-slip constraint at the location of the immersed boundary.The sequence is: The IBPM implemented in PetIBM solves, at every time step, Equations ( 5) to ( 6).(Note: the inverse of the implicit operator A −1 is approximated by a finite Taylor series expansion.) The IBPM requires solving, at each time step, an expensive modified Poisson system, Q 1 A −1 Q 2 , whose non-zero structure changes when the location of the immersed boundary is moving.In the PetIBM implementation of the decoupled IBPM, we apply a second block-LU decomposition to decouple the pressure field from the Lagrangian forces and recover a classical Poisson system.The fully discretized algebraic system can be cast into: The velocity u n+1 and the Lagrangian forces f are coupled together to form a new unknown γ n+1 , as follows: ) where Two successive block-LU decompositions are applied to decouple the Lagrangian forces f from γ n+1 and to decouple the velocity from the pressure field.
The first block-LU decomposition decouples the pressure field from the new unknown γ n+1 : which leads to the following sequence of operations: A second block-LU decomposition is applied to the first equation above: and we end up with the following sequence: The decoupled version of the IBPM implemented in PetIBM solves, at every time step, Equations ( 15) to (17) followed by Equations ( 12) and ( 13).