Krylov deferred correction methods for differential equations with algebraic constraints

In this dissertation, we introduce a new class of spectral time stepping methods for efficient and accurate solutions of ordinary differential equations (ODEs), differential algebraic equations (DAEs), and partial differential equations (PDEs). The methods are based on applying spectral deferred correction techniques as preconditioners to Picard integral collocation formulations, least squares based orthogonal polynomial approximations are computed using Gaussian type quadratures, and spectral integration is used instead of numerically unstable differentiation. For ODE problems, the resulting Krylov deferred correction (KDC) methods solve the preconditioned nonlinear system using Newton-Krylov schemes such as Newton-GMRES method. For PDE systems, method of lines transpose (MoLT ) couples the KDC techniques with fast elliptic equation solvers based on integral equation formulations and fast algorithms. Preliminary numerical results show that the new methods are of arbitrary order of accuracy, extremely stable, and very competitive with existing techniques, particularly when high precision is desired

Method of lines transpose: High order L-stable O(N) schemes for parabolic equations using successive convolution

We present a new solver for nonlinear parabolic problems that is L-stable and achieves high order accuracy in space and time. The solver is built by first constructing a single-dimensional heat equation solver that uses fast O(N) convolution. This fundamental solver has arbitrary order of accuracy in space, and is based on the use of the Green's function to invert a modified Helmholtz equation. Higher orders of accuracy in time are then constructed through a novel technique known as successive convolution (or resolvent expansions). These resolvent expansions facilitate our proofs of stability and convergence, and permit us to construct schemes that have provable stiff decay. The multi-dimensional solver is built by repeated application of dimensionally split independent fundamental solvers. Finally, we solve nonlinear parabolic problems by using the integrating factor method, where we apply the basic scheme to invert linear terms (that look like a heat equation), and make use of Hermite-Birkhoff interpolants to integrate the remaining nonlinear terms. Our solver is applied to several linear and nonlinear equations including heat, Allen-Cahn, and the Fitzhugh-Nagumo system of equations in one and two dimensions

Method of Lines Transpose: Energy Gradient Flows Using Direct Operator Inversion for Phase-Field Models

In this work, we develop an $\mathcal{O}(N)$ implicit real space method in 1D and 2D for the Cahn--Hilliard (CH) and vector Cahn--Hilliard (VCH) equations, based on the method of lines transpose (MOL$^{T}$) formulation. This formulation results in a semidiscrete time stepping algorithm, which we prove is gradient stable in the $H^{-1}$ norm. The spatial discretization follows from dimensional splitting and an $\mathcal{O}(N)$ matrix-free solver, which applies fast convolution to the modified Helmholtz equation. We propose a novel factorization technique, in which fourth-order spatial derivatives are incorporated into the solver. The splitting error is included in the nonlinear fixed point iteration, resulting in a high-order, logically Cartesian (line-by-line) update. Our method is fast but not restricted to periodic boundaries like the fast Fourier transform (FFT). The basic solver is implemented using the backward Euler formulation, and we extend this to both backward difference formula (BDF) stencils, singly diagonal implicit Runge--Kutta (SDIRK), and spectral deferred correction (SDC) frameworks to achieve high orders of temporal accuracy. We demonstrate with numerical results that the CH and VCH equations maintain gradient stability in one and two spatial dimensions. We also explore time-adaptivity, so that meta-stable states and ripening events can be simulated both quickly and efficiently

Method of Lines Transpose: Energy Gradient Flows Using Direct Operator Inversion for Phase-Field Models

In this work, we develop an $\mathcal{O}(N)$ implicit real space method in 1D and 2D for the Cahn--Hilliard (CH) and vector Cahn--Hilliard (VCH) equations, based on the method of lines transpose (MOL$^{T}$) formulation. This formulation results in a semidiscrete time stepping algorithm, which we prove is gradient stable in the $H^{-1}$ norm. The spatial discretization follows from dimensional splitting and an $\mathcal{O}(N)$ matrix-free solver, which applies fast convolution to the modified Helmholtz equation. We propose a novel factorization technique, in which fourth-order spatial derivatives are incorporated into the solver. The splitting error is included in the nonlinear fixed point iteration, resulting in a high-order, logically Cartesian (line-by-line) update. Our method is fast but not restricted to periodic boundaries like the fast Fourier transform (FFT). The basic solver is implemented using the backward Euler formulation, and we extend this to both backward difference formula (BDF) stencils, singly diagonal implicit Runge--Kutta (SDIRK), and spectral deferred correction (SDC) frameworks to achieve high orders of temporal accuracy. We demonstrate with numerical results that the CH and VCH equations maintain gradient stability in one and two spatial dimensions. We also explore time-adaptivity, so that meta-stable states and ripening events can be simulated both quickly and efficiently

An evaluation of solution algorithms and numerical approximation methods for modeling an ion exchange process

The focus of this work is on the modeling of an ion exchange process that occurs in drinking water treatment applications. The model formulation consists of a two-scale model in which a set of microscale diffusion equations representing ion exchange resin particles that vary in size and age are coupled through a boundary condition with a macroscopic ordinary differential equation (ODE), which represents the concentration of a species in a well-mixed reactor. We introduce a new age-averaged model (AAM) that averages all ion exchange particle ages for a given size particle to avoid the expensive Monte-Carlo simulation associated with previous modeling applications. We discuss two different numerical schemes to approximate both the original Monte Carlo algorithm and the new AAM for this two-scale problem. The first scheme is based on the finite element formulation in space coupled with an existing backward-difference-formula-based ODE solver in time. The second scheme uses an integral equation based Krylov deferred correction (KDC) method and a fast elliptic solver (FES) for the resulting elliptic equations. Numerical results are presented to validate the new AAM algorithm, which is also shown to be more computationally efficient than the original Monte Carlo algorithm. We also demonstrate that the higher order KDC scheme is more efficient than the traditional finite element solution approach and this advantage becomes increasingly important as the desired accuracy of the solution increases. We also discuss issues of smoothness, which affect the efficiency of the KDC-FES approach, and outline additional algorithmic changes that would further improve the efficiency of these developing methods for a wide range of applications

Method of Lines Transpose: High Order L-Stable {O}(N) Schemes for Parabolic Equations Using Successive Convolution

We present a new solver for nonlinear parabolic problems that is L-stable and achieves high order accuracy in space and time. The solver is built by first constructing a one-dimensional heat equation solver that uses fast $\mathcal O(N)$ convolution. This fundamental solver has arbitrary order of accuracy in space and is based on the use of the Green\u27s function to invert a modified Helmholtz equation. Higher orders of accuracy in time are then constructed through a novel technique known as successive convolution (or resolvent expansions). These resolvent expansions facilitate our proofs of stability and convergence, and permit us to construct schemes that have provable stiff decay. The multidimensional solver is built by repeated application of dimensionally split independent fundamental solvers. Finally, we solve nonlinear parabolic problems by using the integrating factor method, where we apply the basic scheme to invert linear terms (that look like a heat equation), and make use of Hermite--Birkhoff interpolants to integrate the remaining nonlinear terms. Our solver is applied to several linear and nonlinear equations including heat, Allen--Cahn, and the FitzHugh--Nagumo system of equations in one and two dimensions