Source author record

Jon Wilkening

Jon Wilkening appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

19works
17topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

19 published item(s)

preprint2022arXiv

Numerical Algorithms for Water Waves with Background Flow over Obstacles and Topography

We present two accurate and efficient algorithms for solving the incompressible, irrotational Euler equations with a free surface in two dimensions with background flow over a periodic, multiply-connected fluid domain that includes stationary obstacles and variable bottom topography. One approach is formulated in terms of the surface velocity potential while the other evolves the vortex sheet strength. Both methods employ layer potentials in the form of periodized Cauchy integrals to compute the normal velocity of the free surface, are compatible with arbitrary parameterizations of the free surface and boundaries, and allow for circulation around each obstacle, which leads to multiple-valued velocity potentials but single-valued stream functions. We prove that the resulting second-kind Fredholm integral equations are invertible, possibly after a physically motivated finite-rank correction. In an angle-arclength setting, we show how to avoid curve reconstruction errors that are incompatible with spatial periodicity. We use the proposed methods to study gravity-capillary waves generated by flow around several elliptical obstacles above a flat or variable bottom boundary. In each case, the free surface eventually self-intersects in a splash singularity or collides with a boundary. We also show how to evaluate the velocity and pressure with spectral accuracy throughout the fluid, including near the free surface and solid boundaries. To assess the accuracy of the time evolution, we monitor energy conservation and the decay of Fourier modes and compare the numerical results of the two methods to each other. We implement several solvers for the discretized linear systems and compare their performance. The fastest approach employs a graphics processing unit (GPU) to construct the matrices and carry out iterations of the generalized minimal residual method (GMRES).

preprint2020arXiv

Harmonic Stability of Standing Water Waves

A numerical method is developed to study the stability of standing water waves and other time-periodic solutions of the free-surface Euler equations using Floquet theory. A Fourier truncation of the monodromy operator is computed by solving the linearized Euler equations about the standing wave with initial conditions ranging over all Fourier modes up to a given wave number. The eigenvalues of the truncated monodromy operator are computed and ordered by the mean wave number of the corresponding eigenfunctions, which we introduce as a method of retaining only accurately computed Floquet multipliers. The mean wave number matches up with analytical results for the zero-amplitude standing wave and is helpful in identifying which Floquet multipliers collide and leave the unit circle to form unstable eigenmodes or rejoin the unit circle to regain stability. For standing waves in deep water, most waves with crest acceleration below $A_c=0.889$ are found to be linearly stable to harmonic perturbations; however, we find several bubbles of instability at lower values of $A_c$ that have not been reported previously in the literature. We also study the stability of several new or recently discovered time-periodic gravity-capillary or gravity waves in deep or shallow water, finding several examples of large-amplitude waves that are stable to harmonic perturbations and others that are not. A new method of matching the Floquet multipliers of two nearby standing waves by solving a linear assignment problem is also proposed to track individual eigenvalues via homotopy from the zero-amplitude state to large-amplitude standing waves.

preprint2016arXiv

A Fully Discrete Adjoint Method for Optimization of Flow Problems on Deforming Domains with Time-Periodicity Constraints

A variety of shooting methods for computing fully discrete time-periodic solutions of partial differential equations, including Newton-Krylov and optimization-based methods, are discussed and used to determine the periodic, compressible, viscous flow around a 2D flapping airfoil. The Newton-Krylov method uses matrix-free GMRES to solve the linear systems of equations that arise in the nonlinear iterations, with matrix-vector products computed via the linearized sensitivity evolution equations. The adjoint method is used to compute gradients for the gradient-based optimization shooting methods. The Newton-Krylov method is shown to exhibit superior convergence to the optimal solution for these fluid problems, and fully leverages quality starting data. The central contribution of this work is the derivation of the adjoint equations and the corresponding adjoint method for fully discrete, time-periodically constrained partial differential equations. These adjoint equations constitute a linear, two-point boundary value problem that is provably solvable. The periodic adjoint method is used to compute gradients of quantities of interest along the manifold of time-periodic solutions of the discrete partial differential equation, which is verified against a second-order finite difference approximation. These gradients are then used in a gradient-based optimization framework to determine the energetically optimal flapping motion of a 2D airfoil in compressible, viscous flow over a single cycle, such that the time-averaged thrust is identically zero. In less than 20 optimization iterations, the flapping energy was reduced nearly an order of magnitude and the thrust constraint satisfied to 5 digits of accuracy.

