David Shirokoff

David Shirokoff, Rodolfo Ruben Rosales

Common efficient schemes for the incompressible Navier-Stokes equations, such as projection or fractional step methods, have limited temporal accuracy as a result of matrix splitting errors, or introduce errors near the domain boundaries (which destroy uniform convergence to the solution). In this paper we recast the incompressible (constant density) Navier-Stokes equations (with the velocity prescribed at the boundary) as an equivalent system, for the primary variables velocity and pressure. We do this in the usual way away from the boundaries, by replacing the incompressibility condition on the velocity by a Poisson equation for the pressure. The key difference from the usual approaches occurs at the boundaries, where we use boundary conditions that unequivocally allow the pressure to be recovered from knowledge of the velocity at any fixed time. This avoids the common difficulty of an, apparently, over-determined Poisson problem. Since in this alternative formulation the pressure can be accurately and efficiently recovered from the velocity, the recast equations are ideal for numerical marching methods. The new system can be discretized using a variety of methods, in principle to any desired order of accuracy. In this work we illustrate the approach with a 2-D second order finite difference scheme on a Cartesian grid, and devise an algorithm to solve the equations on domains with curved (non-conforming) boundaries, including a case with a non-trivial topology (a circular obstruction inside the domain). This algorithm achieves second order accuracy (in L-infinity), for both the velocity and the pressure. The scheme has a natural extension to 3-D.

David Shirokoff, Jean-Christophe Nave

The volume penalty method provides a simple, efficient approach for solving the incompressible Navier-Stokes equations in domains with boundaries or in the presence of moving objects. Despite the simplicity, the method is typically limited to first order spatial accuracy. We demonstrate that one may achieve high order accuracy by introducing an active penalty term. One key difference from other works is that we use a sharp, unregularized mask function. We discuss how to construct the active penalty term, and provide numerical examples, in dimensions one and two. We demonstrate second and third order convergence for the heat equation, and second order convergence for the Navier-Stokes equations. In addition, we show that modifying the penalty term does not significantly alter the time step restriction from that of the conventional penalty method.

Steven B. Roberts, David Shirokoff, Abhijit Biswas et al.

Classical convergence theory of Runge-Kutta methods assumes that the time step is small relative to the Lipschitz constant of the ordinary differential equation (ODE). For stiff problems, that assumption is often violated, and a problematic degradation in accuracy, known as order reduction, can arise. Methods with high stage order, e.g., Gauss-Legendre and Radau, are known to avoid order reduction, but they must be fully implicit. For the broad class of semilinear ODEs, which consist of a stiff linear term and non-stiff nonlinear term, we show that weaker conditions suffice. Our new semilinear order conditions are formulated in terms of orthogonality relations and can be enumerated by rooted trees. Finally, we prove global error bounds that hold uniformly with respect to stiffness of the linear term.

Ryan Galagusz, David Shirokoff, Jean-Christophe Nave

We present a high order, Fourier penalty method for the Maxwell's equations in the vicinity of perfect electric conductor boundary conditions. The approach relies on extending the smooth non-periodic domain of the equations to a periodic domain by removing the exact boundary conditions and introducing an analytic forcing term in the extended domain. The forcing, or penalty term is chosen to systematically enforce the boundary conditions to high order in the penalty parameter, which then allows for higher order numerical methods. We present an efficient numerical method for constructing the penalty term, and discretize the resulting equations using a Fourier spectral method. We demonstrate convergence orders of up to 3.5 for the one-dimensional Maxwell's equations, and show that the numerical method does not suffer from dispersion (or pollution) errors. We also illustrate the approach in two dimensions and demonstrate convergence orders of 2.5 for transverse magnetic modes and 1.5 for the transverse electric modes. We conclude the paper with numerous test cases in dimensions two and three including waves traveling in an irregular waveguide, and scattering off of a windmill-like geometry.

Mahdi Bandegi, David Shirokoff

We present a new approach for computing approximate global minimizers to a large class of non-local pairwise interaction problems defined over probability distributions. The approach predicts candidate global minimizers, with a recovery guarantee, that are sometimes exact, and often within a few percent of the optimum energy (under appropriate normalization of the energy). The procedure relies on a convex relaxation of the pairwise energy that exploits translational symmetry, followed by a recovery procedure that minimizes a relative entropy. Numerical discretizations of the convex relaxation yield a linear programming problem over convex cones that can be solved using well-known methods. One advantage of the approach is that it provides sufficient conditions for global minimizers to a non-convex quadratic variational problem, in the form of a linear, convex, optimization problem for the auto-correlation of the probability density. We demonstrate the approach in a periodic domain for examples arising from models in materials, social phenomena and flocking. The approach also exactly recovers the global minimizer when a lattice of Dirac masses solves the convex relaxation. An important by-product of the relaxation is a decomposition of the pairwise energy functional into the sum of a convex functional and non-convex functional. We observe that in some cases, the non-convex component of the decomposition can be used to characterize the support of the recovered minimizers.

Rodolfo Ruben Rosales, Benjamin Seibold, David Shirokoff et al.

This paper presents a new class of high order linear ImEx multistep schemes with large regions of unconditional stability. Unconditional stability is a desirable property of a time stepping scheme, as it allows the choice of time step solely based on accuracy considerations. Of particular interest are problems for which both the implicit and explicit parts of the ImEx splitting are stiff. Such splittings can arise, for example, in variable-coefficient problems, or the incompressible Navier-Stokes equations. To characterize the new ImEx schemes, an unconditional stability region is introduced, which plays a role analogous to that of the stability region in conventional multistep methods. Moreover, computable quantities (such as a numerical range) are provided that guarantee an unconditionally stable scheme for a proposed implicit-explicit matrix splitting. The new approach is illustrated with several examples. Coefficients of the new schemes up to fifth order are provided.

Benjamin Seibold, David Shirokoff, Dong Zhou

This paper focuses on the question of how unconditional stability can be achieved via multistep ImEx schemes, in practice problems where both the implicit and explicit terms are allowed to be stiff. For a class of new ImEx multistep schemes that involve a free parameter, strategies are presented on how to choose the ImEx splitting and the time stepping parameter, so that unconditional stability is achieved under the smallest approximation errors. These strategies are based on recently developed stability concepts, which also provide novel insights into the limitations of existing semi-implicit backward differentiation formulas (SBDF). For instance, the new strategies enable higher order time stepping that is not otherwise possible with SBDF. With specific applications in nonlinear diffusion problems and incompressible channel flows, it is demonstrated how the unconditional stability property can be leveraged to efficiently solve stiff nonlinear or nonlocal problems without the need to solve nonlinear or nonlocal problems implicitly.