Dong Zhou

Department of Mathematics
Temple University

In many computational science and engineering applications, for example, fluid flow and interface propagation, it is extremely challenging to design numerical methods that achieve high order of accuracy. The main theme of my research is the development of numerical schemes that are high-order both in space and time, for partial differential equations arise in these applications.
A fundamental drawback of high-order linear multistep methods is their restrictive stability properties compared to RK schemes. In this work, we develop a new stability theory that characterizes unconditional stability of IMEX multistep methods, and design a new class of IMEX multistep methods that possesses large regions of unconditional stability. Specifically, we consider IMEX multistep methods applied to a system of ordinary differential equations $u_t = Au + Bu$, where $Au$ is treated implicitly and $Bu$ is treated explicitly. We are particularly interested in the problem for which both $A$ and $B$ are stiff. Such splittings can arise, for example, in a pressure Poisson reformulations of the incompressible Navier-Stokes equations. Unconditional stability allows for time step selection solely based on accuracy consideration.

The new IMEX multistep methods have tunable parameters that control the size of the unconditional stability region, and can achieve unconditional stability for problems that are unstable with current popular schemes such as semi-implicit backward differentiation formulas (SBDF). Our new stability theory generalizes the existing stability theories, that make use of stability regions and matrix spectra of the problem, to IMEX multistep methods, and provides rigorous explanations for stability limitations in SBDF.

The figures above illustrate the unconditional stability regions for our new IMEX multistep methods from order 2 to 5, with different choice of the parameter $\delta$. Note that $\delta=1$ recovers SBDF methods.

Related Publications:
Despite the advantages of locality in time and good stability properties, RK methods have a key drawback: in general they are affected by order reduction, meaning that the observed convergence order, when applied to stiff ODEs or initial boundary value problems (IBVP), is less than the formal order of the scheme.

In this work, we analyze the global errors incurred when applying an implicit RK scheme to a general IBVP $u_t = \mathcal{L}u + f(x, t)$ with time dependent boundary conditions. We develop a theoretical framework that yields a fundamental geometric understanding of the mechanism for order reduction. Geometrically, the order reduction stems from the boundary layers, produced by the fact that the scheme is too accurate near the boundary (the approximation error satisfies a singularly perturbed boundary value problem). The amplitude of the boundary layers is limited by the stage order of the scheme, which in general less than the scheme’s order, thus order reduction.

Two remedies are developed to recover full convergence order:
(1) We derive a new condition, the weak stage order (WSO), for constructing new RK schemes devoid of order reduction. The novelty of this condition is that it is compatible with a diagonally implicit Runge-Kutta (DIRK) structure, in contrast to the high stage order which is impossible to achieve with DIRK schemes.
(2) The modified boundary conditions (MBC) achieve order reduction avoidance by suppressing the singular behavior in the numerical solutions computed from existing RK schemes. These boundary conditions are extracted via a formal power series expansion of the solution in the absence of the the singular behavior triggered by the conventional approach of enforcing boundary conditions.

The figures above show the error convergence for 3rd order schemes for a 1D heat equation example. Left: 3rd order WSO 1 DIRK scheme with conventional b.c. exhibits order reduction. Middle: The same WSO 1 DIRK scheme with MBC recovers full 3rd order convergence in $u$ and $u_x$. Right: 3rd order WSO 2 DIRK scheme recovers 3rd order convergence in $u$.

Related Publications:
Pressure Poisson equation (PPE) reformulations of the incompressible Navier-Stokes equations is a class of methods that replace the incompressibility constraint by a Poisson equation for the pressure, with a proper choice of the boundary condition so that the incompressibility is maintained. This approach has important advantages compared to the standard form: (1) the pressure is not implicitly coupled to the velocity and can be recovered by solving a Poisson equation; (2) no numerical boundary layers are generated, in contrast to the projection type methods.

In this work, I proposed numerical methods for solving one PPE reformulation (Shirokoff and Rosales, 2011) using mixed finite element methods and meshfree finite differences in space with high-order implicit-explicit (IMEX) time-stepping in time. This PPE reformulation is special due to its electric boundary conditions ($\hat{n}\times\vec{u}=0$ and $\nabla\cdot\vec{u}=0$). The standard nodal-based FE approach will converge to the wrong solution on domains with reentrant corners or curved boundaries. The former is due to the insufficient regularity of the solution, and the latter is a form of Babuska paradox.

The mixed finite element approach is well-suited for electric boundary conditions, and the meshfree finite difference handles curved boundaries naturally. The IMEX Runge-Kutta methods are used to treat the viscous term implicitly and the pressure explicitly. The method achieves 3rd order convergence in time with a 3rd order IMEX RK method, for problems with time-independent boundary conditions and forcing at the boundary (a special case where order reduction does not occur).

Some Benchmark Test Results (Lid-driven Cavity, Backward-facing Step)

Related Publications:
• D. Zhou, B. Seibold, D. Shirokoff, P. Chidyagwai and R. R. Rosales. Meshfree finite differences for vector Poisson and pressure Poisson equations with electric boundary conditions. Meshfree methods for Partial Differential Equations VII, Lecture Notes in Computational Science and Engineering, Vol. 100. Griebel, M. and Schweitzer, M.A. (ed.), Springer, 2015.
• D. Zhou, B. Seibold, D. Shirokoff, and R. R. Rosales. High-order mixed finite element methods for a pressure Poisson equation reformulation of the Navier-Stokes equations with electric boundary conditions. In preparation.
Ph.D. Thesis:
• D. Zhou. High-order numerical methods for pressure Poisson equation reformulations of the incompressible Navier-Stokes equations. Ph.D. Thesis, Temple University, 2014.
Jet schemes are based on evolving parts of the jet of the solution in time, and solving linear advection equations with high order of accuracy by tracking characteristic curves and using suitable Hermite interpolants. Click here for more on Jet schemes for linear advection problems.

My work is focusing on extending this scheme to nonlinear Hamilton-Jacobi type equations, where the characteristics may collide or emenate from a single point. Solutions to such problems can have kinks even with smooth initial conditions. Tracking characteristics gives rise to minimization problems at each time step. The overall scheme based on bi-cubic Hermite interpolation, approximates the entropy solution correctly and achieves optimal 3rd order accuracy in the smooth parts of the solution, while the low order approximations at the non-smooth parts remain localized.

The animation on the right is a numerical simulation of front propagation with an obstacle in the computational domain. The front is represented implicitly by a levelset function.