preprint2016arXiv

Eigenfunctions and the Dirichlet problem for the Classical Kimura Diffusion Operator

We study the classical Kimura diffusion operator defined on the n-simplex, $$L^{Kim}=\sum_{1\leq i,j\leq n+1}x_ix_j\partial_{x_i}\partial_{x_j}$$ We give novel constructions for the basis of eigenpolynomials, and the solution to the inhomogeneous Dirichlet problem, which are well adapted to numerical applications. Our solution of the Dirichlet problem is quite explicit and provides a precise description of the singularities that arise along the boundary.

preprint2016arXiv

Optimizing intermittent water supply in urban pipe distribution networks

In many urban areas of the developing world, piped water is supplied only intermittently, as valves direct water to different parts of the water distribution system at different times. The flow is transient, and may transition between free-surface and pressurized, resulting in complex dynamical features with important consequences for water suppliers and users. Here, we develop a computational model of transition, transient pipe flow in a network, accounting for a wide variety of realistic boundary conditions. We validate the model against several published data sets, and demonstrate its use on a real pipe network. The model is extended to consider several optimization problems motivated by realistic scenarios. We demonstrate how to infer water flow in a small pipe network from a single pressure sensor, and show how to control water inflow to minimize damaging pressure gradients.

preprint2016arXiv

The instability of Wilton ripples

Wilton ripples are a type of periodic traveling wave solution of the full water wave problem incorporating the effects of surface tension. They are characterized by a resonance phenomenon that alters the order at which the resonant harmonic mode enters in a perturbation expansion. We compute such solutions using non-perturbative numerical methods and investigate their stability by examining the spectrum of the water wave problem linearized about the resonant traveling wave. Instabilities are observed that differ from any previously found in the context of the water wave problem.

preprint2015arXiv

A Spectral Transform Method for Singular Sturm-Liouville Problems with Applications to Energy Diffusion in Plasma Physics

We develop a spectrally accurate numerical method to compute solutions of a model partial differential equation used in plasma physics to describe diffusion in velocity space due to Fokker-Planck collisions. The solution is represented as a discrete and continuous superposition of normalizable and non-normalizable eigenfunctions via the spectral transform associated with a singular Sturm-Liouville operator. We present a new algorithm for computing the spectral density function of the operator that uses Chebyshev polynomials to extrapolate the value of the Titchmarsh-Weyl $m$-function from the complex upper half-plane to the real axis. The eigenfunctions and density function are rescaled and a new formula for the limiting value of the $m$-function is derived to avoid amplification of roundoff errors when the solution is reconstructed. The complexity of the algorithm is also analyzed, showing that the cost of computing the spectral density function at a point grows less rapidly than any fractional inverse power of the desired accuracy. A WKB analysis is used to prove that the spectral density function is real analytic. Using this new algorithm, we highlight key properties of the partial differential equation and its solution that have strong implications on the optimal choice of discretization method in large-scale plasma physics computations.

preprint2015arXiv

Accurate Spectral Numerical Schemes for Kinetic Equations with Energy Diffusion

