Magnetohydrodynamics couples the Navier–Stokes and Maxwell equations to describes flows in electrically conducting fluids. The divergence of the magnetic field must vanish, but numerical algorithms typically do not preserve this condition exactly. Artifacts can then arise in solutions, such as spurious forces parallel to the magnetic field. These artifacts can be alleviated by using extended sets of Maxwell’s equations that include magnetic charges and currents and hence are invariant under duality rotations that interchange the electric and magnetic fields. The eight-wave formulation supposes that the magnetic current arises from magnetic charges being advected by the fluid. The magnetohydrodynamic equations are then Galilean invariant even when the divergence of the magnetic field is nonzero. The evolution equation for the magnetic field then resembles Jeffery’s equation that describes the orientations of a suspension of axisymmetric particles. The ideal electric field is invariant under the extra terms proportional to the divergence of the magnetic field. This property leads to particularly simple lattice Boltzmann formulations for two of the three variants, with different treatments of the momentum equation, that are constructed and compared. The formulation that implements the Lorentz force directly, rather than via the Maxwell stress, proves more stable in numerical experiments.

The discrete time quantum walk is a quantum cellular automaton whose wavefunction comprises pairs of complex numbers assigned to uniformly spaced points on a line. The wavefunction evolves through the application of an alternating sequence of unitary operators: streaming of wavefunction values to adjacent points, and a Hadamard-type unitary matrix to blend pairs of values at individual points. Each operator generates the exact evolution due to part of the Hamiltonian for the one-dimensional Dirac equation over a finite time step. Composing these operators thus creates a discrete approximation to the Dirac equation. However, the composition of two non-commuting operators creates a global splitting error proportional to the length of the time step. The global error can be reduced from first order to second order in the time step by a unitary pre-and post-processing of the initial conditions and final output. The algorithm then becomes equivalent to a symmetric composition, a Strang splitting, between the two operators. This paper describes a fourth-order accurate composition scheme using nine stages, the fewest possible when the lengths of the time steps employed in the different stages are constrained to be integer multiples of some base time step. Each stage is itself a symmetric composition between two operators. This fourth-order scheme produces quantitatively smaller errors for a typical benchmark problem on spatial lattices with 1024 or more points, and shows the expected fourth-order convergence on sufficiently fine lattices. It has greater accuracy, over sufficiently long times, than three better-known fourth-order composition schemes using fewer stages, but with lengths related by irrational coefficients. The truncation error for plane-wave solutions is due to an operator that separates into a resonant part proportional to the Hamiltonian, and a non-resonant part orthogonal to the Hamiltonian. The resonant part commutes with the exact evolution operator, so its error accumulates to grow linearly with time. The orthogonal part produces oscillations that remain bounded over many time steps. The nine-stage integer scheme has the smallest resonant truncation error of the four schemes, despite being the only scheme that can be implemented using local operations. The other schemes implement streaming by irrational fractions of the lattice spacing using discrete Fourier transforms.

The two-relaxation-time collision operator in discrete kinetic theory models collisions between particles by grouping them into pairs with anti parallel velocities. It prescribes a linear relaxation towards equilibrium with one rate for the even combination of distribution functions for each pair, and another rate for the odd combination. We reformulate this collision operator using relaxation rates for the forward-propagating and backward-propagating combinations instead. An optimal pair of relaxation rates sets the forward propagating combination of each pair of distributions to equilibrium. Only the backward-propagating non-equilibrium distributions remain. Applying this result twice gives closed discrete equations for evolving the macroscopic variables alone across three time levels. We split the equivalent equations into a first order system: a conservation law and a kinetic equation for the flux. All other quantities are evaluated at equilibrium. We apply this formalism to the magnetic field in a lattice Boltzmann scheme for magnetohydrodynamics. The antisymmetric part of the kinetic equation matches the Maxwell-Faraday equation and Ohm's law. The symmetric part matches the hyperbolic divergence cleaning model. The discrete divergence of the magnetic field remains zero, to within round-off error, when the initial magnetic field is the discrete curl of a vector potential. We have thus constructed a mimetic or constrained transport scheme for magnetohydrodynamics.

