Quick von Neumann stability analysis

von Neumann stability analysis is a great tool to assess the stabilty of discretizing schemes for PDEs. But too often, imho, the discussion is too convoluted. Here, I try to provide a shortcut.

In a nutshell:

  • The standard approach involves Fourier techniques, involving (of course) complex numbers
  • The real part of these numbers is analysed, with some trigonometric expression resulting, identifying the troublesome modes
  • I claim this mode can be identified in advance, which makes the whole Fourier procedure unnecessary

BTW: it’s pronounced “fon no ee man”


Continue reading


Discrete Fourier transform and Fourier series

This is quite silly, but the relationship between the discrete Fourier transform (DFT) and the Fourier series (FS) is clouded by annoying factors. I will try to connect both in this article. The motivation is to employ DFT techniques in a computer simulation. In the latter, one usually has a finite simulation box, which makes Fourier series the most interesting (a connection to the Fourier transform may also be made, see below).

Continue reading

Van der Waals’ square gradient theory

OK, I have encountered this theory again, after many years. The idea is to describe separation between two phases as the minimum of a free energy with respect to an order parameter \phi:

g(\phi) = \frac a2 \phi^2 + \frac b4 \phi^2

When becomes negative, the minimum of g changes from 0 to other values, \pm \phi_0. The function then has the celebrated double minimum feature, which features prominently in many symmetry-breaking theories, including of course the appearance of particle mass mass and, you know, the Big Bang and the Universe.

But here we are just considering phase separation in materials. The interface between two coexisting phases must have some associated cost, and the simplest (lowest order) way to include it is by introducing a total free energy functional

F = \int dr \,  f(\phi, \nabla\phi)

f(\phi) = g(\phi) + \frac c2 (\nabla\phi)^2 .

This is also called a London-Linzburg-Landau free energy, also appearing in their theory of superconductors.

Now, parameters ab, and c are not easy to measure (or, at least, estimate) experimentally, but they are related to: the surface tension, the width of the interface, and the magnitude of the bulk equilibrium order parameter (i.e. \phi_0 ). Here I show how to obtain it in a slightly more general setting, since I was not able to find it on the internet (it can be found e.g, in the book by Rowlinson & Widom).

Continue reading

Taylor-Green vortex sheet, reduced units

The Taylor-Green vortex sheet is a solution to the 2D Navier-Stokes equations for an incompressible Newtonian fluid:

\frac{d \mathbf{u}}{d t}= \nu \nabla^2 \mathbf{u} - \nabla p/\rho ,

where \mathbf{u} is the velocity field, p is the pressure, \nu=\mu/\rho is the kinematic viscosity, and \rho is the fixed density of the fluid. The time derivative is a total derivative:

\frac{d \mathbf{u}}{d t} =  \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u}

It is common to choose parameters that simplify the equations, but that can obscure the role of the different parameters. In the following, I provide expressions with all relevant parameters included, with their physical dimensions. I later pass to dimensionless, or reduced, units, in terms of the Reynolds and Courant numbers.


Continue reading

Quick and dirty Rayleigh-Taylor instability

The Rayleigh-Taylor instability is a well-known benchmarck for CFD codes. The idea is to start with two phases, on on top of the other, the lighter one being underneath. The interface is slightly perturbed, and this plume appears. I describe a quick and dirty way of getting this instability.

Continue reading

Sound waves with attenuation

Just a simple derivation of the role of attenuation in the standard sound wave equation. Original work: Stokes, 1845.

Starting with the Navier-Stokes momentum equation

\frac{\partial }{\partial t} \mathbf{u} + \mathbf{u} \nabla \mathbf{u} = - \frac{1}{\rho} \nabla p + \frac{\mu}{\rho} \nabla^2 \mathbf{u} + \left(\frac{\lambda+\mu}{\rho}\right)\nabla (\nabla\cdot\mathbf{u}) ,

where \lambda is a Lamé viscosity coefficient. The bulk viscosity coeficient  is defined as  \zeta = \lambda + (2/3) \mu. The last term  is often neglected, even in compressible flow, but sound attenuation is one of the few cases where it may have some influence. All viscosities are assumed to be constant, but in this case this is a safe assumption, since we are going to assume small departures about equilibrium values.

Continue reading

Sparse Poisson problem in eigen

OwlgenBack to scientific computing. Lately, I have been using the Eigen libraries for linear algebra. They were recommended by the CGAL project, and they indeed share a common philosophy. Thanks to the rather high C++ level they can accomplish this sort of magic:

  int n = 100;
  VectorXd x(n), b(n);
  SpMat A(n,n);



  // solve Ax = b
  Eigen::BiCGSTAB<SpMat> solver;
  //ConjugateGradient<SpMat> solver;


  x = solver.solveWithGuess(b,x0);

Notice that A is a sparse matrix! I am next describing how to use this in order to solve the 1D Poisson equation.

Continue reading