We examine the merits of using a family of polynomials that are orthogonal with respect to a non-classical weight function to discretize the speed variable in continuum kinetic calculations. We consider a model one-dimensional partial differential equation describing energy diffusion in velocity space due to Fokker-Planck collisions. This relatively simple case allows us to compare the results of the projected dynamics with an expensive but highly accurate spectral transform approach. It also allows us to integrate in time exactly, and to focus entirely on the effectiveness of the discretization of the speed variable. We show that for a fixed number of modes or grid points, the non-classical polynomials can be many orders of magnitude more accurate than classical Hermite polynomials or finite-difference solvers for kinetic equations in plasma physics. We provide a detailed analysis of the difference in behavior and accuracy of the two families of polynomials. For the non-classical polynomials, if the initial condition is not smooth at the origin when interpreted as a three-dimensional radial function, the exact solution leaves the polynomial subspace for a time, but returns (up to roundoff accuracy) to the same point evolved to by the projected dynamics in that time. By contrast, using classical polynomials, the exact solution differs significantly from the projected dynamics solution when it returns to the subspace. We also explore the connection between eigenfunctions of the projected evolution operator and (non-normalizable) eigenfunctions of the full evolution operator, as well as the effect of truncating the computational domain.

preprint2015arXiv

Parameter estimation by implicit sampling

Implicit sampling is a weighted sampling method that is used in data assimilation, where one sequentially updates estimates of the state of a stochastic model based on a stream of noisy or incomplete data. Here we describe how to use implicit sampling in parameter estimation problems, where the goal is to find parameters of a numerical model, e.g.~a partial differential equation (PDE), such that the output of the numerical model is compatible with (noisy) data. We use the Bayesian approach to parameter estimation, in which a posterior probability density describes the probability of the parameter conditioned on data and compute an empirical estimate of this posterior with implicit sampling. Our approach generates independent samples, so that some of the practical difficulties one encounters with Markov Chain Monte Carlo methods, e.g.~burn-in time or correlations among dependent samples, are avoided. We describe a new implementation of implicit sampling for parameter estimation problems that makes use of multiple grids (coarse to fine) and BFGS optimization coupled to adjoint equations for the required gradient calculations. The implementation is "dimension independent", in the sense that a well-defined finite dimensional subspace is sampled as the mesh used for discretization of the PDE is refined. We illustrate the algorithm with an example where we estimate a diffusion coefficient in an elliptic equation from sparse and noisy pressure measurements. In the example, dimension\slash mesh-independence is achieved via Karhunen-Loève expansions.

preprint2014arXiv

Boltzmann Equation Solver Adapted to Emergent Chemical Non-equilibrium

We present a novel method to solve the spatially homogeneous and isotropic relativistic Boltzmann equation. We employ a basis set of orthogonal polynomials dynamically adapted to allow for emergence of chemical non-equilibrium. Two time dependent parameters characterize the set of orthogonal polynomials, the effective temperature $T(t)$ and phase space occupation factor $Υ(t)$. In this first paper we address (effectively) massless fermions and derive dynamical equations for $T(t)$ and $Υ(t)$ such that the zeroth order term of the basis alone captures the particle number density and energy density of each particle distribution. We validate our method and illustrate the reduced computational cost and the ability to easily represent final state chemical non-equilibrium by studying a model problem that is motivated by the physics of the neutrino freeze-out processes in the early Universe, where the essential physical characteristics include reheating from another disappearing particle component ($e^\pm$-annihilation).

preprint2014arXiv

Comparison of five methods of computing the Dirichlet-Neumann operator for the water wave problem

We compare the effectiveness of solving Dirichlet-Neumann problems via the Craig-Sulem (CS) expansion, the Ablowitz-Fokas-Musslimani (AFM) implicit formulation, the dual AFM formulation (AFM*), a boundary integral collocation method (BIM), and the transformed field expansion (TFE) method. The first three methods involve highly ill-conditioned intermediate calculations that we show can be overcome using multiple-precision arithmetic. The latter two methods avoid catastrophic cancellation of digits in intermediate results, and are much better suited to numerical computation. For the Craig-Sulem expansion, we explore the cancellation of terms at each order (up to 150th) for three types of wave profiles, namely band-limited, real-analytic, or smooth. For the AFM and AFM* methods, we present an example in which representing the Dirichlet or Neumann data as a series using the AFM basis functions is impossible, causing the methods to fail. The example involves band-limited wave profiles of arbitrarily small amplitude, with analytic Dirichlet data. We then show how to regularize the AFM and AFM* methods by over-sampling the basis functions and using the singular value decomposition or QR-factorization to orthogonalize them. Two additional examples are used to compare all five methods in the context of water waves, namely a large-amplitude standing wave in deep water, and a pair of interacting traveling waves in finite depth.