Magnetohydrodynamics couples the Navier-Stokes and Maxwell's equations to describe the flow of electrically conducting fluids in magnetic fields. Maxwell's equations require the divergence of the magnetic field to vanish, but this condition is typically not preserved exactly by numerical algorithms. Solutions can develop artifacts because structural properties of the magnetohydrodynamic equations then fail to hold. Magnetohydrodynamics with hyperbolic divergence cleaning permits a nonzero divergence that evolves under a telegraph equation, designed to both damp the divergence, and propagate it away from any sources, such as poorly resolved regions with large spatial gradients, without significantly increasing the computational cost. We show that existing lattice Boltzmann algorithms for magnetohydrodynamics already incorporate hyperbolic divergence cleaning, though they typically use parameter values for which it reduces to parabolic divergence cleaning under a slowly-varying approximation. We recover hyperbolic divergence cleaning by adjusting the relaxation rate for the trace of the tensor that represents the electric field, and absorb the contribution from the symmetric-traceless part of this tensor using a change of variables. Numerical experiments confirm that the qualitative behaviour changes from parabolic to hyperbolic when the relaxation time for the trace of the electric field tensor is increased.

We present a lattice Boltzmann algorithm for simulating magnetohydrodynamics, and extend it to simulate the Jeffery equation that describes the rotating orientations of axisymmetric particles in a dilute suspension. Both systems involve material vector fields that evolve through the curl of another vector field. Both systems thus require an underlying kinetic formulation using vector fields, in contrast to the scalar fields used in the Boltzmann equation, and in lattice Boltzmann algorithms for hydrodynamics. Simulating Jeffery’s equation requires extra gradient terms that cannot be written in conservation form. These gradients are obtained locally at grid points using the non-equilibrium parts of the kinetic vector fields representing the particle orientations, and the kinetic scalar fields representing the suspending fluid. The kinetic formulation is discretised using a Strang splitting between advection to neighbouring grid points and local algebraic operations at grid points.

There are a large number of industrial applications which involve the extrusion of aerated composite materials through an orifice. Often such delicate structures change due to shear and pressure fields experienced on extrusion. We formulate and solve a two-phase flow model for the motion of aerated composite materials. We find that the main resistance to the flow was getting the material out through the orifice. In one-dimensional motion, we find that the material merely translates with no change in composition. In a model describing the motion of air bubbles in a slender block of material, we find that the geometry is critical and that the volume fraction of air is not constant throughout the block. Finally, we formulate a model to describe the flow of a one-phase power-law fluid in a cone and we calculate the power required to move the material down the cone.

The Du Fort–Frankel scheme for the one-dimensional Schrödinger equation is shown to be equivalent, under a time-dependent unitary transformation, to the Ablowitz–Kruskal–Ladik scheme for the Klein–Gordon equation. The Schrödinger equation describes a non-relativistic quantum particle, while the Klein–Gordon equation describes a relativistic particle. The conditional convergence of the Du Fort–Frankel scheme to solutions of the Schrödinger equation arises because solutions of the Klein–Gordon equation only approximate solutions of the Schrödinger equation in the non-relativistic limit. The time-dependent unitary transformation is the discrete analog of the transformation that arises from seeking a non-relativistic limit using the interaction picture of quantum mechanics to decompose the Klein–Gordon Hamiltonian into the relativistic rest energy and a remainder. The Ablowitz–Kruskal–Ladik scheme is in turn decomposed into a quantum lattice gas automaton for the one-dimensional Dirac equation, which is also the one-dimensional discrete time quantum walk. This relativistic interpretation clarifies the origin of the known discrete invariant of the Du Fort–Frankel scheme as expressing conservation of probability for the 2-component wavefunction in the one-dimensional Dirac equation under discrete unitary evolution. It also leads to a second invariant, the matrix element of the evolution operator, whose imaginary part gives a discrete approximation to the expectation of the non-relativistic Schrödinger Hamiltonian.

Numerical simulations of the shallow water equations on rotating spheres produce mixtures of robust vortices and alternating zonal jets, as seen in the atmospheres of the gas giant planets. However, simulations that include Rayleigh friction invariably produce a sub-rotating (retrograde) equatorial jet for Jovian parameter regimes, whilst observations of Jupiter show a super-rotating (prograde) equatorial jet that has persisted over several decades. Super-rotating equatorial jets have recently been obtained in shallow water simulations that include a Newtonian relaxation of perturbations to the layer thickness to model radiative cooling to space, and in simulations of the thermal shallow water equations that include a similar relaxation term in their temperature equation. Simulations of global quasigeostrophic forms of these different models produce equatorial jets in the same directions as the parent models, suggesting that the mechanism responsible for setting the direction lies within quasigeostrophic theory. We provide such a mechanism by calculating the effective force acting on the thickness-weighted zonal mean flow due to the decay of an equatorially trapped Rossby wave. Decay due to Newtonian cooling creates an eastward zonal mean flow at the equator, consistent with the formation of a super-rotating equatorial jet, while decay due to Rayleigh friction leads to a westward zonal mean flow at the equator, consistent with the formation of a sub-rotating equatorial jet. In both cases the meridionally integrated zonal mean of the absolute zonal momentum is westward, consistent with the standard result that Rossby waves carry westward pseudomomentum, but this does not preclude the zonal mean flow being eastward on and close to the equator.

