Parallel Spectral Numerical Methods/Incompressible Magnetohydrodynamics

= Incompressible Magnetohydrodynamics =

Background
The equations of incompressible magnetohydrodynamics are similar to the Navier-Stokes equations, but can show a variety of other phenomena. This section assumes that the reader is familiar with the material covered for the Navier-Stokes case here. Their mathematical investigation of the incompressible magnetohydrodynamics equations is also in as poor a state as that of the Navier Stokes equations. A physical introduction to these equations can be found in Carbone and Pouquet, Chandrasekhar , Davidson , Gerbeau, Le Bris and Lelièvre , Schnak and Verma. A more general overview can be found here. Under the assumption that fluid velocities are significantly smaller than the speed of light, and the fluid is nearly incompressible, dimensionless equations are given as

In these equations, $$\mathbf{u}(x,y,z)=(u,v,w)$$ is the fluid velocity with components in the $$x$$, $$y$$ and $$z$$ directions, $$\mathbf{b}(x,y,z)=(bx,by,bz)$$ is the magnetic field with components in the $$x$$, $$y$$ and $$z$$ directions, $$p$$ is pressure field. Equation $$ represents conservation of momentum, eq. $$ corresponds to Maxwell's equations for nonrelativistic flows for which the time evolution of the electric field is neglected (see Verma ), eq. $$ is the continuity equation which represents conservation of mass for an incompressible fluid and eq. $$ represents the fact that there are no magnetic monopole sources.

The Two Dimensional Case
We will first consider the two-dimensional case. A difficulty in simulating the incompressible magnetohydrodynamic equations is the numerical satisfaction of the constraints in eqs. $$ and $$, these are sometimes referred to as a divergence free conditions or a solenoidal constraints. To automatically satisfy these constraints in two dimensions, where

$$\mathbf u(x,y)=\left(u(x,y),v(x,y)\right)$$

and

$$\mathbf b(x,y)=\left(bx(x,y),by(x,y)\right)$$

it is possible to re-write the equations using a potential formulation. In this case, we let

$$u=\frac{\partial \psi}{\partial y}$$

$$v=-\frac{\partial \psi}{\partial x}$$

$$bx=\frac{\partial \phi}{\partial y}$$

and

$$by=-\frac{\partial \phi}{\partial x}$$

where $$\psi(x,y)$$ is the streamfunction and $$\phi(x,y)$$ the magnetic potential. Note that

$$\nabla\cdot\mathbf u =\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}= \frac{\partial^2\psi}{\partial x \partial y} - \frac{\partial^2\psi}{\partial y \partial x}=0$$,

and

$$\nabla\cdot\mathbf b =\frac{\partial bx}{\partial x}+\frac{\partial by}{\partial y}= \frac{\partial^2\phi}{\partial x \partial y} - \frac{\partial^2\phi}{\partial y \partial x}=0$$,

so eqs. $$ and $$ are automatically satisfied. Making this change of variables, we obtain a reduced system of two coupled partial differential equation by taking the curl of the momentum equation, eq. $$ and the magnetic transport eq. $$. We define the fluid vorticity $$\omega$$, so that

$$\omega=\nabla\times\mathbf u= \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}=-\Delta \psi$$

and the magnetic vorticity so that

$$\alpha=\nabla\times\mathbf b= \frac{\partial by}{\partial x}-\frac{\partial bx}{\partial y}=-\Delta \phi$$

and eqs. $$ and $$ become

and

where $$fx$$ and $$fy$$ represent the $$x$$ and $$y$$ components of the force $$\mathbf f$$.

The Three Dimensional Case
In the three dimensional case, a potential formulation of the equations is not as convenient. Consequently, it is convenient to introduce the notion of a magnetic pressure to allow for a simple algorithm for enforcing the divergence free magnetic field constraint. The equations then become,

and

both of which can be discretized by the implicit midpoint rule. Here $$Re$$ denotes the fluids Reynolds number and $$Rem$$ the magnetic Reynolds number. The fluid pressure is given by

$$ \Delta^{-1}\left[\nabla\cdot \left(\mathbf u\cdot\nabla\mathbf u\right) -\left(\mathbf{b} \cdot \nabla \mathbf{b}\right) \right] - \mathbf{b}\cdot \mathbf{b} $$

and the magnetic ``pressure" is given by