preprint2014arXiv

Relative-Periodic Elastic Collisions of Water Waves

We compute time-periodic and relative-periodic solutions of the free-surface Euler equations that take the form of overtaking collisions of unidirectional solitary waves of different amplitude on a periodic domain. As a starting guess, we superpose two Stokes waves offset by half the spatial period. Using an overdetermined shooting method, the background radiation generated by collisions of the Stokes waves is tuned to be identical before and after each collision. In some cases, the radiation is effectively eliminated in this procedure, yielding smooth soliton-like solutions that interact elastically forever. We find examples in which the larger wave subsumes the smaller wave each time they collide, and others in which the trailing wave bumps into the leading wave, transferring energy without fully merging. Similarities notwithstanding, these solutions are found quantitatively to lie outside of the Korteweg-de Vries regime. We conclude that quasi-periodic elastic collisions are not unique to integrable model water wave equations when the domain is periodic.

preprint2012arXiv

Overdetermined Shooting Methods for Computing Standing Water Waves with Spectral Accuracy

A high-performance shooting algorithm is developed to compute time-periodic solutions of the free-surface Euler equations with spectral accuracy in double and quadruple precision. The method is used to study resonance and its effect on standing water waves. We identify new nucleation mechanisms in which isolated large-amplitude solutions, and closed loops of such solutions, suddenly exist for depths below a critical threshold. We also study degenerate and secondary bifurcations related to Wilton's ripples in the traveling case, and explore the breakdown of self-similarity at the crests of extreme standing waves. In shallow water, we find that standing waves take the form of counter-propagating solitary waves that repeatedly collide quasi-elastically. In deep water with surface tension, we find that standing waves resemble counter-propagating depression waves. We also discuss existence and non-uniqueness of solutions, and smooth versus erratic dependence of Fourier modes on wave amplitude and fluid depth. In the numerical method, robustness is achieved by posing the problem as an overdetermined nonlinear system and using either adjoint-based minimization techniques or a quadratically convergent trust-region method to minimize the objective function. Accuracy is maintained using spectral collocation with optional mesh refinement in space, a high order Runge-Kutta or spectral deferred correction method in time, and quadruple-precision for improved navigation of delicate regions of parameter space as well as validation of double-precision results. Implementation issues for GPU acceleration are briefly discussed, and the performance of the algorithm is tested for a number of hardware configurations.

preprint2011arXiv

Breakdown of self-similarity at the crests of large amplitude standing water waves

We study the limiting behavior of large-amplitude standing waves on deep water using high-resolution numerical simulations in double and quadruple precision. While periodic traveling waves approach Stokes's sharply crested extreme wave in an asymptotically self-similar manner, we find that standing waves behave differently. Instead of sharpening to a corner or cusp as previously conjectured, the crest tip develops a variety of oscillatory structures. This causes the bifurcation curve that parametrizes these waves to fragment into disjoint branches corresponding to the different oscillation patterns that occur. In many cases, a vertical jet of fluid pushes these structures upward, leading to wave profiles commonly seen in wave tank experiments. Thus, we observe a rich array of dynamic behavior at small length scales in a regime previously thought to be self-similar.

preprint2010arXiv

A local construction of the Smith normal form of a matrix polynomial

We present an algorithm for computing a Smith form with multipliers of a regular matrix polynomial over a field. This algorithm differs from previous ones in that it computes a local Smith form for each irreducible factor in the determinant separately and then combines them into a global Smith form, whereas other algorithms apply a sequence of unimodular row and column operations to the original matrix. The performance of the algorithm in exact arithmetic is reported for several test cases.

preprint2010arXiv

Computation of Time-Periodic Solutions of the Benjamin-Ono Equation