Transfer of free energy from large to small velocity-space scales by phase mixing leads to Landau damping in a linear plasma. In a turbulent drift-kinetic plasma, this transfer is statistically nearly canceled by an inverse transfer from small to large velocity-space scales due to “anti-phase-mixing” modes excited by a stochastic form of plasma echo. Fluid moments (density, velocity, and temperature) are thus approximately energetically isolated from the higher moments of the distribution function, so phase mixing is ineffective as a dissipation mechanism when the plasma collisionality is small.

We present an energy- and potential enstrophy-conserving scheme for the non-traditional shallow water equations that include the complete Coriolis force and topography. These integral conservation properties follow from material conservation of potential vorticity in the continuous shallow water equations. The latter property cannot be preserved by a discretisation on a fixed Eulerian grid, but exact conservation of a discrete energy and a discrete potential enstrophy seems to be an effective substitute that prevents any distortion of the forward and inverse cascades in quasi-two dimensional turbulence through spurious sources and sinks of energy and potential enstrophy, and also increases the robustness of the scheme against nonlinear instabilities. We exploit the existing Arakawa–Lamb scheme for the traditional shallow water equations, reformulated by Salmon as a discretisation of the Hamiltonian and Poisson bracket for this system. The non-rotating, traditional, and our non-traditional shallow water equations all share the same continuous Hamiltonian structure and Poisson bracket, provided one distinguishes between the particle velocity and the canonical momentum per unit mass. We have determined a suitable discretisation of the non-traditional canonical momentum, which includes additional coupling between the layer thickness and velocity fields, and modified the discrete kinetic energy to suppress an internal symmetric computational instability that otherwise arises for multiple layers. The resulting scheme exhibits the expected second-order convergence under spatial grid refinement. We also show that the drifts in the discrete total energy and potential enstrophy due to temporal truncation error may be reduced to machine precision under suitable refinement of the timestep using the third-order Adams–Bashforth or fourth-order Runge–Kutta integration schemes.

Quantum lattice algorithms originated with the Feynman checkerboard model for the one-dimensional Dirac equation. They offer discrete models of quantum mechanics in which the complex numbers representing wavefunction values on a discrete spatial lattice evolve through discrete unitary operations. This paper draws together some of the identical, or at least unitarily equivalent, algorithms that have appeared in three largely disconnected strands of research. Treated as conventional numerical algorithms, they are all only first order accurate under refinement of the discrete space/time grid, but may be raised to second order by a unitary change of variables. Much more efficient implementations arise from replacing the evolution through a sequence of unitary intermediate steps with a short path integral formulation that expresses the wavefunction at each spatial point on the most recent time level as a linear combination of values at immediately preceding time levels and neighbouring spatial points. In one dimension, a particularly elegant reformulation replaces two variables at two time levels with a single variable over three time levels. The resulting algorithm is a variational integrator arising from a discrete action principle, and coincides with the Ablowitz–Kruskal–Ladik finite difference scheme for the Klein–Gordon equation.

Coastal currents flowing along continental shelves are a complex dynamical feature of the global ocean. This chapter describes the evolution of a coastal current in terms of large-amplitude shelf waves. It also describes the experimental setup and procedure and characterizes the evolution of large amplitude waves generated by retrograde flow past a headland. The chapter briefly reviews the quasi-geostrophic (QG) equations that underlie the nonlinear wave theory and numerical simulations. It adapts nonlinear shelf wave theory to the annular channel. The chapter describes the numerical scheme for the QG equations. It compares the skill of the theory, numerical solutions, and experiments in predicting the characteristics of large-amplitude wave breaking. Finally, the chapter summarizes the findings and relates the results to previous studies of topographic Rossby waves and coastal currents.

We study Landau damping in the 1+1D Vlasov–Poisson system using a Fourier–Hermite spectral representation. We describe the propagation of free energy in Fourier–Hermite phase space using forwards and backwards propagating Hermite modes recently developed for gyrokinetic theory. We derive a free energy equation that relates the change in the electric field to the net Hermite flux out of the zeroth Hermite mode. In linear Landau damping, decay in the electric field corresponds to forward propagating Hermite modes; in nonlinear damping, the initial decay is followed by a growth phase characterized by the generation of backwards propagating Hermite modes by the nonlinear term. The free energy content of the backwards propagating modes increases exponentially until balancing that of the forward propagating modes. Thereafter there is no systematic net Hermite flux, so the electric field cannot decay and the nonlinearity effectively suppresses Landau damping. These simulations are performed using the fully-spectral 5D gyrokinetics code S pectro GK, modified to solve the 1+1D Vlasov–Poisson system. This captures Landau damping via Hou–Li filtering in velocity space. Therefore the code is applicable even in regimes where phase mixing and filamentation are dominant.