$$\Delta^{-1}\left[\nabla\cdot \left(\mathbf u\cdot\nabla\mathbf b\right) -\left(\mathbf{b} \cdot \nabla \mathbf{u} \right) \right] $$

Serial Programs
We first write Matlab programs to demonstrate how to solve these equations on a single processor. The first program uses implicit midpoint rule timestepping to solve the two-dimensional magnetohydrodynamic equations and is in listing $$. In this program, an additional equation is solved for a passive scalar with a diffusion constant $$D_{\theta}$$

A Matlab program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download

A Python program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download

The second program uses the implicit midpoint rule to do timestepping for the three-dimensional incompressible Magnetohydrodynamics equations and it is in listing $$.

A Matlab program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download

A Python program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download

Fortran Programs
The programs in this section are based on those for the Navier-Stokes equations in https://en.wikibooks.org/wiki/Parallel_Spectral_Numerical_Methods/The_Two-_and_Three-Dimensional_Navier-Stokes_Equations.

A serial Fortran program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download A makefile to compile the program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download

Parallel Fortran Programs
To make a parallel program, the library 2DECOMP&FFT is used for the domain decomposition and the Fourier transforms. 2DECOMP&FFT is available from http://www.2decomp.org/index.html. The website includes some examples which indicate how this library should be used, in particular the sample code at http://www.2decomp.org/case_study1.html is a very helpful indication of how one converts a code that uses FFTW to one that uses MPI and the aforementioned library. To get the programs to work, first compile 2DECOMP&FFT, update the link in the makefile, and then compile the parallel Fortran program.

A parallel Fortran program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download

A subroutine to save real array data for the parallel MPI Fortran program to solve the 3D incompressible magnetohydrodynamics equations Code download

A makefile to compile the program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download

Other Publicly Available Fourier Spectral Programs for Solving the Magnetohydrodynamics Equations
One other publicly available Fast Fourier Transform based program for solving the incompressible magnetohydrodynamics equations is:

Tarang available from http://turbulencehub.org/index.php/codes/tarang/

Exercises

 * 1) The three dimensional program we have written is not the most efficient since one can use a real to complex transform to halve the work done. Implement a real to complex transform in one of the Navier-Stokes programs.
 * 2) The programs we have written can also introduce some aliasing errors. By reading a book on spectral methods, such as Canuto et al. , find out what aliasing errors are. Explain why the strategy explained in Johnstone can reduce aliasing errors.
 * 3) See if you can find other open source Fast Fourier Transform based solvers for the incompressible magnetohydrodynamics equations and add them to the list above.
 * 4) There are suggestions that the incompressible magnetohydrodynamics equations can exhibit singular behavior by having one of the norms of the derivatives of the magnetic or velocity field becoming infinite. Experiment with simulations to see if you can find evidence for this. One of the first papers reporting computational experiments in this area is Orszag and Tang.
 * 5) Translate one of the programs above to another language. The following site may help, http://rosettacode.org/wiki/Rosetta_Code
 * 6) The full magnetohydrodynamic equations are

where $$\mathbf{U}$$ is the velocity, $$\mathbf{p}$$ is the pressure, $$\mathbf{J}$$ is the current density, $$\mathbf{B}$$ is the magnetic field and $$\mathbf{H}$$ is the electric field. Simulate these full equations in two and three dimensions - it is at present unknown whether solutions to these equations exist for all time in three dimensions see Germain, Ibrahim and Masmoudi, though Masmoudi shows that some solutions are well behaved in two dimensions. Analysis of the reduced equations is in Duvaut and Lions and can also be found in Gerbeau, Le Bris, and Lelièvre. An example program is in listing $$ A Matlab program which finds a numerical solution to the 3D full magnetohydrodynamic equations Code download
 * 1) A compressible magnetohydrodynamic model is

where $$\mathbf{U}$$ is the velocity, $$\mathbf{p}$$ is the pressure, $$\mathbf{J}$$ is the current density, $$\mathbf{B}$$ is the magnetic field, $$\mathbf{H}$$ is the electric field, $$\mathbf{E}$$ is the energy density, $$\mathbf{\tau}$$ is the viscous fluid stress, $$T$$ is the temperature field, $$C_v$$ is the constant volume specific heat capacity, $$R$$ is the ideal gas constant, $$\mathbf{q}$$ is the heat flux. Simulate these full equations in two and three dimensions.