We present a spectrally accurate numerical method for finding non-trivial time-periodic solutions of non-linear partial differential equations. The method is based on minimizing a functional (of the initial condition and the period) that is positive unless the solution is periodic, in which case it is zero. We solve an adjoint PDE to compute the gradient of this functional with respect to the initial condition. We include additional terms in the functional to specify the free parameters, which, in the case of the Benjamin-Ono equation, are the mean, a spatial phase, a temporal phase and the real part of one of the Fourier modes at $t=0$. We use our method to study global paths of non-trivial time-periodic solutions connecting stationary and traveling waves of the Benjamin-Ono equation. As a starting guess for each path, we compute periodic solutions of the linearized problem by solving an infinite dimensional eigenvalue problem in closed form. We then use our numerical method to continue these solutions beyond the realm of linear theory until another traveling wave is reached. By experimentation with data fitting, we identify the analytical form of the solutions on the path connecting the one-hump stationary solution to the two-hump traveling wave. We then derive exact formulas for these solutions by explicitly solving the system of ODE's governing the evolution of solitons using the ansatz suggested by the numerical simulations.

preprint2010arXiv

Global paths of time-periodic solutions of the Benjamin-Ono equation connecting pairs of traveling waves

We classify all bifurcations from traveling waves to non-trivial time-periodic solutions of the Benjamin-Ono equation that are predicted by linearization. We use a spectrally accurate numerical continuation method to study several paths of non-trivial solutions beyond the realm of linear theory. These paths are found to either re-connect with a different traveling wave or to blow up. In the latter case, as the bifurcation parameter approaches a critical value, the amplitude of the initial condition grows without bound and the period approaches zero. We then prove a theorem that gives the mapping from one bifurcation to its counterpart on the other side of the path and exhibits exact formulas for the time-periodic solutions on this path. The Fourier coefficients of these solutions are power sums of a finite number of particle positions whose elementary symmetric functions execute simple orbits (circles or epicycles) in the unit disk of the complex plane. We also find examples of interior bifurcations from these paths of already non-trivial solutions, but we do not attempt to analyze their analytic structure.

preprint2010arXiv

Practical Error Estimates for Reynolds' Lubrication Approximation and its Higher Order Corrections

Reynolds' lubrication approximation is used extensively to study flows between moving machine parts, in narrow channels, and in thin films. The solution of Reynolds' equation may be thought of as the zeroth order term in an expansion of the solution of the Stokes equations in powers of the aspect ratio $ε$ of the domain. In this paper, we show how to compute the terms in this expansion to arbitrary order on a two-dimensional, $x$-periodic domain and derive rigorous, a-priori error bounds for the difference between the exact solution and the truncated expansion solution. Unlike previous studies of this sort, the constants in our error bounds are either independent of the function $h(x)$ describing the geometry, or depend on $h$ and its derivatives in an explicit, intuitive way. Specifically, if the expansion is truncated at order $2k$, the error is $O(ε^{2k+2})$ and $h$ enters into the error bound only through its first and third inverse moments $\int_0^1 h(x)^{-m} dx$, $m=1,3$ and via the max norms $\big\|\frac{1}{\ell!} h^{\ell-1} \partial_x^\ell h\big\|_\infty$, $1\le\ell\le2k+2$. We validate our estimates by comparing with finite element solutions and present numerical evidence that suggests that even when $h$ is real analytic and periodic, the expansion solution forms an asymptotic series rather than a convergent series.

preprint2010arXiv

Variational Particle Schemes for the Porous Medium Equation and for the System of Isentropic Euler Equations

Both the porous medium equation and the system of isentropic Euler equations can be considered as steepest descents on suitable manifolds of probability measures in the framework of optimal transport theory. By discretizing these variational characterizations instead of the partial differential equations themselves, we obtain new schemes with remarkable stability properties. We show that they capture successfully the nonlinear features of the flows, such as shocks and rarefaction waves for the isentropic Euler equations. We also show how to design higher order methods for these problems in the optimal transport setting using backward differentiation formula (BDF) multi-step methods or diagonally implicit Runge-Kutta methods.