Conventional shallow water theory successfully reproduces many key features of the Jovian atmosphere: a mixture of coherent vortices and stable, large-scale, zonal jets whose amplitude decreases with distance from the equator. However, both freely decaying and forced-dissipative simulations of the shallow water equations in Jovian parameter regimes invariably yield retrograde equatorial jets, while Jupiter itself has a strong prograde equatorial jet. Simulations by Scott and Polvani [“Equatorial superrotation in shallow atmospheres,” Geophys. Res. Lett. 35, L24202 (2008)] have produced prograde equatorial jets through the addition of a model for radiative relaxation in the shallow water height equation. However, their model does not conserve mass or momentum in the active layer, and produces mid-latitude jets much weaker than the equatorial jet. We present the thermal shallow water equations as an alternative model for Jovian atmospheres. These equations permit horizontal variations in the thermodynamic properties of the fluid within the active layer. We incorporate a radiative relaxation term in the separate temperature equation, leaving the mass and momentum conservation equations untouched. Simulations of this model in the Jovian regime yield a strong prograde equatorial jet, and larger amplitude mid-latitude jets than the Scott and Polvani model. For both models, the slope of the non-zonal energy spectra is consistent with the classic Kolmogorov scaling, and the slope of the zonal energy spectra is consistent with the much steeper spectrum observed for Jupiter. We also perform simulations of the thermal shallow water equations for Neptunian parameter values, with a radiative relaxation time scale calculated for the same 25 mbar pressure level we used for Jupiter. These Neptunian simulations reproduce the broad, retrograde equatorial jet and prograde mid-latitude jets seen in observations. The much longer radiative time scale for the colder planet Neptune explains the transition from a prograde to a retrograde equatorial jet, while the broader jets are due to the deformation radius being a larger fraction of the planetary radius.

The vast majority of lattice Boltzmann algorithms produce a non-Galilean invariant viscous stress. This defect arises from the absence of a term in the third moment, the equilibrium heat flow tensor, proportional to the cube of the fluid velocity. This moment cannot be specified independently of the lower moments on the standard lattices such as D2Q9, D3Q15, D3Q19 or D3Q27. A partial correction has recently been demonstrated that restores some of these missing cubic terms on the D2Q9 and D3Q27 tensor product lattices. This correction restores Galilean invariance for shear flows aligned with the coordinate axes, but flows inclined at arbitrary angles may show larger errors than before. These remaining errors are due to the diagonal terms of the equilibrium heat flow tensor, which cannot be corrected on standard lattices. However, the remaining errors may be largely absorbed by introducing a matrix collision operator with velocity-dependent collision rates for the diagonal components of the momentum flux tensor. This completely restores Galilean invariance for flows with uniform density, and in general reduces the magnitude of the defect in Galilean invariance from Mach number cubed to Mach number to the fifth power. The effectiveness of the resulting algorithm is demonstrated by comparisons with the standard and partially corrected lattice Boltzmann algorithms for two- and three-dimensional flows.

The kinetic theory of gases implies an independent evolutionequation for the momentum flux tensor that closely resembles anevolution equation for the elastic stress in continuum descriptions ofviscoelastic liquids. However, kinetic theory leads to anonobjective convected derivative for the evolution of thedeviatoric stress, and a fixed relation between thestress relaxation rate and the viscosity. We show that simulations of freelydecaying shear flow using the standard two-dimensional latticeBoltzmann kinetic model develop a tangential stress consistent with thisnonobjective convected derivative, and this fixed relation betweenparameters. By contrast, viscoelastic liquids are typically modeled byan upper convected derivative, and with two independent parametersfor the viscosity and stress relaxation rate. Although we are unableto obtain an upper convected derivative from kinetic theory with ascalar distribution function, we show that introducing a generallinear coupling to a second stress tensor yields the linear Jeffreysviscoelastic model with three independent parameters in the incompressible limit. Unlike previouswork, we do not attempt to represent the additional stress throughmoments of additional distribution functions, but treat it only asan abstract tensor that couples to the corresponding tensorialmoment of the hydrodynamic distribution functions. This greatlysimplifies the derivation, and the implementation of flows driven bybody forces. The utility of the approach is demonstrated throughsimulations of Stokes' second problem for an oscillating boundary,of the four-roller mill, and of three-dimensional Arnold--Beltrami--Childress and Taylor--Green flows.