# 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[1], Chandrasekhar[2], Davidson[3], Gerbeau, Le Bris and Lelièvre[4], Schnak[5] and Verma[6]. A more general overview can be found here [7]. 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

$\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p + \mathbf{b} \cdot \nabla \mathbf{b}-\nabla\left(\mathbf{b}\cdot\mathbf{b}\right) + \Delta \mathbf{u}$

( 1)

$\frac{\partial \mathbf{b}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{b} = \mathbf{b} \cdot \nabla \mathbf{u} + \Delta \mathbf{b}$

( 2)

$\nabla \cdot \mathbf{u} = 0$.

( 3)

$\nabla \cdot \mathbf{b} = 0$.

( 4)

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 1 represents conservation of momentum, eq. 2 corresponds to Maxwell's equations for nonrelativistic flows for which the time evolution of the electric field is neglected (see Verma[8]), eq. 3 is the continuity equation which represents conservation of mass for an incompressible fluid and eq. 4 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. 3 and 4 , 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. 3 and 4 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. 1 and the magnetic transport eq. 2 . 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. 1 and 2 become

$\frac{\partial \omega}{\partial t} + u\frac{\partial \omega}{\partial x} + v\frac{\partial \omega }{\partial y} = bx\frac{\partial \alpha}{\partial x} + by\frac{\partial \alpha }{\partial y} + \mu\Delta\omega+\frac{\partial fy}{\partial x}- \frac{\partial fx}{\partial y}$

( 5)

$\frac{\partial \alpha}{\partial t} + u\frac{\partial \alpha}{\partial x} + v\frac{\partial \alpha }{\partial y} =bx\frac{\partial \omega}{\partial x} + by\frac{\partial \omega }{\partial y} + \eta\Delta\alpha$

( 6)

$\Delta \psi = -\omega$.

( 7)

and

$\Delta \phi = -\alpha$.

( 8)

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,

$\frac{\partial \mathbf u}{\partial t}=\frac{1}{\text{Re}}\Delta \mathbf u + \mathbf{b} \cdot \nabla \mathbf{b} - \mathbf u\cdot\nabla\mathbf u - \nabla\left(\mathbf{b}\cdot \mathbf{b}\right) +\nabla \Delta^{-1}\left[\nabla\cdot \left(\mathbf u\cdot\nabla\mathbf u\right) - \left(\mathbf{b} \cdot \nabla \mathbf{b}\right) \right]$,

( 9)

and

$\frac{\partial \mathbf b}{\partial t}=\frac{1}{\text{Rem}}\Delta \mathbf b + \mathbf{b} \cdot \nabla \mathbf{u} - \mathbf u\cdot\nabla\mathbf b +\nabla \Delta^{-1}\left[\nabla\cdot \left(\mathbf u\cdot\nabla\mathbf b\right) -\left(\mathbf{b} \cdot \nabla \mathbf{u} \right) \right]$,

( 10)

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 A . In this program, an additional equation is solved for a passive scalar with a diffusion constant $D_{\theta}$

$\frac{\partial \theta}{\partial t} + \mathbf{u} \cdot \nabla \theta = D_{\theta} \Delta \theta$

( 11)

( A)

A Matlab program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download
% Numerical solution of the 2D incompressible Magnetohydrodynamics equations
% on a Square Domain [0,1]x[0,1] using a Fourier pseudo-spectral method
% and Implicit midpoint rule timestepping. A passive scalar is also advected
% by the flow
%
%Periodic free-slip boundary conditions and Initial conditions:
%u(x,y,0)=sin(2*pi*x)cos(2*pi*y)
%v(x,y,0)=-cos(2*pi*x)sin(2*pi*y)
%Analytical Solution:
%u(x,y,t)=sin(2*pi*x)cos(2*pi*y)exp(-8*pi^2*t/Re)
%v(x,y,t)=-cos(2*pi*x)sin(2*pi*y)exp(-8*pi^2*t/Re)
% also consider an advection field
% \theta_t+u\theta_x+v\theta_y=Dtheta(\theta_{xx}+\theta_{yy})
clear all; format compact; format short; clc; clf;

Re=100;%Reynolds number
Rem=100; % Magnetic Reynolds number
Dtheta=0.01;%scalar diffusion constant
%grid
Nx=64; h=1/Nx; x=h*(1:Nx);
Ny=64; h=1/Ny; y=h*(1:Ny)';
[xx,yy]=meshgrid(x,y);

%initial conditions for velocity field
u=sin(2*pi*xx).*cos(2*pi*yy);
v=-cos(2*pi*xx).*sin(2*pi*yy);
u_y=-2*pi*sin(2*pi*xx).*sin(2*pi*yy);
v_x=2*pi*sin(2*pi*xx).*sin(2*pi*yy);
omega=v_x-u_y;
% initial magnetic potential field
alpha=sin(4*pi*xx).*cos(6*pi*yy);
% initial passive scalar field
theta=exp(-2.0*((4*pi*(xx-0.25)).^2+(4*pi*(yy-0.25)).^2));
dt=0.0025; t(1)=0; tmax=1.0;
nplots=ceil(tmax/dt);

%wave numbers for derivatives
k_x=2*pi*(1i*[(0:Nx/2-1)  0 1-Nx/2:-1]');
k_y=2*pi*(1i*[(0:Ny/2-1)  0 1-Ny/2:-1]);
k2x=k_x.^2;
k2y=k_y.^2;

%wave number grid for multiplying matricies
[kxx,kyy]=meshgrid(k2x,k2y);
[kx,ky]=meshgrid(k_x,k_y);

% use a high tolerance so time stepping errors
% are not dominated by errors in solution to nonlinear
% system
tol=10^(-10);
%compute \hat{Phi}^{n+1,k+1}
alphahat=fft2(alpha);
phihat=-alphahat./(kxx+kyy);

%NOTE: kxx+kyy has to be zero at the following points to avoid a
% discontinuity. However, we suppose that the streamfunction has
% mean value zero, so we set them equal to zero
phihat(1,1)=0;
phihat(Nx/2+1,Ny/2+1)=0;
phihat(Nx/2+1,1)=0;
phihat(1,Ny/2+1)=0;

%computes {\psi}_x by differentiation via FFT
dphix = real(ifft2(phihat.*kx));
%computes {\psi}_y by differentiation via FFT
dphiy = real(ifft2(phihat.*ky));
% components of magnetic field
bx=dphiy;
by=-dphix;

%compute \hat{\omega}^{n+1,k}
omegahat=fft2(omega); alphahat=fft2(alpha);
thetatotal(1)=sum(sum(theta.^2))*h*h;
for i=1:nplots
chg=1;
% save old values
uold=u; vold=v;  omegaold=omega; omegacheck=omega;
omegahatold=omegahat; thetaold=theta; thetacheck=theta;
alphahatold=alphahat; alphaold=alpha; alphacheck=alpha;
bxold=bx; byold=by;
while chg>tol
% Fluid field
%nonlinear {n+1,k}
nonlinhat=0.25*fft2((u+uold).*ifft2((omegahat+omegahatold).*kx)+...
(v+vold).*ifft2((omegahat+omegahatold).*ky)-...
(bx+bxold).*ifft2((alphahat+alphahatold).*kx)-...
(by+byold).*ifft2((alphahat+alphahatold).*ky));

%Implicit midpoint rule timestepping
omegahat=((1/dt + 0.5*(1/Re)*(kxx+kyy)).*omegahatold-nonlinhat)...
./(1/dt -0.5*(1/Re)*(kxx+kyy));

%compute \hat{\psi}^{n+1,k+1}
psihat=-omegahat./(kxx+kyy);

%NOTE: kxx+kyy has to be zero at the following points to avoid a
% discontinuity. However, we suppose that the streamfunction has
% mean value zero, so we set them equal to zero
psihat(1,1)=0;
psihat(Nx/2+1,Ny/2+1)=0;
psihat(Nx/2+1,1)=0;
psihat(1,Ny/2+1)=0;

%computes {\psi}_x by differentiation via FFT
dpsix = real(ifft2(psihat.*kx));
%computes {\psi}_y by differentiation via FFT
dpsiy = real(ifft2(psihat.*ky));

u=dpsiy;    %u^{n+1,k+1}
v=-dpsix;   %v^{n+1,k+1}

% magnetic field
nonlinhat=0.25*fft2((u+uold).*ifft2((alphahat+alphahatold).*kx)+...
(v+vold).*ifft2((alphahat+alphahatold).*ky)-...
(bx+bxold).*ifft2((omegahat+omegahatold).*kx)-...
(by+byold).*ifft2((omegahat+omegahatold).*ky));

%Implicit midpoint rule timestepping
alphahat=((1/dt + 0.5*(1/Re)*(kxx+kyy)).*alphahatold-nonlinhat)...
./(1/dt -0.5*(1/Rem)*(kxx+kyy));

%compute \hat{\psi}^{n+1,k+1}
phihat=-alphahat./(kxx+kyy);

%NOTE: kxx+kyy has to be zero at the following points to avoid a
% discontinuity. However, we suppose that the streamfunction has
% mean value zero, so we set them equal to zero
phihat(1,1)=0;
phihat(Nx/2+1,Ny/2+1)=0;
phihat(Nx/2+1,1)=0;
phihat(1,Ny/2+1)=0;

%computes {\psi}_x by differentiation via FFT
dphix = real(ifft2(phihat.*kx));
%computes {\psi}_y by differentiation via FFT
dphiy = real(ifft2(phihat.*ky));

bx=dphiy;    %u^{n+1,k+1}
by=-dphix;   %v^{n+1,k+1}

thetax=0.5*ifft2(kx.*fft2(theta+thetaold));
thetay=0.5*ifft2(ky.*fft2(theta+thetaold));
theta=ifft2((fft2(thetaold-dt*0.5*((uold+u).*thetax+(vold+v).*thetay))+...
dt*0.5*Dtheta*(kxx+kyy).*fft2(thetaold))./(1-dt*0.5*Dtheta*(kxx+kyy)));

%\omega^{n+1,k+1}
omega=ifft2(omegahat);
% check for convergence
chg=max(max(abs(omega-omegacheck))) + max(max(abs(theta-thetacheck)))+...
max(max(abs(alpha-alphacheck)))

% store omega and theta to check for convergence of next iteration
omegacheck=omega; thetacheck=theta; alphacheck=alpha;
end
t(i+1)=t(i)+dt;
thetatotal(i+1)=sum(sum(theta.^2))*h*h;
figure(1);
subplot(2,2,1);
pcolor(xx,yy,omega);  shading interp; xlabel x; ylabel y;
title(['Fluid Vorticity, Time ',num2str(t(i+1))]); colorbar; drawnow;
subplot(2,2,2);
pcolor(xx,yy,alpha); shading interp;  xlabel x; ylabel y;
title(['Alpha, Time ',num2str(t(i+1))]); colorbar; drawnow;
subplot(2,2,3);
pcolor(xx,yy,real(ifft2(psihat))); shading interp; xlabel x; ylabel y;
title(['Psi, Time ',num2str(t(i+1))]); colorbar; drawnow;
subplot(2,2,4);
pcolor(xx,yy,real(theta)); shading interp;  xlabel x; ylabel y;
title(['Theta, Time ',num2str(t(i+1))]); colorbar; drawnow;
end
figure(2); clf;
subplot(2,1,1); plot(t,thetatotal,'r-'); xlabel time; ylabel('Total theta');
subplot(2,1,2); semilogy(t,abs(1-thetatotal./thetatotal(1)),'r-');
xlabel time; ylabel('Total theta Error');


( A1)

A Python program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download
#!/usr/bin/env python
"""
Numerical solution of the 2D incompressible Magnetohydrodynamics equations on a
Square Domain [0,1]x[0,1] using a Fourier pseudo-spectral method
and Implicit Midpoint rule timestepping.

A scalar field is also advected by the flow
\theta_t+u\theta_x+v\theta_y=Dtheta(\theta_{xx}+\theta_{yy})

Periodic free-slip boundary conditions and Initial conditions:
u(x,y,0)=sin(2*pi*x)cos(2*pi*y)
v(x,y,0)=-cos(2*pi*x)sin(2*pi*y)

"""

import math
import numpy
import matplotlib.pyplot as plt
from mayavi import mlab
import time

# Grid
N=64; h=1.0/N
x = [h*i for i in xrange(1,N+1)]
y = [h*i for i in xrange(1,N+1)]
numpy.savetxt('x.txt',x)

xx=numpy.zeros((N,N), dtype=float)
yy=numpy.zeros((N,N), dtype=float)

for i in xrange(N):
for j in xrange(N):
xx[i,j] = x[i]
yy[i,j] = y[j]

dt=0.0025; t=0.0; tmax=1.0
#nplots=int(tmax/dt)
Re=100
Rem=100
Dtheta=0.001

u=numpy.zeros((N,N), dtype=float)
v=numpy.zeros((N,N), dtype=float)
u_y=numpy.zeros((N,N), dtype=float)
v_x=numpy.zeros((N,N), dtype=float)
omega=numpy.zeros((N,N), dtype=float)
theta=numpy.zeros((N,N), dtype=float)
thetaold=numpy.zeros((N,N), dtype=float)
thetax=numpy.zeros((N,N), dtype=float)
thetay=numpy.zeros((N,N), dtype=float)
thetacheck=numpy.zeros((N,N), dtype=float)
alpha=numpy.zeros((N,N), dtype=float)
alphaold=numpy.zeros((N,N), dtype=float)
bx=numpy.zeros((N,N), dtype=float)
bxold=numpy.zeros((N,N), dtype=float)
by=numpy.zeros((N,N), dtype=float)
byold=numpy.zeros((N,N), dtype=float)

# Initial conditions
for i in range(len(x)):
for j in range(len(y)):
u[i][j]=numpy.sin(2*math.pi*x[i])*numpy.cos(2*math.pi*y[j])
v[i][j]=-numpy.cos(2*math.pi*x[i])*numpy.sin(2*math.pi*y[j])
u_y[i][j]=-2*math.pi*numpy.sin(2*math.pi*x[i])*numpy.sin(2*math.pi*y[j])
v_x[i][j]=2*math.pi*numpy.sin(2*math.pi*x[i])*numpy.sin(2*math.pi*y[j])
omega[i][j]=v_x[i][j]-u_y[i][j]
theta[i][j]=numpy.exp(-2.0*((4*math.pi*(x[i]-0.3))**2+(4*math.pi*(y[j]-0.3))**2))
alpha[i][j]=numpy.sin(4*math.pi*x[i])*numpy.cos(6*math.pi*y[j])

src = mlab.imshow(xx,yy,theta,colormap='jet')
mlab.scalarbar(object=src)
mlab.xlabel('x',object=src)
mlab.ylabel('y',object=src)

# Wavenumber
k_x = 2*math.pi*numpy.array([complex(0,1)*n for n in range(0,N/2) \
+ [0] + range(-N/2+1,0)])
k_y=k_x

kx=numpy.zeros((N,N), dtype=complex)
ky=numpy.zeros((N,N), dtype=complex)
kxx=numpy.zeros((N,N), dtype=complex)
kyy=numpy.zeros((N,N), dtype=complex)

for i in xrange(N):
for j in xrange(N):
kx[i,j] = k_x[i]
ky[i,j] = k_y[j]
kxx[i,j] = k_x[i]**2
kyy[i,j] = k_y[j]**2

tol=10**(-10)
psihat=numpy.zeros((N,N), dtype=complex)
phihat=numpy.zeros((N,N), dtype=complex)
alphahat=numpy.zeros((N,N), dtype=complex)
alphahatold=numpy.zeros((N,N), dtype=complex)
omegahat=numpy.zeros((N,N), dtype=complex)
omegahatold=numpy.zeros((N,N), dtype=complex)
thetahat=numpy.zeros((N,N), dtype=complex)
thetahatold=numpy.zeros((N,N), dtype=complex)
nlhat=numpy.zeros((N,N), dtype=complex)
dpsix=numpy.zeros((N,N), dtype=float)
dpsiy=numpy.zeros((N,N), dtype=float)
dphix=numpy.zeros((N,N), dtype=float)
dphiy=numpy.zeros((N,N), dtype=float)
omegacheck=numpy.zeros((N,N), dtype=float)
omegaold=numpy.zeros((N,N), dtype=float)
temp=numpy.zeros((N,N), dtype=float)
omegahat=numpy.fft.fft2(omega)
thetahat=numpy.fft.fft2(theta)
alphahat=numpy.fft.fft2(alpha)

while (t<=tmax):
chg=1.0

# Save old values
uold=u
vold=v
omegaold=omega
omegacheck=omega
omegahatold = omegahat
thetaold=theta
thetahatold=thetahat
thetacheck=theta
bxold=bx
byold=by
alphahatold=alphahat
alphaold=alpha
alphacheck=alpha
while(chg>tol):
# fluid field
# nolinear {n+1,k}
nlhat=0.25*numpy.fft.fft2((u+uold)*numpy.fft.ifft2((omegahat+omegahatold)*kx)+\
(v+vold)*numpy.fft.ifft2((omegahat+omegahatold)*ky)-\
(bx+bxold)*numpy.fft.ifft2((alphahat+alphahatold)*kx)-\
(by+byold)*numpy.fft.ifft2((alphahat+alphahatold)*ky))

# Implicit midpoint rule timestepping
omegahat=((1/dt + 0.5*(1/Re)*(kxx+kyy))*omegahatold \
-nlhat) \
/(1/dt -0.5*(1/Re)*(kxx+kyy))

psihat=-omegahat/(kxx+kyy)
psihat[0][0]=0
psihat[N/2][N/2]=0
psihat[N/2][0]=0
psihat[0][N/2]=0

dpsix = numpy.real(numpy.fft.ifft2(psihat*kx))
dpsiy = numpy.real(numpy.fft.ifft2(psihat*ky))
u=dpsiy
v=-1.0*dpsix

omega=numpy.real(numpy.fft.ifft2(omegahat))
# magnetic field
nlhat=0.25*numpy.fft.fft2((u+uold)*numpy.fft.ifft2((alphahat+alphahatold)*kx)+\
(v+vold)*numpy.fft.ifft2((alphahat+alphahatold)*ky)-\
(bx+bxold)*numpy.fft.ifft2((omegahat+omegahatold)*kx)-\
(by+byold)*numpy.fft.ifft2((omegahat+omegahatold)*ky))

# Implicit midpoint rule timestepping
alphhat=((1/dt + 0.5*(1/Rem)*(kxx+kyy))*alphahatold \
-nlhat) \
/(1/dt -0.5*(1/Rem)*(kxx+kyy))

phihat=-alphahat/(kxx+kyy)
phihat[0][0]=0
phihat[N/2][N/2]=0
phihat[N/2][0]=0
phihat[0][N/2]=0

dphix = numpy.real(numpy.fft.ifft2(phihat*kx))
dphiy = numpy.real(numpy.fft.ifft2(phihat*ky))
bx=dphiy
by=-1.0*dphix

alpha=numpy.real(numpy.fft.ifft2(alphahat))

# passive scalar
thetax=0.5*numpy.real(numpy.fft.ifft2(kx*(thetahat+thetahatold)))
thetay=0.5*numpy.real(numpy.fft.ifft2(ky*(thetahat+thetahatold)))
thetahat= numpy.fft.fft2(thetaold-dt*0.5*((uold+u)*thetax+(vold+v)*thetay) )
theta=numpy.real(numpy.fft.ifft2( thetahat))

temp=abs(omega-omegacheck)
chg=numpy.max(temp)
temp=abs(theta-thetacheck)
chg=chg+numpy.max(temp)
print(chg)
omegacheck=omega
thetacheck=theta
t+=dt
src.mlab_source.scalars = theta


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

( B)

A Matlab program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download
% A program to solve the 3D incompressible magnetohydrodynamic equations
% with periodic boundary
% conditions. The program is based on the Orszag-Patterson algorithm as
% documented on pg. 98 of C. Canuto, M.Y. Hussaini, A. Quarteroni and
% T.A. Zhang "Spectral Methods: Evolution to Complex Geometries and
% Applications to Fluid Dynamics" Springer (2007)
%
% The Helmholtz decomposition is used to project the magnetic field onto a
% divergence free subspace. Initial work on this has been done with Damian
% San Roman Alerigi
%
% For plotting, the program uses vol3d, which is available at:
%
% http://www.mathworks.com/matlabcentral/fileexchange/22940-vol3d-v2
%
clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
'defaultaxesfontweight','bold')

% set up grid
tic
Lx = 1;         % period  2*pi*L
Ly = 1;         % period  2*pi*L
Lz = 1;         % period  2*pi*L
Nx = 32;        % number of harmonics
Ny = 32;        % number of harmonics
Nz = 32;        % number of harmonics
Nt = 50;        % number of time slices
dt = 2/Nt;    % time step
t=0;            % initial time
Re = 100.0;       % Reynolds number
Rem = 100.0;      % Magnetic Reynolds number
tol=10^(-10);
% initialise variables
x = (2*pi/Nx)*(-Nx/2:Nx/2 -1)'*Lx;          % x coordinate
kx = 1i*[0:Nx/2-1 0 -Nx/2+1:-1]'/Lx;        % wave vector
y = (2*pi/Ny)*(-Ny/2:Ny/2 -1)'*Ly;          % y coordinate
ky = 1i*[0:Ny/2-1 0 -Ny/2+1:-1]'/Ly;        % wave vector
z = (2*pi/Nz)*(-Nz/2:Nz/2 -1)'*Lz;          % y coordinate
kz = 1i*[0:Nz/2-1 0 -Nz/2+1:-1]'/Lz;        % wave vector
[xx,yy,zz]=meshgrid(x,y,z);
[kxm,kym,kzm]=meshgrid(kx,ky,kz);
[k2xm,k2ym,k2zm]=meshgrid(kx.^2,ky.^2,kz.^2);

% initial conditions for Taylor-Green vortex
% theta=0;
% u=(2/sqrt(3))*sin(theta+2*pi/3)*sin(xx).*cos(yy).*cos(zz);
% v=(2/sqrt(3))*sin(theta-2*pi/3)*cos(xx).*sin(yy).*cos(zz);
% w=(2/sqrt(3))*sin(theta)*cos(xx).*cos(yy).*sin(zz);

% initial condition
sl=1; sk=1; sm=1; lamlkm=sqrt(sl.^2+sk.^2+sm.^2);
u=-0.5*(lamlkm*sl*cos(sk*xx).*sin(sl*yy).*sin(sm.*zz)...
+sm*sk*sin(sk*xx).*cos(sl*yy).*cos(sm.*zz))...
.*exp(-t*(lamlkm^2)/Re);

v=0.5*(lamlkm*sk*sin(sk*xx).*cos(sl*yy).*sin(sm.*zz)...
-sm*sl*cos(sk*xx).*sin(sl*yy).*cos(sm.*zz))...
*exp(-t*(lamlkm^2)/Re);

w=cos(sk*xx).*cos(sl*yy).*sin(sm*zz)*exp(-t*(lamlkm^2)/Re);

uhat=fftn(u);
vhat=fftn(v);
what=fftn(w);

ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm);
vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm);
wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm);

bx=sin(sl*yy).*sin(sm.*zz);
by=sin(sm.*zz);
bz=cos(sk*xx).*cos(sl*yy);

bxhat=fftn(bx);
byhat=fftn(by);
bzhat=fftn(bz);

bxx=ifftn(bxhat.*kxm);bxy=ifftn(bxhat.*kym);bxz=ifftn(bxhat.*kzm);
byx=ifftn(byhat.*kxm);byy=ifftn(byhat.*kym);byz=ifftn(byhat.*kzm);
bzx=ifftn(bzhat.*kxm);bzy=ifftn(bzhat.*kym);bzz=ifftn(bzhat.*kzm);

% calculate fluid vorticity for plotting
omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy;
omegatot=omegax.^2+omegay.^2+omegaz.^2;
% calculate magnetic vorticity for plotting
omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy;
omegabtot=omegabx.^2+omegaby.^2+omegabz.^2;
figure(1); clf; n=0;
subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]);
H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]);
H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

for n=1:Nt
uold=u; uxold=ux; uyold=uy; uzold=uz;
vold=v; vxold=vx; vyold=vy; vzold=vz;
wold=w; wxold=wx; wyold=wy; wzold=wz;

bxold=bx; bxxold=bxx; bxyold=bxy; bxzold=bxz;
byold=by; byxold=byx; byyold=byy; byzold=byz;
bzold=bz; bzxold=bzx; bzyold=bzy; bzzold=bzz;

rhsuhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*uhat;
rhsvhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*vhat;
rhswhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*what;

rhsbxhatfix=(1/dt+(0.5/Rem)*(k2xm+k2ym+k2zm)).*bxhat;
rhsbyhatfix=(1/dt+(0.5/Rem)*(k2xm+k2ym+k2zm)).*byhat;
rhsbzhatfix=(1/dt+(0.5/Rem)*(k2xm+k2ym+k2zm)).*bzhat;

chg=1; t=t+dt;
while (chg>tol)
nonlinu=0.25*((u+uold).*(ux+uxold)...
+(v+vold).*(uy+uyold)...
+(w+wold).*(uz+uzold)...
-(bx+bxold).*(bxx+bxxold)...
-(by+byold).*(bxy+bxyold)...
-(bz+bzold).*(bxz+bxzold));
nonlinv=0.25*((u+uold).*(vx+vxold)...
+(v+vold).*(vy+vyold)...
+(w+wold).*(vz+vzold)...
-(bx+bxold).*(byx+byxold)...
-(by+byold).*(byy+byyold)...
-(bz+bzold).*(byz+byzold));
nonlinw=0.25*((u+uold).*(wx+wxold)...
+(v+vold).*(wy+wyold)...
+(w+wold).*(wz+wzold)...
-(bx+bxold).*(bzx+bzxold)...
-(by+byold).*(bzy+bzyold)...
-(bz+bzold).*(bzz+bzzold));
nonlinuhat=fftn(nonlinu);
nonlinvhat=fftn(nonlinv);
nonlinwhat=fftn(nonlinw);
phat=-1.0*(kxm.*nonlinuhat+kym.*nonlinvhat+kzm.*nonlinwhat)...
./(k2xm+k2ym+k2zm+0.1^13);
uhat=(rhsuhatfix-nonlinuhat-kxm.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
vhat=(rhsvhatfix-nonlinvhat-kym.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
what=(rhswhatfix-nonlinwhat-kzm.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm);
vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm);
wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm);
utemp=u; vtemp=v; wtemp=w;
u=ifftn(uhat); v=ifftn(vhat); w=ifftn(what);

nonlinu=0.25*((u+uold).*(bxx+bxxold)...
+(v+vold).*(bxy+bxyold)...
+(w+wold).*(bxz+bxzold)...
-(bx+bxold).*(ux+uxold)...
-(by+byold).*(uy+uyold)...
-(bz+bzold).*(uz+uzold));
nonlinv=0.25*((u+uold).*(byx+byxold)...
+(v+vold).*(byy+byyold)...
+(w+wold).*(byz+byzold)...
-(bx+bxold).*(vx+vxold)...
-(by+byold).*(vy+vyold)...
-(bz+bzold).*(vz+vzold));
nonlinw=0.25*((u+uold).*(bzx+bzxold)...
+(v+vold).*(bzy+bzyold)...
+(w+wold).*(bzz+bzzold)...
-(bx+bxold).*(wx+wxold)...
-(by+byold).*(wy+wyold)...
-(bz+bzold).*(wz+wzold));
nonlinuhat=fftn(nonlinu);
nonlinvhat=fftn(nonlinv);
nonlinwhat=fftn(nonlinw);
phat=-1.0*(kxm.*nonlinuhat+kym.*nonlinvhat+kzm.*nonlinwhat)...
./(k2xm+k2ym+k2zm+0.1^13);
bxhat=(rhsbxhatfix-nonlinuhat-kxm.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
byhat=(rhsbyhatfix-nonlinvhat-kym.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
bzhat=(rhsbzhatfix-nonlinwhat-kzm.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
bxx=ifftn(bxhat.*kxm);bxy=ifftn(bxhat.*kym);bxz=ifftn(bxhat.*kzm);
byx=ifftn(byhat.*kxm);byy=ifftn(byhat.*kym);byz=ifftn(byhat.*kzm);
bzx=ifftn(bzhat.*kxm);bzy=ifftn(bzhat.*kym);bzz=ifftn(bzhat.*kzm);
bxtemp=bx; bytemp=by; bztemp=bz;
bx=ifftn(bxhat); by=ifftn(byhat); bz=ifftn(bzhat);

chg=max(abs(utemp-u))   +max(abs(vtemp-v))  +max(abs(wtemp-w))+...
max(abs(bxtemp-bx)) +max(abs(bytemp-by))+max(abs(bztemp-bz));
end
% calculate vorticity for plotting
omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy;
omegatot=omegax.^2+omegay.^2+omegaz.^2;

% calculate magnetic vorticity for plotting
omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy;
omegabtot=omegabx.^2+omegaby.^2+omegabz.^2;
figure(1); clf;
subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]);
H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]);
H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

end


( B1)

A Python program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download
#!/usr/bin/env python
"""
A program to solve the 3D incompressible magnetohydrodynamics equations using the
implicit midpoint rule

The program is based on the Orszag-Patterson algorithm as documented on pg. 98
of C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zhang
"Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics"
Springer (2007)

The Helmholtz decomposition is used to project the magnetic field onto a divergence
free subspace. Initial work on this has been done with Damian San Roman Alerigi

website, in particular:
http://github.enthought.com/mayavi/mayavi/mlab.html
which was last checked on 6 April 2012

"""

import math
import numpy
from mayavi import mlab
import matplotlib.pyplot as plt
import time

# Grid
Lx=1.0 	    	# Period 2*pi*Lx
Ly=1.0 	    	# Period 2*pi*Ly
Lz=1.0 	    	# Period 2*pi*Lz
Nx=64 	    	# Number of harmonics
Ny=64	    	# Number of harmonics
Nz=64 	    	# Number of harmonics
Nt=20 	    	# Number of time slices
tmax=0.2   	# Maximum time
dt=tmax/Nt  	# time step
t=0.0	    	# initial time
Re=1.0      	# Reynolds number
Rem=1.0		# Magnetic Reynolds number
tol=0.1**(10) 	# tolerance for fixed point iterations

x = [i*2.0*math.pi*(Lx/Nx) for i in xrange(-Nx/2,1+Nx/2)]
y = [i*2.0*math.pi*(Ly/Ny) for i in xrange(-Ny/2,1+Ny/2)]
z = [i*2.0*math.pi*(Lz/Nz) for i in xrange(-Nz/2,1+Nz/2)]
k_x = (1.0/Lx)*numpy.array([complex(0,1)*n for n in range(0,Nx/2) \
+ [0] + range(-Nx/2+1,0)])
k_y = (1.0/Ly)*numpy.array([complex(0,1)*n for n in range(0,Ny/2) \
+ [0] + range(-Ny/2+1,0)])
k_z = (1.0/Lz)*numpy.array([complex(0,1)*n for n in range(0,Nz/2) \
+ [0] + range(-Nz/2+1,0)])

kxm=numpy.zeros((Nx,Ny,Nz), dtype=complex)
kym=numpy.zeros((Nx,Ny,Nz), dtype=complex)
kzm=numpy.zeros((Nx,Ny,Nz), dtype=complex)
k2xm=numpy.zeros((Nx,Ny,Nz), dtype=float)
k2ym=numpy.zeros((Nx,Ny,Nz), dtype=float)
k2zm=numpy.zeros((Nx,Ny,Nz), dtype=float)
xx=numpy.zeros((Nx,Ny,Nz), dtype=float)
yy=numpy.zeros((Nx,Ny,Nz), dtype=float)
zz=numpy.zeros((Nx,Ny,Nz), dtype=float)

for i in xrange(Nx):
for j in xrange(Ny):
for k in xrange(Nz):
kxm[i,j,k] = k_x[i]
kym[i,j,k] = k_y[j]
kzm[i,j,k] = k_z[k]
k2xm[i,j,k] = numpy.real(k_x[i]**2)
k2ym[i,j,k] = numpy.real(k_y[j]**2)
k2zm[i,j,k] = numpy.real(k_z[k]**2)
xx[i,j,k] = x[i]
yy[i,j,k] = y[j]
zz[i,j,k] = z[k]

# allocate arrays
u=numpy.zeros((Nx,Ny,Nz), dtype=float)
uold=numpy.zeros((Nx,Ny,Nz), dtype=float)
v=numpy.zeros((Nx,Ny,Nz), dtype=float)
vold=numpy.zeros((Nx,Ny,Nz), dtype=float)
w=numpy.zeros((Nx,Ny,Nz), dtype=float)
wold=numpy.zeros((Nx,Ny,Nz), dtype=float)

bx=numpy.zeros((Nx,Ny,Nz), dtype=float)
bxold=numpy.zeros((Nx,Ny,Nz), dtype=float)
by=numpy.zeros((Nx,Ny,Nz), dtype=float)
byold=numpy.zeros((Nx,Ny,Nz), dtype=float)
bz=numpy.zeros((Nx,Ny,Nz), dtype=float)
bzold=numpy.zeros((Nx,Ny,Nz), dtype=float)

utemp=numpy.zeros((Nx,Ny,Nz), dtype=float)
vtemp=numpy.zeros((Nx,Ny,Nz), dtype=float)
wtemp=numpy.zeros((Nx,Ny,Nz), dtype=float)

bxtemp=numpy.zeros((Nx,Ny,Nz), dtype=float)
bytemp=numpy.zeros((Nx,Ny,Nz), dtype=float)
bztemp=numpy.zeros((Nx,Ny,Nz), dtype=float)

omegax=numpy.zeros((Nx,Ny,Nz), dtype=float)
omegay=numpy.zeros((Nx,Ny,Nz), dtype=float)
omegaz=numpy.zeros((Nx,Ny,Nz), dtype=float)
omegatot=numpy.zeros((Nx,Ny,Nz), dtype=float)

omegabx=numpy.zeros((Nx,Ny,Nz), dtype=float)
omegaby=numpy.zeros((Nx,Ny,Nz), dtype=float)
omegabz=numpy.zeros((Nx,Ny,Nz), dtype=float)
omegabtot=numpy.zeros((Nx,Ny,Nz), dtype=float)

ux=numpy.zeros((Nx,Ny,Nz), dtype=float)
uy=numpy.zeros((Nx,Ny,Nz), dtype=float)
uz=numpy.zeros((Nx,Ny,Nz), dtype=float)
vx=numpy.zeros((Nx,Ny,Nz), dtype=float)
vy=numpy.zeros((Nx,Ny,Nz), dtype=float)
vz=numpy.zeros((Nx,Ny,Nz), dtype=float)
wx=numpy.zeros((Nx,Ny,Nz), dtype=float)
wy=numpy.zeros((Nx,Ny,Nz), dtype=float)
wz=numpy.zeros((Nx,Ny,Nz), dtype=float)

uxold=numpy.zeros((Nx,Ny,Nz), dtype=float)
uyold=numpy.zeros((Nx,Ny,Nz), dtype=float)
uzold=numpy.zeros((Nx,Ny,Nz), dtype=float)
vxold=numpy.zeros((Nx,Ny,Nz), dtype=float)
vyold=numpy.zeros((Nx,Ny,Nz), dtype=float)
vzold=numpy.zeros((Nx,Ny,Nz), dtype=float)
wxold=numpy.zeros((Nx,Ny,Nz), dtype=float)
wyold=numpy.zeros((Nx,Ny,Nz), dtype=float)
wzold=numpy.zeros((Nx,Ny,Nz), dtype=float)

bxx=numpy.zeros((Nx,Ny,Nz), dtype=float)
bxy=numpy.zeros((Nx,Ny,Nz), dtype=float)
bxz=numpy.zeros((Nx,Ny,Nz), dtype=float)
byx=numpy.zeros((Nx,Ny,Nz), dtype=float)
byy=numpy.zeros((Nx,Ny,Nz), dtype=float)
byz=numpy.zeros((Nx,Ny,Nz), dtype=float)
bzx=numpy.zeros((Nx,Ny,Nz), dtype=float)
bzy=numpy.zeros((Nx,Ny,Nz), dtype=float)
bzz=numpy.zeros((Nx,Ny,Nz), dtype=float)

bxxold=numpy.zeros((Nx,Ny,Nz), dtype=float)
bxyold=numpy.zeros((Nx,Ny,Nz), dtype=float)
bxzold=numpy.zeros((Nx,Ny,Nz), dtype=float)
byxold=numpy.zeros((Nx,Ny,Nz), dtype=float)
byyold=numpy.zeros((Nx,Ny,Nz), dtype=float)
byzold=numpy.zeros((Nx,Ny,Nz), dtype=float)
bzxold=numpy.zeros((Nx,Ny,Nz), dtype=float)
bzyold=numpy.zeros((Nx,Ny,Nz), dtype=float)
bzzold=numpy.zeros((Nx,Ny,Nz), dtype=float)

nonlinu=numpy.zeros((Nx,Ny,Nz), dtype=float)
nonlinv=numpy.zeros((Nx,Ny,Nz), dtype=float)
nonlinw=numpy.zeros((Nx,Ny,Nz), dtype=float)

uhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)
what=numpy.zeros((Nx,Ny,Nz), dtype=complex)
vhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)

bxhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)
byhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)
bzhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)

phat=numpy.zeros((Nx,Ny,Nz), dtype=complex)
temphat=numpy.zeros((Nx,Ny,Nz), dtype=complex)

rhsuhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)
rhsvhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)
rhswhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)

rhsbxhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)
rhsbyhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)
rhsbzhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)

nonlinuhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)
nonlinvhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)
nonlinwhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)

tdata=numpy.zeros((Nt), dtype=float)

# initial conditions for Taylor-Green Vortex
theta=0.0
u=(2.0/(3.0**0.5))*numpy.sin(theta+2.0*math.pi/3.0)*numpy.sin(xx)*numpy.cos(yy)*numpy.cos(zz)
v=(2.0/(3.0**0.5))*numpy.sin(theta-2.0*math.pi/3.0)*numpy.cos(xx)*numpy.sin(yy)*numpy.cos(zz)
w=(2.0/(3.0**0.5))*numpy.sin(theta)*numpy.cos(xx)*numpy.cos(yy)*numpy.sin(zz)

# Exact solution
#sl=1
#sk=1
#sm=1
#lamlkm=(sl**2.0+sk**2.0+sm**2.0)**0.5
#u=-0.5*(lamlkm*sl*numpy.cos(sk*xx)*numpy.sin(sl*yy)*numpy.sin(sm*zz) \
#+sm*sk*numpy.sin(sk*xx)*numpy.cos(sl*yy)*numpy.cos(sm*zz))*numpy.exp(-t*(lamlkm**2.0)/Rey)
#v= 0.5*(lamlkm*sk*numpy.sin(sk*xx)*numpy.cos(sl*yy)*numpy.sin(sm*zz) \
#-sm*sl*numpy.cos(sk*xx)*numpy.sin(sl*yy)*numpy.cos(sm*zz))*numpy.exp(-t*(lamlkm**2.0)/Rey)
#w= numpy.cos(sk*xx)*numpy.cos(sl*yy)*numpy.sin(sm*zz)*numpy.exp(-t*(lamlkm**2.0)/Rey)

# initial fluid field terms
uhat=numpy.fft.fftn(u)
vhat=numpy.fft.fftn(v)
what=numpy.fft.fftn(w)

temphat=kxm*uhat
ux=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*uhat
uy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*uhat
uz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*vhat
vx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*vhat
vy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*vhat
vz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*what
wx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*what
wy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*what
wz=numpy.real(numpy.fft.ifftn(temphat))

# Calculate fluid vorticity for plotting

omegax=wy-vz
omegay=uz-wx
omegaz=vx-uy
omegatot=omegax**2.0 + omegay**2.0 + omegaz**2.0

# initial magnetic field terms
bxhat=numpy.fft.fftn(bx)
byhat=numpy.fft.fftn(by)
bzhat=numpy.fft.fftn(bz)

temphat=kxm*bxhat
bxx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*bxhat
bxy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*bxhat
bxz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*byhat
byx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*byhat
byy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*byhat
byz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*bzhat
bzx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*bzhat
bzy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*bzhat
bzz=numpy.real(numpy.fft.ifftn(temphat))

# Calculate magnetic vorticity for plotting

omegabx=bzy-byz
omegaby=bxz-bzx
omegabz=byx-bxy
omegabtot=omegabx**2.0 + omegaby**2.0 + omegabz**2.0

#src=mlab.contour3d(xx,yy,zz,u,colormap='jet',opacity=0.1,contours=4)
src = mlab.pipeline.scalar_field(xx,yy,zz,omegatot,colormap='YlGnBu')
mlab.pipeline.iso_surface(src, contours=[omegatot.min()+0.1*omegatot.ptp(), ], \
colormap='YlGnBu',opacity=0.85)
mlab.pipeline.iso_surface(src, contours=[omegatot.max()-0.1*omegatot.ptp(), ], \
colormap='YlGnBu',opacity=1.0)

#src = mlab.pipeline.scalar_field(xx,yy,zz,omegabtot,colormap='YlGnBu')
#mlab.pipeline.iso_surface(src, contours=[omegabtot.min()+0.1*omegatot.ptp(), ], \
#   colormap='YlGnBu',opacity=0.85)
#mlab.pipeline.iso_surface(src, contours=[omegabtot.max()-0.1*omegatot.ptp(), ], \
#   colormap='YlGnBu',opacity=1.0)

mlab.pipeline.image_plane_widget(src,plane_orientation='z_axes', \
slice_index=Nz/2,colormap='YlGnBu', \
opacity=0.01)
mlab.pipeline.image_plane_widget(src,plane_orientation='y_axes', \
slice_index=Ny/2,colormap='YlGnBu', \
opacity=0.01)
mlab.pipeline.image_plane_widget(src,plane_orientation='x_axes', \
slice_index=Nx/2,colormap='YlGnBu', \
opacity=0.01)
mlab.scalarbar()
mlab.xlabel('x',object=src)
mlab.ylabel('y',object=src)
mlab.zlabel('z',object=src)

t=0.0
tdata[0]=t
#solve pde and plot results
for n in xrange(Nt):
uold=u
uxold=ux
uyold=uy
uzold=uz
vold=v
vxold=vx
vyold=vy
vzold=vz
wold=w
wxold=wx
wyold=wy
wzold=wz
bxold=bx
bxxold=bxx
bxyold=bxy
bxzold=bxz
byold=by
byxold=byx
byyold=byy
byzold=byz
bzold=bz
bzxold=bzx
bzyold=bzy
bzzold=bzz

rhsuhatfix=(1.0/dt + (0.5/Re)*(k2xm+k2ym+k2zm))*uhat
rhsvhatfix=(1.0/dt + (0.5/Re)*(k2xm+k2ym+k2zm))*vhat
rhswhatfix=(1.0/dt + (0.5/Re)*(k2xm+k2ym+k2zm))*what

rhsbxhatfix=(1.0/dt + (0.5/Rem)*(k2xm+k2ym+k2zm))*bxhat
rhsbyhatfix=(1.0/dt + (0.5/Rem)*(k2xm+k2ym+k2zm))*byhat
rhsbzhatfix=(1.0/dt + (0.5/Rem)*(k2xm+k2ym+k2zm))*bzhat

chg=1.0
t=t+dt
while(chg>tol):
# Fluid field
nonlinu=0.25*((u+uold)*(ux+uxold)+(v+vold)*(uy+uyold)+(w+wold)*(uz+uzold)-\
(bx+bxold)*(bxx+bxxold)-(by+byold)*(bxy+bxyold)-(bz+bzold)*(bxz+bxzold))
nonlinv=0.25*((u+uold)*(vx+vxold)+(v+vold)*(vy+vyold)+(w+wold)*(vz+vzold)-\
(bx+bxold)*(byx+byxold)-(by+byold)*(byy+byyold)-(bz+bzold)*(byz+byzold))
nonlinw=0.25*((u+uold)*(wx+wxold)+(v+vold)*(wy+wyold)+(w+wold)*(wz+wzold)-\
(bx+bxold)*(bzx+bzxold)-(by+byold)*(bzy+bzyold)-(bz+bzold)*(bzz+bzzold))
nonlinuhat=numpy.fft.fftn(nonlinu)
nonlinvhat=numpy.fft.fftn(nonlinv)
nonlinwhat=numpy.fft.fftn(nonlinw)
phat=-1.0*(kxm*nonlinuhat+kym*nonlinvhat+kzm*nonlinwhat)/(k2xm+k2ym+k2zm+0.1**13)
uhat=(rhsuhatfix-nonlinuhat-kxm*phat)/(1.0/dt - (0.5/Re)*(k2xm+k2ym+k2zm))
vhat=(rhsvhatfix-nonlinvhat-kym*phat)/(1.0/dt - (0.5/Re)*(k2xm+k2ym+k2zm))
what=(rhswhatfix-nonlinwhat-kzm*phat)/(1.0/dt - (0.5/Re)*(k2xm+k2ym+k2zm))

temphat=kxm*uhat
ux=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*uhat
uy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*uhat
uz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*vhat
vx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*vhat
vy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*vhat
vz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*what
wx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*what
wy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*what
wz=numpy.real(numpy.fft.ifftn(temphat))

utemp=u
vtemp=v
wtemp=w
u=numpy.real(numpy.fft.ifftn(uhat))
v=numpy.real(numpy.fft.ifftn(vhat))
w=numpy.real(numpy.fft.ifftn(what))

# Magnetic field
nonlinu=0.25*((u+uold)*(bxx+bxxold)+(v+vold)*(bxy+bxyold)+(w+wold)*(bxz+bxzold)-\
(bx+bxold)*(ux+uxold)-(by+byold)*(uy+uyold)-(bz+bzold)*(uz+uzold))
nonlinv=0.25*((u+uold)*(byx+byxold)+(v+vold)*(byy+byyold)+(w+wold)*(byz+byzold)-\
(bx+bxold)*(vx+vxold)-(by+byold)*(vy+vyold)-(bz+bzold)*(vz+vzold))
nonlinw=0.25*((u+uold)*(bzx+bzxold)+(v+vold)*(bzy+bzyold)+(w+wold)*(bzz+bzzold)-\
(bx+bxold)*(wx+wxold)-(by+byold)*(wy+wyold)-(bz+bzold)*(wz+wzold))
nonlinuhat=numpy.fft.fftn(nonlinu)
nonlinvhat=numpy.fft.fftn(nonlinv)
nonlinwhat=numpy.fft.fftn(nonlinw)
phat=-1.0*(kxm*nonlinuhat+kym*nonlinvhat+kzm*nonlinwhat)/(k2xm+k2ym+k2zm+0.1**13)
bxhat=(rhsbxhatfix-nonlinuhat-kxm*phat)/(1.0/dt - (0.5/Rem)*(k2xm+k2ym+k2zm))
byhat=(rhsbyhatfix-nonlinvhat-kym*phat)/(1.0/dt - (0.5/Rem)*(k2xm+k2ym+k2zm))
bzhat=(rhsbzhatfix-nonlinwhat-kzm*phat)/(1.0/dt - (0.5/Rem)*(k2xm+k2ym+k2zm))

temphat=kxm*bxhat
bxx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*bxhat
bxy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*bxhat
bxz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*byhat
byx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*byhat
byy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*byhat
byz=numpy.real(numpy.fft.ifftn(temphat))

temphat=kxm*bzhat
bzx=numpy.real(numpy.fft.ifftn(temphat))
temphat=kym*bzhat
bzy=numpy.real(numpy.fft.ifftn(temphat))
temphat=kzm*bzhat
bzz=numpy.real(numpy.fft.ifftn(temphat))

bxtemp=bx
bytemp=by
bztemp=bz
bx=numpy.real(numpy.fft.ifftn(bxhat))
by=numpy.real(numpy.fft.ifftn(byhat))
bz=numpy.real(numpy.fft.ifftn(bzhat))

chg=numpy.max(abs(u-utemp))+numpy.max(abs(v-vtemp))+numpy.max(abs(w-wtemp))+\
numpy.max(abs(bx-bxtemp))+numpy.max(abs(by-bytemp))+numpy.max(abs(bz-bztemp))
# calculate vorticity for plotting
omegax=wy-vz
omegay=uz-wx
omegaz=vx-uy
omegatot=omegax**2.0 + omegay**2.0 + omegaz**2.0
src.mlab_source.scalars = omegatot
tdata[n]=t
omegabx=bzy-byz
omegaby=bxz-bzx
omegabz=byx-bxy
omegabtot=omegabx**2.0 + omegaby**2.0 + omegabz**2.0
#src.mlab_source.scalars = omegatot


## 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.

( C)

A serial Fortran program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download
	PROGRAM main
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program numerically solves the 2D incompressible magnetohydrodynamics
! equations on a Square Domain [0,1]x[0,1] using pseudo-spectral methods and
! implicit midpoint rule timestepping. A passive scalar is also advected
!
! Periodic free-slip boundary conditions and Initial conditions:
!	u(x,y,0)=sin(2*pi*x)cos(2*pi*y)
!	v(x,y,0)=-cos(2*pi*x)sin(2*pi*y)
!
! .. Parameters ..
!  Nx				= number of modes in x - power of 2 for FFT
!  Ny				= number of modes in y - power of 2 for FFT
!  Nt				= number of timesteps to take
!  Tmax				= maximum simulation time
!  FFTW_IN_PLACE 		= value for FFTW input
!  FFTW_MEASURE 		= value for FFTW input
!  FFTW_EXHAUSTIVE 		= value for FFTW input
!  FFTW_PATIENT 		= value for FFTW input
!  FFTW_ESTIMATE 		= value for FFTW input
!  FFTW_FORWARD    	 	= value for FFTW input
!  FFTW_BACKWARD		= value for FFTW input
!  pi = 3.14159265358979323846264338327950288419716939937510d0
!  mu				= viscosity
!  rho				= density
! .. Scalars ..
!  i				= loop counter in x direction
!  j				= loop counter in y direction
!  n				= loop counter for timesteps direction
!  allocatestatus		= error indicator during allocation
!  count			= keep track of information written to disk
!  iol				= size of array to write to disk
!  start			= variable to record start time of program
!  finish			= variable to record end time of program
!  count_rate			= variable for clock count rate
!  planfx			= Forward 1d fft plan in x
!  planbx			= Backward 1d fft plan in x
!  planfy			= Forward 1d fft plan in y
!  planby			= Backward 1d fft plan in y
!  dt				= timestep
! .. Arrays ..
!  u				= velocity in x direction
!  uold				= velocity in x direction at previous timestep
!  v				= velocity in y direction
!  vold				= velocity in y direction at previous timestep
!  u_y				= y derivative of velocity in x direction
!  v_x				= x derivative of velocity in y direction
!  bx				= x component of magnetic field
!  by				= y component of magnetic field
!  bx				= x component of magnetic field at previous time step
!  by				= y component of magnetic field at previous time step
!  omeg				= vorticity	in real space
!  omegold			= vorticity in real space at previous
!						iterate
!  omegcheck			= store of vorticity at previous iterate
!  omegoldhat			= 2D Fourier transform of vorticity at previous
!						iterate
!  omegoldhat_x			= x-derivative of vorticity in Fourier space
!						at previous iterate
!  omegold_x			= x-derivative of vorticity in real space
!						at previous iterate
!  omegoldhat_y 		= y-derivative of vorticity in Fourier space
!						at previous iterate
!  omegold_y			= y-derivative of vorticity in real space
!						 at previous iterate
!  nlold			= nonlinear term in real space
!						at previous iterate
!  nloldhat			= nonlinear term in Fourier space
!						at previous iterate
!  omeghat			= 2D Fourier transform of vorticity
!						at next iterate
!  omeghat_x			= x-derivative of vorticity in Fourier space
!						at next timestep
!  omeghat_y			= y-derivative of vorticity in Fourier space
!						at next timestep
!  omeg_x			= x-derivative of vorticity in real space
!						at next timestep
!  omeg_y			= y-derivative of vorticity in real space
!						at next timestep
!  alpha			= magnetic vorticity in real space
!  alphaold			= magnetic vorticity in real space at previous
!						iterate
!  alphacheck			= store of vorticity potential at previous iterate
!  alphaoldhat			= 2D Fourier transform of vorticity potential at previous
!					iterate
!  alphaoldhat_x		= x-derivative of magnetic vorticity in Fourier space
!						at previous iterate
!  alphaold_x			= x-derivative of magnetic vorticity in real space
!						at previous iterate
!  alphaoldhat_y 		= y-derivative of magnetic vorticity in Fourier space
!						at previous iterate
!  alphaold_y			= y-derivative of  magnetic vorticity in real space
!  alphahat			= 2D Fourier transform of magnetic vorticity
!						at next iterate
!  alphahat_x			= x-derivative of magnetic vorticity in Fourier space
!						at next timestep
!  alphahat_y			= y-derivative of magnetic vorticity in Fourier space
!						at next timestep
!  alpha_x			= x-derivative of magnetic vorticity in real space
!						at next timestep
!  alpha_y			= y-derivative of magnetic vorticity in real space
!						at next timestep
!  psihat			= streamfunction at next timestep
!  psi_x			= x-derivative of streamfunction in real space
!						at next timestep
!  psi_y			= y-derivative of streamfunction in real space
!						at next timestep
!  phi				= Magnetic potential at next timestep
!  phihat			= 2D Fourier transform of magnetic potential at next timestep
!  phi_x			= x-derivative of streamfunction in real space
!						at next timestep
!  phi_y			= y-derivative of streamfunction in real space
!						at next timestep
! .. Vectors ..
!  kx				= fourier frequencies in x direction
!  ky				= fourier frequencies in y direction
!  kxx				= square of fourier frequencies in x direction
!  kyy				= square of fourier frequencies in y direction
!  x				= x locations
!  y				= y locations
!  time				= times at which save data
!  name_config			= array to store filename for data to be saved
!  fftfx			= array to setup x Fourier transform
!  fftbx			= array to setup y Fourier transform
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! This program has not been optimized to use the least amount of memory
! but is intended as an example only for which all states can be saved
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! FFTW3	 -- Fast Fourier Transform in the West Library
!			(http://www.fftw.org/)
! declare variables

IMPLICIT NONE
INTEGER(kind=4), PARAMETER 	::  Nx=64
INTEGER(kind=4), PARAMETER 	::  Ny=64
REAL(kind=8), PARAMETER		::  dt=0.00125
REAL(kind=8), PARAMETER	&
::  pi=3.14159265358979323846264338327950288419716939937510
REAL(kind=8), PARAMETER		::  Re=1.0d0
REAL(kind=8), PARAMETER		::  Rem=1.0d0
REAL(kind=8), PARAMETER		::	tol=0.1d0**10
REAL(kind=8)				::	chg
INTEGER(kind=4), PARAMETER	::	nplots=50
REAL(kind=8), DIMENSION(:), ALLOCATABLE			::	time
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE		::  kx,kxx
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE		::  ky,kyy
REAL(kind=8), DIMENSION(:), ALLOCATABLE			::  x
REAL(kind=8), DIMENSION(:), ALLOCATABLE			::  y
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE		:: &
u,uold,v,vold,u_y,v_x,omegold, omegcheck, omeg,&
omegoldhat, omegold_x, omegold_y, &
omeghat,  omeg_x, omeg_y,&
nl, nlhat, psihat,  psi_x, psi_y,&
theta, thetaold, thetax, thetay, thetacheck, thetahat, thetaoldhat, &
alpha, alphaold, alpha_x, alpha_y, bx, by, bxold, byold, alphahat, alphaoldhat, alphacheck, &
phi, phihat,  phi_x, phi_y, temp
INTEGER(kind=4)						::  i,j,k,n, allocatestatus, count, iol
INTEGER(kind=4)			 			::	start, finish, count_rate
INTEGER(kind=4), PARAMETER      			::  FFTW_IN_PLACE = 8, FFTW_MEASURE = 0, &
FFTW_EXHAUSTIVE = 8, FFTW_PATIENT = 32,    &
FFTW_ESTIMATE = 64
INTEGER(kind=4),PARAMETER       			::  FFTW_FORWARD = -1, FFTW_BACKWARD=1
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE		::  fftfx,fftbx
INTEGER(kind=8)						::  planfxy,planbxy
CHARACTER*100			 			::  name_config

CALL system_clock(start,count_rate)
ALLOCATE(time(1:nplots),kx(1:Nx),kxx(1:Nx),ky(1:Ny),kyy(1:Ny),x(1:Nx),y(1:Ny),&
u(1:Nx,1:Ny),uold(1:Nx,1:Ny),v(1:Nx,1:Ny),vold(1:Nx,1:Ny),u_y(1:Nx,1:Ny),&
v_x(1:Nx,1:Ny),omegold(1:Nx,1:Ny),omegcheck(1:Nx,1:Ny), omeg(1:Nx,1:Ny),&
omegoldhat(1:Nx,1:Ny), omegold_x(1:Nx,1:Ny), omegold_y(1:Nx,1:Ny), &
omeghat(1:Nx,1:Ny),  omeg_x(1:Nx,1:Ny), omeg_y(1:Nx,1:Ny),&
nl(1:Nx,1:Ny), nlhat(1:Nx,1:Ny), psihat(1:Nx,1:Ny), &
psi_x(1:Nx,1:Ny),  psi_y(1:Nx,1:Ny),&
theta(1:Nx,1:Ny), thetaold(1:Nx,1:Ny), thetax(1:Nx,1:Ny), thetay(1:Nx,1:Ny), &
thetacheck(1:Nx,1:Ny), thetahat(1:Nx,1:Ny), thetaoldhat(1:Nx,1:Ny), &
alpha(1:Nx,1:Ny), alphaold(1:Nx,1:Ny), alpha_x(1:Nx,1:Ny), alpha_y(1:Nx,1:Ny),&
bx(1:Nx,1:Ny), by(1:Nx,1:Ny), bxold(1:Nx,1:Ny), byold(1:Nx,1:Ny),&
alphahat(1:Nx,1:Ny), alphaoldhat(1:Nx,1:Ny),alphacheck(1:Nx,1:Ny), &
phi(1:Nx,1:Ny), phihat(1:Nx,1:Ny),  &
phi_x(1:Nx,1:Ny), phi_y(1:Nx,1:Ny), temp(1:Nx,1:Ny),&
fftfx(1:Nx,1:Ny),fftbx(1:Nx,1:Ny),stat=AllocateStatus)
IF (AllocateStatus .ne. 0) STOP
PRINT *,'allocated space'

! set up ffts
CALL dfftw_plan_dft_2d_(planfxy,Nx,Ny,fftfx(1:Nx,1:Ny),fftbx(1:Nx,1:Ny),&
FFTW_FORWARD,FFTW_EXHAUSTIVE)
CALL dfftw_plan_dft_2d_(planbxy,Nx,Ny,fftbx(1:Nx,1:Ny),fftfx(1:Nx,1:Ny),&
FFTW_BACKWARD,FFTW_EXHAUSTIVE)

! setup fourier frequencies in x-direction
DO i=1,1+Nx/2
kx(i)= 2.0d0*pi*cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))
END DO
kx(1+Nx/2)=0.0d0
DO i = 1,Nx/2 -1
kx(i+1+Nx/2)=-kx(1-i+Nx/2)
END DO
DO i=1,Nx
kxx(i)=kx(i)*kx(i)
END DO
DO i=1,Nx
x(i)=REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0))
END DO

! setup fourier frequencies in y-direction
DO j=1,1+Ny/2
ky(j)= 2.0d0*pi*cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0))
END DO
ky(1+Ny/2)=0.0d0
DO j = 1,Ny/2 -1
ky(j+1+Ny/2)=-ky(1-j+Ny/2)
END DO
DO j=1,Ny
kyy(j)=ky(j)*ky(j)
END DO
DO j=1,Ny
y(j)=REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0))
END DO
PRINT *,'Setup grid and fourier frequencies'

DO j=1,Ny
DO i=1,Nx
u(i,j)=sin(2.0d0*pi*x(i))*cos(2.0d0*pi*y(j))
v(i,j)=-cos(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j))
u_y(i,j)=-2.0d0*pi*sin(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j))
v_x(i,j)=2.0d0*pi*sin(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j))
omeg(i,j)=v_x(i,j)-u_y(i,j)
theta(i,j)=exp(-2.0*((4.0*pi*(x(i)-0.3))**2+(4.0*pi*(y(j)-0.3))**2))
alpha(i,j)=sin(4.0*pi*x(i))*cos(6.0*pi*y(j))
END DO
END DO

! Vorticity to Fourier Space
CALL dfftw_execute_dft_(planfxy,omeg(1:Nx,1:Ny),omeghat(1:Nx,1:Ny))

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
time(1)=0.0d0
PRINT *,'Got initial data, starting timestepping'
DO n=1,nplots
chg=1
! save old values
uold(1:Nx,1:Ny)=u(1:Nx,1:Ny)
vold(1:Nx,1:Ny)=v(1:Nx,1:Ny)
omegold(1:Nx,1:Ny)=omeg(1:Nx,1:Ny)
omegcheck(1:Nx,1:Ny)=omeg(1:Nx,1:Ny)
omegoldhat(1:Nx,1:Ny)=omeghat(1:Nx,1:Ny)
thetaold(1:Nx,1:Ny)=theta(1:Nx,1:Ny)
thetaoldhat(1:Nx,1:Ny)=thetahat(1:Nx,1:Ny)
thetacheck(1:Nx,1:Ny)=theta(1:Nx,1:Ny)
bxold(1:Nx,1:Ny)=bx(1:Nx,1:Ny)
byold(1:Nx,1:Ny)=by(1:Nx,1:Ny)
alphaoldhat(1:Nx,1:Ny)=alphahat(1:Nx,1:Ny)
alphaold(1:Nx,1:Ny)=alpha(1:Nx,1:Ny)
alphacheck(1:Nx,1:Ny)=alpha(1:Nx,1:Ny)
DO WHILE (chg>tol)
! Fluid flow field
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!nonlinear fixed (n,k+1)!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! obtain \hat{\omega}_x^{n+1,k}
DO i=1,Nx; DO j=1,Ny;
temp(i,j)=0.5d0*(omeghat(i,j)+omegoldhat(i,j))*kx(i)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
! convert back to real space
CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),omeg_x(1:Nx,1:Ny))
! obtain \hat{\omega}_y^{n+1,k}
DO i=1,Nx; DO j=1,Ny;
temp(i,j)=0.5d0*(omeghat(i,j)+omegoldhat(i,j))*ky(j)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
! convert back to real space
CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),omeg_y(1:Nx,1:Ny))
! calculate nonlinear term in real space
DO i=1,Nx; DO j=1,Ny;
nl(i,j)=0.5d0*( (u(i,j) + uold(i,j))*omeg_x(i,j)+&
(v(i,j) + vold(i,j))*omeg_y(i,j)-&
(bx(i,j) + bxold(i,j))*alpha_x(i,j)-&
(by(i,j) + byold(i,j))*alpha_y(i,j) )
END DO; END DO
! convert back to fourier
CALL dfftw_execute_dft_(planfxy,nl(1:Nx,1:Ny),nlhat(1:Nx,1:Ny))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! obtain \hat{\omega}^{n+1,k+1} with implicit midpoint rule timestepping
DO i=1,Nx; DO j=1,Ny;
omeghat(i,j)=( (1.0d0/dt+0.5d0*(1.0d0/Re)*(kxx(i)+kyy(j)))&
*omegoldhat(i,j) - nlhat(i,j))/&
(1.0d0/dt-0.5d0*(1.0d0/Re)*(kxx(i)+kyy(j)))
END DO; END DO

! calculate \hat{\psi}^{n+1,k+1}
DO i=1,Nx; DO j=1,Ny;
psihat(i,j)=-omeghat(i,j)/(kxx(i)+kyy(j))
END DO; END DO
psihat(1,1)=0.0d0
psihat(Nx/2+1,Ny/2+1)=0.0d0
psihat(Nx/2+1,1)=0.0d0
psihat(1,Ny/2+1)=0.0d0

! obtain \psi_x^{n+1,k+1} and \psi_y^{n+1,k+1}
DO i=1,Nx; DO j=1,Ny;
temp(i,j)=psihat(i,j)*kx(i)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),psi_x(1:Nx,1:Ny))
DO i=1,Nx; DO j=1,Ny
temp(i,j)=psihat(i,j)*ky(j)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
CALL dfftw_execute_dft_(planbxy,temp(1:Ny,1:Ny),psi_y(1:Ny,1:Ny))

! obtain \omega^{n+1,k+1}
CALL dfftw_execute_dft_(planbxy,omeghat(1:Nx,1:Ny),omeg(1:Nx,1:Ny))
DO i=1,Nx; DO j=1,Ny
omeg(i,j)=omeg(i,j)/REAL(Nx*Ny,kind(0d0))
END DO; END DO

! obtain u^{n+1,k+1} and v^{n+1,k+1} using stream function (\psi) in real space
DO i=1,Nx; DO j=1,Ny
u(i,j)=psi_y(i,j)
v(i,j)=-psi_x(i,j)
END DO; END DO

! Magnetic field
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!nonlinear fixed (n,k+1)!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! obtain \hat{\omega}_x^{n+1,k}
DO i=1,Nx; DO j=1,Ny;
temp(i,j)=0.5*(alphahat(i,j)+alphaoldhat(i,j))*kx(i)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
! convert back to real space
CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),alpha_x(1:Nx,1:Ny))
! obtain \hat{\omega}_y^{n+1,k}
DO i=1,Nx; DO j=1,Ny;
temp(i,j)=0.5*(alphahat(i,j)+alphaoldhat(i,j))*ky(j)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
! convert back to real space
CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),omeg_y(1:Nx,1:Ny))
! calculate nonlinear term in real space
DO i=1,Nx; DO j=1,Ny
nl(i,j)=0.5d0*( (u(i,j) + uold(i,j))*alpha_x(i,j)+&
(v(i,j) + vold(i,j))*alpha_y(i,j)-&
(bx(i,j) + bxold(i,j))*omeg_x(i,j)-&
(by(i,j) + byold(i,j))*omeg_y(i,j) )
END DO; END DO
! convert back to fourier
CALL dfftw_execute_dft_(planfxy,nl(1:Nx,1:Ny),nlhat(1:Nx,1:Ny))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! obtain \hat{\omega}^{n+1,k+1} with implicit midpoint rule timestepping
DO i=1,Nx; DO j=1,Ny
alphahat(i,j)=( (1.0d0/dt+0.5d0*(1.0d0/Rem)*(kxx(i)+kyy(j)))&
*alphaoldhat(i,j) - nlhat(i,j))/&
(1.0d0/dt-0.5d0*(1.0d0/Rem)*(kxx(i)+kyy(j)))
END DO; END DO

! calculate \hat{\psi}^{n+1,k+1}
DO i=1,Nx; DO j=1,Ny
phihat(i,j)=-alphahat(i,j)/(kxx(i)+kyy(j))
END DO; END DO
phihat(1,1)=0.0d0
phihat(Nx/2+1,Ny/2+1)=0.0d0
phihat(Nx/2+1,1)=0.0d0
phihat(1,Ny/2+1)=0.0d0

! obtain \psi_x^{n+1,k+1} and \psi_y^{n+1,k+1}
DO i=1,Nx; DO j=1,Ny;
temp(i,j)=phihat(i,j)*kx(i)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
CALL dfftw_execute_dft_(planbxy,temp(1:Nx,1:Ny),phi_x(1:Nx,1:Ny))
DO i=1,Nx; DO j=1,Ny
temp(i,j)=phihat(i,j)*ky(j)/REAL(Nx*Ny,kind(0d0))
END DO; END DO
CALL dfftw_execute_dft_(planbxy,temp(1:Ny,1:Ny),phi_y(1:Ny,1:Ny))

! obtain \omega^{n+1,k+1}
CALL dfftw_execute_dft_(planbxy,alphahat(1:Nx,1:Ny),alpha(1:Nx,1:Ny))
DO i=1,Nx; DO j=1,Ny
alpha(i,j)=alpha(i,j)/REAL(Nx*Ny,kind(0d0))
END DO; END DO

! obtain u^{n+1,k+1} and v^{n+1,k+1} using stream function (\psi) in real space
DO i=1,Nx; DO j=1,Ny
bx(i,j)=phi_y(i,j)
END DO; END DO
DO i=1,Nx; DO j=1,Ny
by(i,j)=-phi_x(i,j)
END DO; END DO

! check for convergence
chg=maxval(abs(omeg-omegcheck)) + maxval(abs(alpha-alphacheck))
! saves {n+1,k+1} to {n,k} for next iteration
DO i=1,Nx; DO j=1,Ny
omegcheck(i,j)=omeg(i,j)
END DO; END DO
DO i=1,Nx; DO j=1,Ny
alphacheck(i,j)=alpha(i,j)
END DO; END DO
END DO
time(n+1)=time(n)+dt
PRINT *,'TIME ',time(n+1)
END DO

name_config = 'omegafinal.datbin'
INQUIRE(iolength=iol) omeg(1,1)
OPEN(unit=11,FILE=name_config,form="unformatted", access="direct",recl=iol)
count = 1
DO j=1,Ny
DO i=1,Nx
WRITE(11,rec=count) REAL(omeg(i,j),KIND(0d0))
count=count+1
END DO
END DO
CLOSE(11)

name_config = 'alphfinal.datbin'
INQUIRE(iolength=iol) alpha(1,1)
OPEN(unit=11,FILE=name_config,form="unformatted", access="direct",recl=iol)
count = 1
DO j=1,Ny
DO i=1,Nx
WRITE(11,rec=count) REAL(alpha(i,j),KIND(0d0))
count=count+1
END DO
END DO
CLOSE(11)

name_config = 'xcoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO i=1,Nx
WRITE(11,*) x(i)
END DO
CLOSE(11)

name_config = 'ycoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Ny
WRITE(11,*) y(j)
END DO
CLOSE(11)

CALL dfftw_destroy_plan_(planfxy)
CALL dfftw_destroy_plan_(planbxy)
CALL dfftw_cleanup_()

DEALLOCATE(time,kx,kxx,ky,kyy,x,y,&
u,uold,v,vold,u_y,&
v_x,omegold,omegcheck, omeg,&
omegoldhat, omegold_x, &
omegold_y, &
omeghat,  omeg_x,&
omeg_y, nl, nlhat, psihat, &
psi_x,  psi_y,&
theta, thetaold, thetax, thetay, &
thetacheck, thetahat, thetaoldhat, &
alpha, alphaold, alpha_x, alpha_y,&
bx, by,&
bxold, byold, alphahat, alphaoldhat, &
alphacheck, phi, phihat,  &
phi_x, phi_y, temp,&
fftfx,fftbx,stat=AllocateStatus)
IF (AllocateStatus .ne. 0) STOP
PRINT *,'Program execution complete'
END PROGRAM main


( D)

A makefile to compile the program which finds a numerical solution to the 2D magnetohydrodynamic equations Code download
COMPILER =  gfortran
FLAGS = -O0

LIBS = -L${FFTW_LINK} -lfftw3 -lm SOURCES = magnetohydrodynamics2D.f90 mhd2d:$(SOURCES)
${COMPILER} -o mhd2d$(FLAGS) $(SOURCES)$(LIBS)

clean:
rm -f *.o
rm -f *.mod
clobber:
rm -f mhd2d


## 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.

( E)

A parallel Fortran program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download
	PROGRAM main
!-----------------------------------------------------------------------------------
!
!
! PURPOSE
!
! This program numerically solves the 3D incompressible magnetohydrodynamic equations
! on a Cubic Domain [0,2pi]x[0,2pi]x[0,2pi] using fourier pseudo-spectral methods and
! Implicit Midpoint rule timestepping.
!
! Analytical Solution:
!	u(x,y,z,t)=-0.25*(cos(x)sin(y)sin(z)+sin(x)cos(y)cos(z))exp(-t/Re)
!	v(x,y,z,t)= 0.25*(sin(x)cos(y)sin(z)-cos(x)sin(y)cos(z))exp(-t/Re)
!	w(x,y,z,t)= 0.5*cos(x)cos(y)sin(z)exp(-t/Re)
!
! .. Parameters ..
!  Nx				= number of modes in x - power of 2 for FFT
!  Ny				= number of modes in y - power of 2 for FFT
!  Nz				= number of modes in z - power of 2 for FFT
!  Nt				= number of timesteps to take
!  Tmax				= maximum simulation time
!  FFTW_IN_PLACE 	= value for FFTW input
!  FFTW_MEASURE 	= value for FFTW input
!  FFTW_EXHAUSTIVE 	= value for FFTW input
!  FFTW_PATIENT 	= value for FFTW input
!  FFTW_ESTIMATE 	= value for FFTW input
!  FFTW_FORWARD     = value for FFTW input
!  FFTW_BACKWARD	= value for FFTW input
!  pi = 3.14159265358979323846264338327950288419716939937510d0
!  Re				= Reynolds number
! .. Scalars ..
!  i				= loop counter in x direction
!  j				= loop counter in y direction
!  k				= loop counter in z direction
!  n				= loop counter for timesteps direction
!  allocatestatus	= error indicator during allocation
!  count			= keep track of information written to disk
!  iol				= size of array to write to disk
!  start			= variable to record start time of program
!  finish			= variable to record end time of program
!  count_rate		= variable for clock count rate
!  planfxyz			= Forward 3d fft plan
!  planbxyz			= Backward 3d fft plan
!  dt				= timestep
! .. Arrays ..
!  u				= velocity in x direction
!  v				= velocity in y direction
!  w				= velocity in z direction
!  uold				= velocity in x direction at previous timestep
!  vold				= velocity in y direction at previous timestep
!  wold				= velocity in z direction at previous timestep
!  ux				= x derivative of velocity in x direction
!  uy				= y derivative of velocity in x direction
!  uz				= z derivative of velocity in x direction
!  vx				= x derivative of velocity in y direction
!  vy				= y derivative of velocity in y direction
!  vz				= z derivative of velocity in y direction
!  wx				= x derivative of velocity in z direction
!  wy				= y derivative of velocity in z direction
!  wz				= z derivative of velocity in z direction
!  uxold			= x derivative of velocity in x direction
!  uyold			= y derivative of velocity in x direction
!  uzold			= z derivative of velocity in x direction
!  vxold			= x derivative of velocity in y direction
!  vyold			= y derivative of velocity in y direction
!  vzold			= z derivative of velocity in y direction
!  wxold			= x derivative of velocity in z direction
!  wyold			= y derivative of velocity in z direction
!  wzold			= z derivative of velocity in z direction
!  utemp			= temporary storage of u to check convergence
!  vtemp			= temporary storage of v to check convergence
!  wtemp			= temporary storage of w to check convergence
!  temp_r			= temporary storage for untransformed variables
!  uhat				= Fourier transform of u
!  vhat				= Fourier transform of v
!  what				= Fourier transform of w
!  rhsuhatfix			= Fourier transform of righthand side for u for timestepping
!  rhsvhatfix			= Fourier transform of righthand side for v for timestepping
!  rhswhatfix			= Fourier transform of righthand side for w for timestepping
!  bx				= x component of magnetic field
!  by				= y component of magnetic field
!  bz				= z component of magnetic field
!  bxold			= x component of magnetic field at previous timestep
!  byold			= y component of magnetic field at previous timestep
!  bzold			= z component of magnetic field at previous timestep
!  bxx				= x derivative of x component of magnetic field
!  bxy				= y derivative of x component of magnetic field
!  bxz				= z derivative of x component of magnetic field
!  byx				= x derivative of y component of magnetic field
!  byy				= y derivative of y component of magnetic field
!  byz				= z derivative of y component of magnetic field
!  bzx				= x derivative of z component of magnetic field
!  bzy				= y derivative of z component of magnetic field
!  bzz				= z derivative of z component of magnetic field
!  bxxold			= x derivative of x component of magnetic field
!  bxyold			= y derivative of x component of magnetic field
!  bxzold			= z derivative of x component of magnetic field
!  byxold			= x derivative of y component of magnetic field
!  byyold			= y derivative of y component of magnetic field
!  byzold			= z derivative of y component of magnetic field
!  bzxold			= x derivative of z component of magnetic field
!  bzyold			= y derivative of z component of magnetic field
!  bzzold			= z derivative of z component of magnetic field
!  bxtemp			= temporary storage of bx to check convergence
!  bytemp			= temporary storage of by to check convergence
!  bztemp			= temporary storage of bz to check convergence
!  bxhat			= Fourier transform of bx
!  byhat			= Fourier transform of by
!  bzhat			= Fourier transform of bz
!  rhsbxhatfix			= Fourier transform of righthand side for bx for timestepping
!  rhsbyhatfix			= Fourier transform of righthand side for by for timestepping
!  rhsbzhatfix			= Fourier transform of righthand side for bz for timestepping
!  nonlinuhat			= Fourier transform of nonlinear term for u
!  nonlinvhat  			= Fourier transform of nonlinear term for u
!  nonlinwhat			= Fourier transform of nonlinear term for u
!  phat				= Fourier transform of nonlinear term for pressure, p
!  temp_c			= temporary storage for Fourier transforms
!  realtemp			= Real storage
!
! .. Vectors ..
!  kx				= fourier frequencies in x direction
!  ky				= fourier frequencies in y direction
!  kz				= fourier frequencies in z direction
!  x				= x locations
!  y				= y locations
!  z				= y locations
!  time				= times at which save data
!  name_config			= array to store filename for data to be saved
!
! REFERENCES
!
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
!
! This program has not been optimized to use the least amount of memory
! but is intended as an example only for which all states can be saved
!
!--------------------------------------------------------------------------------
! External routines required
!
! External libraries required
! 2DECOMP&FFT -- Fast Fourier Transform in the West Library
!			(http://2decomp.org/)

USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
USE MPI
IMPLICIT NONE
! declare variables
INTEGER(kind=4), PARAMETER 		:: Nx=64
INTEGER(kind=4), PARAMETER 		:: Ny=64
INTEGER(kind=4), PARAMETER 		:: Nz=64
INTEGER(kind=4), PARAMETER 		:: Lx=1
INTEGER(kind=4), PARAMETER 		:: Ly=1
INTEGER(kind=4), PARAMETER 		:: Lz=1
INTEGER(kind=4), PARAMETER		:: Nt=5
REAL(kind=8), PARAMETER			:: dt=0.05d0/Nt
REAL(kind=8), PARAMETER			:: Re=1.0d0
REAL(kind=8), PARAMETER			:: Rem=1.0d0
REAL(kind=8), PARAMETER			:: tol=0.1d0**10
REAL(kind=8), PARAMETER			:: theta=0.0d0

REAL(kind=8), PARAMETER	&
::  pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8), PARAMETER		::	ReInv=1.0d0/REAL(Re,kind(0d0))
REAL(kind=8), PARAMETER		::	RemInv=1.0d0/REAL(Rem,kind(0d0))
REAL(kind=8), PARAMETER		::  	dtInv=1.0d0/REAL(dt,kind(0d0))
REAL(kind=8)			:: 	scalemodes,chg,factor
REAL(kind=8), DIMENSION(:), ALLOCATABLE		:: x, y, z, time,mychg,allchg
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE	:: u, v, w,&
ux, uy, uz,&
vx, vy, vz,&
wx, wy, wz,&
uold, uxold, uyold, uzold,&
vold, vxold, vyold, vzold,&
wold, wxold, wyold, wzold,&
utemp, vtemp, wtemp, temp_r, &
bx, by, bz,&
bxx, bxy, bxz,&
byx, byy, byz,&
bzx, bzy, bzz,&
bxold, bxxold, bxyold, bxzold,&
byold, byxold, byyold, byzold,&
bzold, bzxold, bzyold, bzzold,&
bxtemp, bytemp, bztemp

COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx, ky, kz
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE	:: uhat, vhat, what, &
bxhat, byhat, bzhat, &
rhsuhatfix, rhsvhatfix, &
rhswhatfix, rhsbxhatfix, &
rhsbyhatfix, rhsbzhatfix, &
nonlinuhat, nonlinvhat, &
nonlinwhat, phat, temp_c
REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE 	::  realtemp
! MPI and 2DECOMP variables
TYPE(DECOMP_INFO)				::  decomp
INTEGER(kind=MPI_OFFSET_KIND) 			::  filesize, disp
INTEGER(kind=4)					::  p_row=0, p_col=0, numprocs, myid, ierr

! variables used for saving data and timing
INTEGER(kind=4)					:: count, iol
INTEGER(kind=4)					:: i,j,k,n,t,allocatestatus
INTEGER(kind=4)					:: ind, numberfile
CHARACTER*100			 		:: name_config
INTEGER(kind=4)					:: start, finish, count_rate

! initialisation of 2DECOMP&FFT
CALL MPI_INIT(ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
! do automatic domain decomposition
CALL decomp_2d_init(Nx,Ny,Nz,p_row,p_col)
! get information about domain decomposition choosen
CALL decomp_info_init(Nx,Ny,Nz,decomp)
! initialise FFT library
CALL decomp_2d_fft_init
IF (myid.eq.0) THEN
PRINT *,'Grid:',Nx,'X',Ny,'Y',Nz,'Z'
PRINT *,'dt:',dt
END IF
ALLOCATE(x(1:Nx),y(1:Ny),z(1:Nz),time(1:Nt+1),mychg(1:6),allchg(1:6),&
u(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
v(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
w(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
ux(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
uy(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
uz(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vx(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vy(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vz(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wx(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wy(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wz(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
uold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
uxold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
uyold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
uzold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vxold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vyold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vzold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wxold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wyold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wzold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
utemp(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
vtemp(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
wtemp(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
temp_r(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bx(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
by(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bz(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxx(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxy(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxz(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
byx(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
byy(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
byz(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bzx(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bzy(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bzz(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxxold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxyold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxzold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
byold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
byxold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
byyold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
byzold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bzold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bzxold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bzyold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bzzold(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bxtemp(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bytemp(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
bztemp(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)),&
kx(1:Nx),ky(1:Ny),kz(1:Nz),&
uhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
vhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
what(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
bxhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
byhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
bzhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
rhsuhatfix(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
rhsvhatfix(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
rhswhatfix(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
rhsbxhatfix(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
rhsbyhatfix(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
rhsbzhatfix(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
nonlinuhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
nonlinvhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
nonlinwhat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
phat(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
temp_c(decomp%zst(1):decomp%zen(1),&
decomp%zst(2):decomp%zen(2),&
decomp%zst(3):decomp%zen(3)),&
realtemp(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)), stat=AllocateStatus)
IF (AllocateStatus .ne. 0) STOP
IF (myid.eq.0) THEN
PRINT *,'allocated space'
END IF

! setup fourier frequencies in x-direction
DO i=1,Nx/2+1
kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/Lx
END DO
kx(1+Nx/2)=0.0d0
DO i = 1,Nx/2 -1
kx(i+1+Nx/2)=-kx(1-i+Nx/2)
END DO
ind=1
DO i=-Nx/2,Nx/2-1
x(ind)=2.0d0*pi*REAL(i,kind(0d0))*Lx/REAL(Nx,kind(0d0))
ind=ind+1
END DO
! setup fourier frequencies in y-direction
DO j=1,Ny/2+1
ky(j)= cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0))/Ly
END DO
ky(1+Ny/2)=0.0d0
DO j = 1,Ny/2 -1
ky(j+1+Ny/2)=-ky(1-j+Ny/2)
END DO
ind=1
DO j=-Ny/2,Ny/2-1
y(ind)=2.0d0*pi*REAL(j,kind(0d0))*Ly/REAL(Ny,kind(0d0))
ind=ind+1
END DO
! setup fourier frequencies in z-direction
DO k=1,Nz/2+1
kz(k)= cmplx(0.0d0,1.0d0)*REAL(k-1,kind(0d0))/Lz
END DO
kz(1+Nz/2)=0.0d0
DO k = 1,Nz/2 -1
kz(k+1+Nz/2)=-kz(1-k+Nz/2)
END DO
ind=1
DO k=-Nz/2,Nz/2-1
z(ind)=2.0d0*pi*REAL(k,kind(0d0))*Lz/REAL(Nz,kind(0d0))
ind=ind+1
END DO
scalemodes=1.0d0/REAL(Nx*Ny*Nz,kind(0d0))
IF (myid.eq.0) THEN
PRINT *,'Setup grid and fourier frequencies'
END IF

! Initial conditions for fluid flow field
!initial conditions for Taylor-Green vortex
!	factor=2.0d0/sqrt(3.0d0)
!	DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
!		u(i,j,k)=factor*sin(theta+2.0d0*pi/3.0d0)*sin(x(i))*cos(y(j))*cos(z(k))
!	END DO; END DO; END DO
!	DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
!		v(i,j,k)=factor*sin(theta-2.0d0*pi/3.0d0)*cos(x(i))*sin(y(j))*cos(z(k))
!	END DO ; END DO ; END DO
!	DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
!		w(i,j,k)=factor*sin(theta)*cos(x(i))*cos(y(j))*sin(z(k))
!	END DO ; END DO ; END DO

time(1)=0.0d0
factor=sqrt(3.0d0)
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
u(i,j,k)=-0.5*( factor*cos(x(i))*sin(y(j))*sin(z(k))&
+sin(x(i))*cos(y(j))*cos(z(k)) )*exp(-(factor**2)*time(1)/Re)
END DO; END DO; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
v(i,j,k)=0.5*(  factor*sin(x(i))*cos(y(j))*sin(z(k))&
-cos(x(i))*sin(y(j))*cos(z(k)) )*exp(-(factor**2)*time(1)/Re)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
w(i,j,k)=cos(x(i))*cos(y(j))*sin(z(k))*exp(-(factor**2)*time(1)/Re)
END DO ; END DO ; END DO

CALL decomp_2d_fft_3d(u,uhat,DECOMP_2D_FFT_FORWARD)
CALL decomp_2d_fft_3d(v,vhat,DECOMP_2D_FFT_FORWARD)
CALL decomp_2d_fft_3d(w,what,DECOMP_2D_FFT_FORWARD)

! derivative of u with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=uhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,ux,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=uhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,uy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,uz,DECOMP_2D_FFT_BACKWARD)

! derivative of v with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=vhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,vx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=vhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,vy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=vhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,vz,DECOMP_2D_FFT_BACKWARD)

! derivative of w with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=what(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=what(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=what(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wz,DECOMP_2D_FFT_BACKWARD)
! save initial data
n=0
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(wy(i,j,k)-vz(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegax'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegay
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(uz(i,j,k)-wx(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegay'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegaz
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(vx(i,j,k)-uy(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegaz'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)

! Initial conditions for magnetic field
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bx(i,j,k)=-0.5*( sin(y(j))*sin(z(k)))
END DO; END DO; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
by(i,j,k)=0.5*(  sin(x(i))*sin(z(k)))
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bz(i,j,k)=cos(x(i))*cos(y(j))
END DO ; END DO ; END DO

CALL decomp_2d_fft_3d(bx,bxhat,DECOMP_2D_FFT_FORWARD)
CALL decomp_2d_fft_3d(by,byhat,DECOMP_2D_FFT_FORWARD)
CALL decomp_2d_fft_3d(bz,bzhat,DECOMP_2D_FFT_FORWARD)

! derivative of u with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bxhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bxx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bxhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bxy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bxhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bxz,DECOMP_2D_FFT_BACKWARD)

! derivative of v with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=byhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,byx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=byhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,byy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=byhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,byz,DECOMP_2D_FFT_BACKWARD)

! derivative of w with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bzhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bzhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bzhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bzz,DECOMP_2D_FFT_BACKWARD)
! save initial data
n=0
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(bzy(i,j,k)-byz(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegabx'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegay
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(bxz(i,j,k)-bzx(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegaby'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegaz
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(byx(i,j,k)-bxy(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegabz'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)

!start timer
CALL system_clock(start,count_rate)
DO n=1,Nt
!fixed point iteration
! store old fluid field terms
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
uold(i,j,k)=u(i,j,k)
uxold(i,j,k)=ux(i,j,k)
uyold(i,j,k)=uy(i,j,k)
uzold(i,j,k)=uz(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
vold(i,j,k)=v(i,j,k)
vxold(i,j,k)=vx(i,j,k)
vyold(i,j,k)=vy(i,j,k)
vzold(i,j,k)=vz(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
wold(i,j,k)=w(i,j,k)
wxold(i,j,k)=wx(i,j,k)
wyold(i,j,k)=wy(i,j,k)
wzold(i,j,k)=wz(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
rhsuhatfix(i,j,k) = (dtInv+(0.5*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*uhat(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
rhsvhatfix(i,j,k) = (dtInv+(0.5*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*vhat(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
rhswhatfix(i,j,k) = (dtInv+(0.5*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*what(i,j,k)
END DO ; END DO ; END DO
! Store old Magnetic field terms
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bxold(i,j,k)=bx(i,j,k)
bxxold(i,j,k)=bxx(i,j,k)
bxyold(i,j,k)=bxy(i,j,k)
bxzold(i,j,k)=bxz(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
byold(i,j,k)=by(i,j,k)
byxold(i,j,k)=byx(i,j,k)
byyold(i,j,k)=byy(i,j,k)
byzold(i,j,k)=byz(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bzold(i,j,k)=bz(i,j,k)
bzxold(i,j,k)=bzx(i,j,k)
bzyold(i,j,k)=bzy(i,j,k)
bzzold(i,j,k)=bzz(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
rhsbxhatfix(i,j,k) = (dtInv+(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*bxhat(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
rhsbyhatfix(i,j,k) = (dtInv+(0.5*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*byhat(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
rhsbzhatfix(i,j,k) = (dtInv+(0.5*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*bzhat(i,j,k)
END DO ; END DO ; END DO

chg=1
DO WHILE (chg .gt. tol)
! Fluid field
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(ux(i,j,k)+uxold(i,j,k))&
+(v(i,j,k)+vold(i,j,k))*(uy(i,j,k)+uyold(i,j,k))&
+(w(i,j,k)+wold(i,j,k))*(uz(i,j,k)+uzold(i,j,k))&
-(bx(i,j,k)+bxold(i,j,k))*(bxx(i,j,k)+bxxold(i,j,k))&
-(by(i,j,k)+byold(i,j,k))*(bxy(i,j,k)+bxyold(i,j,k))&
-(bz(i,j,k)+bzold(i,j,k))*(bxz(i,j,k)+bxzold(i,j,k)))
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_r,nonlinuhat,DECOMP_2D_FFT_FORWARD)
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(vx(i,j,k)+vxold(i,j,k))&
+(v(i,j,k)+vold(i,j,k))*(vy(i,j,k)+vyold(i,j,k))&
+(w(i,j,k)+wold(i,j,k))*(vz(i,j,k)+vzold(i,j,k))&
-(bx(i,j,k)+bxold(i,j,k))*(byx(i,j,k)+byxold(i,j,k))&
-(by(i,j,k)+byold(i,j,k))*(byy(i,j,k)+byyold(i,j,k))&
-(bz(i,j,k)+bzold(i,j,k))*(byz(i,j,k)+byzold(i,j,k)))
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_r,nonlinvhat,DECOMP_2D_FFT_FORWARD)
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(wx(i,j,k)+wxold(i,j,k))&
+(v(i,j,k)+vold(i,j,k))*(wy(i,j,k)+wyold(i,j,k))&
+(w(i,j,k)+wold(i,j,k))*(wz(i,j,k)+wzold(i,j,k))&
-(bx(i,j,k)+bxold(i,j,k))*(bzx(i,j,k)+bzxold(i,j,k))&
-(by(i,j,k)+byold(i,j,k))*(bzy(i,j,k)+bzyold(i,j,k))&
-(bz(i,j,k)+bzold(i,j,k))*(bzz(i,j,k)+bzzold(i,j,k)))
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_r,nonlinwhat,DECOMP_2D_FFT_FORWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
phat(i,j,k)=-1.0d0*( kx(i)*nonlinuhat(i,j,k)&
+ky(j)*nonlinvhat(i,j,k)&
+kz(k)*nonlinwhat(i,j,k))&
/(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)+0.1d0**13)
END DO ; END DO ; END DO

DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
uhat(i,j,k)=(rhsuhatfix(i,j,k)-nonlinuhat(i,j,k)-kx(i)*phat(i,j,k))/&
(dtInv-(0.5d0*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
vhat(i,j,k)=(rhsvhatfix(i,j,k)-nonlinvhat(i,j,k)-ky(j)*phat(i,j,k))/&
(dtInv-(0.5d0*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
what(i,j,k)=(rhswhatfix(i,j,k)-nonlinwhat(i,j,k)-kz(k)*phat(i,j,k))/&
(dtInv-(0.5d0*ReInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes
END DO ; END DO ; END DO

! derivative of u with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=uhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,ux,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=uhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,uy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,uz,DECOMP_2D_FFT_BACKWARD)

! derivative of v with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=vhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,vx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=vhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,vy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=vhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,vz,DECOMP_2D_FFT_BACKWARD)

! derivative of w with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=what(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=what(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=what(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,wz,DECOMP_2D_FFT_BACKWARD)

DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
utemp(i,j,k)=u(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
vtemp(i,j,k)=v(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
wtemp(i,j,k)=w(i,j,k)
END DO ; END DO ; END DO

CALL decomp_2d_fft_3d(uhat,u,DECOMP_2D_FFT_BACKWARD)
CALL decomp_2d_fft_3d(vhat,v,DECOMP_2D_FFT_BACKWARD)
CALL decomp_2d_fft_3d(what,w,DECOMP_2D_FFT_BACKWARD)

DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
u(i,j,k)=u(i,j,k)*scalemodes
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
v(i,j,k)=v(i,j,k)*scalemodes
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
w(i,j,k)=w(i,j,k)*scalemodes
END DO ; END DO ; END DO

! Magnetic field
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(bxx(i,j,k)+bxxold(i,j,k))&
+(v(i,j,k)+vold(i,j,k))*(bxy(i,j,k)+bxyold(i,j,k))&
+(w(i,j,k)+wold(i,j,k))*(bxz(i,j,k)+bxzold(i,j,k))&
-(bx(i,j,k)+bxold(i,j,k))*(ux(i,j,k)+uxold(i,j,k))&
-(by(i,j,k)+byold(i,j,k))*(uy(i,j,k)+uyold(i,j,k))&
-(bz(i,j,k)+bzold(i,j,k))*(uz(i,j,k)+uzold(i,j,k)))
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_r,nonlinuhat,DECOMP_2D_FFT_FORWARD)
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(byx(i,j,k)+byxold(i,j,k))&
+(v(i,j,k)+vold(i,j,k))*(byy(i,j,k)+byyold(i,j,k))&
+(w(i,j,k)+wold(i,j,k))*(byz(i,j,k)+byzold(i,j,k))&
-(bx(i,j,k)+bxold(i,j,k))*(vx(i,j,k)+vxold(i,j,k))&
-(by(i,j,k)+byold(i,j,k))*(vy(i,j,k)+vyold(i,j,k))&
-(bz(i,j,k)+bzold(i,j,k))*(vz(i,j,k)+vzold(i,j,k)))
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_r,nonlinvhat,DECOMP_2D_FFT_FORWARD)
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
temp_r(i,j,k)=0.25d0*((u(i,j,k)+uold(i,j,k))*(bzx(i,j,k)+bzxold(i,j,k))&
+(v(i,j,k)+vold(i,j,k))*(bzy(i,j,k)+bzyold(i,j,k))&
+(w(i,j,k)+wold(i,j,k))*(bzz(i,j,k)+bzzold(i,j,k))&
-(bx(i,j,k)+bxold(i,j,k))*(wx(i,j,k)+wxold(i,j,k))&
-(by(i,j,k)+byold(i,j,k))*(wy(i,j,k)+wyold(i,j,k))&
-(bz(i,j,k)+bzold(i,j,k))*(wz(i,j,k)+wzold(i,j,k)))
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_r,nonlinwhat,DECOMP_2D_FFT_FORWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
phat(i,j,k)=-1.0d0*( kx(i)*nonlinuhat(i,j,k)&
+ky(j)*nonlinvhat(i,j,k)&
+kz(k)*nonlinwhat(i,j,k))&
/(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)+0.1d0**13)
END DO ; END DO ; END DO

DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
bxhat(i,j,k)=(rhsbxhatfix(i,j,k)-nonlinuhat(i,j,k)-kx(i)*phat(i,j,k))/&
(dtInv-(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
byhat(i,j,k)=(rhsbyhatfix(i,j,k)-nonlinvhat(i,j,k)-ky(j)*phat(i,j,k))/&
(dtInv-(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes
END DO ; END DO ; END DO
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
bzhat(i,j,k)=(rhsbzhatfix(i,j,k)-nonlinwhat(i,j,k)-kz(k)*phat(i,j,k))/&
(dtInv-(0.5d0*RemInv)*(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k))) !*scalemodes
END DO ; END DO ; END DO

! derivative of u with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bxhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bxx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bxhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bxy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bxz,DECOMP_2D_FFT_BACKWARD)

! derivative of v with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=byhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,byx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3); DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=byhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,byy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=byhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,byz,DECOMP_2D_FFT_BACKWARD)

! derivative of w with respect to x, y, and z
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bzhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bzx,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bzhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bzy,DECOMP_2D_FFT_BACKWARD)
DO k=decomp%zst(3),decomp%zen(3) ; DO j=decomp%zst(2),decomp%zen(2) ; DO i=decomp%zst(1),decomp%zen(1)
temp_c(i,j,k)=bzhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL decomp_2d_fft_3d(temp_c,bzz,DECOMP_2D_FFT_BACKWARD)

DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bxtemp(i,j,k)=bx(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bytemp(i,j,k)=by(i,j,k)
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bztemp(i,j,k)=bz(i,j,k)
END DO ; END DO ; END DO

CALL decomp_2d_fft_3d(bxhat,bx,DECOMP_2D_FFT_BACKWARD)
CALL decomp_2d_fft_3d(byhat,by,DECOMP_2D_FFT_BACKWARD)
CALL decomp_2d_fft_3d(bzhat,bz,DECOMP_2D_FFT_BACKWARD)

DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bx(i,j,k)=bx(i,j,k)*scalemodes
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
by(i,j,k)=by(i,j,k)*scalemodes
END DO ; END DO ; END DO
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
bz(i,j,k)=bz(i,j,k)*scalemodes
END DO ; END DO ; END DO

mychg(1) =maxval(abs(utemp-u))
mychg(2) =maxval(abs(vtemp-v))
mychg(3) =maxval(abs(wtemp-w))
mychg(4) =maxval(abs(bxtemp-bx))
mychg(5) =maxval(abs(bytemp-by))
mychg(6) =maxval(abs(bztemp-bz))
CALL MPI_ALLREDUCE(mychg,allchg,6,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,ierr)
chg=allchg(1)+allchg(2)+allchg(3)+allchg(4)+allchg(5)+allchg(6)
IF (myid.eq.0) THEN
PRINT *,'chg:',chg
END IF
END DO
time(n+1)=n*dt

IF (myid.eq.0) THEN
PRINT *,'time',n*dt
END IF

!save omegax, omegay, omegaz, omegabx, omegaby, and omegabz
!omegax
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(wy(i,j,k)-vz(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegax'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegay
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(uz(i,j,k)-wx(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegay'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegaz
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(vx(i,j,k)-uy(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegaz'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegabx
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(bzy(i,j,k)-byz(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegabx'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegaby
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(bxz(i,j,k)-bzx(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegaby'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)
!omegabz
DO k=decomp%xst(3),decomp%xen(3); DO j=decomp%xst(2),decomp%xen(2); DO i=decomp%xst(1),decomp%xen(1)
realtemp(i,j,k)=REAL(byx(i,j,k)-bxy(i,j,k),KIND=8)
END DO ; END DO ; END DO
name_config='./data/omegabz'
CALL savedata(Nx,Ny,Nz,n,name_config,realtemp,decomp)

END DO

CALL system_clock(finish,count_rate)

IF (myid.eq.0) then
PRINT *, 'Program took', REAL(finish-start)/REAL(count_rate), 'for main timestepping loop'
END IF

IF (myid.eq.0) THEN
name_config = './data/tdata.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO n=1,1+Nt
WRITE(11,*) time(n)
END DO
CLOSE(11)

name_config = './data/xcoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO i=1,Nx
WRITE(11,*) x(i)
END DO
CLOSE(11)

name_config = './data/ycoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Ny
WRITE(11,*) y(j)
END DO
CLOSE(11)

name_config = './data/zcoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO k=1,Nz
WRITE(11,*) z(k)
END DO
CLOSE(11)
PRINT *,'Saved data'
END IF

! clean up
CALL decomp_2d_fft_finalize
CALL decomp_2d_finalize

DEALLOCATE(x,y,z,time,mychg,allchg,&
u,v,w,ux,uy,uz,vx,vy,vz,&
wx,wy,wz,uold,uxold,&
uyold,uzold,vold,vxold,&
vyold,vzold,wold,wxold,&
wyold,wzold,utemp,vtemp,&
wtemp,temp_r,&
bx,by,bz,bxx,bxy,bxz,&
byx,byy,byz,bzx,bzy,bzz,&
bxold,bxxold,bxyold,bxzold,&
byold,byxold,byyold,byzold,&
bzold,bzxold,bzyold,bzzold,&
bxtemp,bytemp,bztemp,kx,ky,kz,&
uhat,vhat,what,bxhat,&
byhat,bzhat,rhsuhatfix,rhsvhatfix,&
rhswhatfix,rhsbxhatfix,rhsbyhatfix,rhsbzhatfix,&
nonlinuhat,nonlinvhat,nonlinwhat,phat,&
temp_c,realtemp,stat=AllocateStatus)
IF (AllocateStatus .ne. 0) STOP
IF (myid.eq.0) THEN
PRINT *,'Program execution complete'
END IF
CALL MPI_FINALIZE(ierr)

END PROGRAM main


( F)

A subroutine to save real array data for the parallel MPI Fortran program to solve the 3D incompressible magnetohydrodynamics equations Code download
SUBROUTINE savedata(Nx,Ny,Nz,plotnum,name_config,field,decomp)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine saves a three dimensional real array in binary
! format
!
! INPUT
!
! .. Scalars ..
! Nx = number of modes in x - power of 2 for FFT
! Ny = number of modes in y - power of 2 for FFT
! Nz = number of modes in z - power of 2 for FFT
! plotnum = number of plot to be made
! .. Arrays ..
! field = real data to be saved
! name_config = root of filename to save to
!
! .. Output ..
! plotnum = number of plot to be saved
! .. Special Structures ..
! decomp = contains information on domain decomposition
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
! k = loop counter in z direction
! count = counter
! iol = size of file
! .. Arrays ..
! number_file = array to hold the number of the plot
!
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! 2DECOMP&FFT -- Domain decomposition and Fast Fourier Library
! (http://www.2decomp.org/index.html)
! MPI library
USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
IMPLICIT NONE
INCLUDE 'mpif.h'
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny,Nz
INTEGER(KIND=4), INTENT(IN)	:: plotnum
TYPE(DECOMP_INFO), INTENT(IN)	:: decomp
REAL(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
decomp%xst(2):decomp%xen(2),&
decomp%xst(3):decomp%xen(3)), &
INTENT(IN) :: field
CHARACTER*100, INTENT(IN)	:: name_config
INTEGER(kind=4)	:: i,j,k,iol,count,ind
CHARACTER*100	:: number_file

! create character array with full filename
ind = index(name_config,' ') - 1
WRITE(number_file,'(i0)') 10000000+plotnum
number_file = name_config(1:ind)//number_file
ind = index(number_file,' ') - 1
number_file = number_file(1:ind)//'.datbin'
CALL decomp_2d_write_one(1,field,number_file)

END SUBROUTINE savedata


( G)

A makefile to compile the program which finds a numerical solution to the 3D magnetohydrodynamic equations Code download
COMPILER =  mpif90
decompdir= ../2decomp_fft
FLAGS = -O0

DECOMPLIB = -I${decompdir}/include -L${decompdir}/lib -l2decomp_fft
LIBS = #-L${FFTW_LINK} -lfftw3 -lm SOURCES = Magnetohydrodynamics3DfftIMR.f90 savedata.f90 ns3d:$(SOURCES)
${COMPILER} -o mhd3d$(FLAGS) $(SOURCES)$(LIBS) \$(DECOMPLIB)

clean:
rm -f *.o
rm -f *.mod
clobber:
rm -f mhd3d


## 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:

## 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.[9][10], find out what aliasing errors are. Explain why the strategy explained in Johnstone[11] 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[12].
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

$\frac{\partial \mathbf{U}}{\partial t} + \mathbf{U} \cdot \nabla \mathbf{U} = -\nabla p + \Delta \mathbf{U} + \mathbf{J} \times \mathbf{B}$

( 12)

$\frac{\partial \mathbf{H}}{\partial t} - \nabla \times \mathbf{B} = - \mathbf{J}$

( 13)

$\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{H} = 0$

( 14)

$\nabla\cdot\mathbf{B}=0$

( 15)

$\nabla\cdot\mathbf{U}=0$

( 16)

$\mathbf{J}=\sigma\left(\mathbf{H} + \mathbf{U}\times\mathbf{B}\right)$

( 17)

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[13], though Masmoudi[14] shows that some solutions are well behaved in two dimensions. Analysis of the reduced equations is in Duvaut and Lions[15] and can also be found in Gerbeau, Le Bris, and Lelièvre[16]. An example program is in listing H

( H)

A Matlab program which finds a numerical solution to the 3D full magnetohydrodynamic equations Code download
% A program to solve the Full 3D incompressible magnetohydrodynamic equations
% with periodic boundary
% conditions. The program is based on the Orszag-Patterson algorithm as
% documented on pg. 98 of C. Canuto, M.Y. Hussaini, A. Quarteroni and
% T.A. Zhang "Spectral Methods: Evolution to Complex Geometries and
% Applications to Fluid Dynamics" Springer (2007)
%
% The Helmholtz decomposition is used to project the magnetic field onto a
% divergence free subspace. Initial work on this has been done with Damian
% San Roman Alerigi
%
% For plotting, the program uses vol3d, which is available at:
%
% http://www.mathworks.com/matlabcentral/fileexchange/22940-vol3d-v2
%
clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
'defaultaxesfontweight','bold')

% set up grid
tic
Lx = 1;         % period  2*pi*L
Ly = 1;         % period  2*pi*L
Lz = 1;         % period  2*pi*L
Nx = 32;        % number of harmonics
Ny = 32;        % number of harmonics
Nz = 32;        % number of harmonics
Nt = 50;        % number of time slices
dt = 1/Nt;      % time step
t=0;            % initial time
Re = 1.0;       % Reynolds number
Rem = 1.0;      % Magnetic Reynolds number
tol=10^(-10);
% initialise variables
x = (2*pi/Nx)*(-Nx/2:Nx/2 -1)'*Lx;          % x coordinate
kx = 1i*[0:Nx/2-1 0 -Nx/2+1:-1]'/Lx;        % wave vector
y = (2*pi/Ny)*(-Ny/2:Ny/2 -1)'*Ly;          % y coordinate
ky = 1i*[0:Ny/2-1 0 -Ny/2+1:-1]'/Ly;        % wave vector
z = (2*pi/Nz)*(-Nz/2:Nz/2 -1)'*Lz;          % y coordinate
kz = 1i*[0:Nz/2-1 0 -Nz/2+1:-1]'/Lz;        % wave vector
[xx,yy,zz]=meshgrid(x,y,z);
[kxm,kym,kzm]=meshgrid(kx,ky,kz);
[k2xm,k2ym,k2zm]=meshgrid(kx.^2,ky.^2,kz.^2);

% initial conditions for Taylor-Green vortex
% theta=0;
% u=(2/sqrt(3))*sin(theta+2*pi/3)*sin(xx).*cos(yy).*cos(zz);
% v=(2/sqrt(3))*sin(theta-2*pi/3)*cos(xx).*sin(yy).*cos(zz);
% w=(2/sqrt(3))*sin(theta)*cos(xx).*cos(yy).*sin(zz);

% initial condition
sl=1; sk=1; sm=1; lamlkm=sqrt(sl.^2+sk.^2+sm.^2);
u=-0.5*(lamlkm*sl*cos(sk*xx).*sin(sl*yy).*sin(sm.*zz)...
+sm*sk*sin(sk*xx).*cos(sl*yy).*cos(sm.*zz))...
.*exp(-t*(lamlkm^2)/Re);

v=0.5*(lamlkm*sk*sin(sk*xx).*cos(sl*yy).*sin(sm.*zz)...
-sm*sl*cos(sk*xx).*sin(sl*yy).*cos(sm.*zz))...
*exp(-t*(lamlkm^2)/Re);

w=cos(sk*xx).*cos(sl*yy).*sin(sm*zz)*exp(-t*(lamlkm^2)/Re);

uhat=fftn(u);
vhat=fftn(v);
what=fftn(w);

ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm);
vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm);
wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm);

% Magnetic field
bx=sin(sl*yy).*sin(sm.*zz);
by=sin(sm.*zz);
bz=cos(sk*xx).*cos(sl*yy);

bxhat=fftn(bx);
byhat=fftn(by);
bzhat=fftn(bz);

bxx=ifftn(bxhat.*kxm);bxy=ifftn(bxhat.*kym);bxz=ifftn(bxhat.*kzm);
byx=ifftn(byhat.*kxm);byy=ifftn(byhat.*kym);byz=ifftn(byhat.*kzm);
bzx=ifftn(bzhat.*kxm);bzy=ifftn(bzhat.*kym);bzz=ifftn(bzhat.*kzm);

% Electric field
ex=sin(sl*yy).*sin(sm.*zz);
ey=sin(sm.*zz);
ez=cos(sk*xx).*cos(sl*yy);

exhat=fftn(ex);
eyhat=fftn(ey);
ezhat=fftn(ez);

exx=ifftn(exhat.*kxm);exy=ifftn(exhat.*kym);exz=ifftn(exhat.*kzm);
eyx=ifftn(eyhat.*kxm);eyy=ifftn(eyhat.*kym);eyz=ifftn(eyhat.*kzm);
ezx=ifftn(ezhat.*kxm);ezy=ifftn(ezhat.*kym);ezz=ifftn(ezhat.*kzm);

% Current density
jx=ex+uy.*bz-uz.*by;
jy=ey+uz.*bx-ux.*bz;
jz=ez+ux.*by-uy.*bx;

jxhat=fftn(jx);
jyhat=fftn(jy);
jzhat=fftn(jz);

% calculate fluid vorticity for plotting
omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy;
omegatot=omegax.^2+omegay.^2+omegaz.^2;
% calculate magnetic vorticity for plotting
omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy;
omegabtot=omegabx.^2+omegaby.^2+omegabz.^2;
figure(1); clf; n=0;
subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]);
H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]);
H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

for n=1:Nt
uold=u; uxold=ux; uyold=uy; uzold=uz;
vold=v; vxold=vx; vyold=vy; vzold=vz;
wold=w; wxold=wx; wyold=wy; wzold=wz;

bxold=bx; byold=by; bzold=bz;

exold=ex; exxold=exx; exyold=exy; exzold=exz;
eyold=ey; eyxold=eyx; eyyold=eyy; eyzold=eyz;
ezold=ez; ezxold=ezx; ezyold=ezy; ezzold=ezz;

jxold=jx; jyold=jy; jzold=jz;

rhsuhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*uhat;
rhsvhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*vhat;
rhswhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*what;

rhsbxhatfix=bxhat-0.5*dt*(kym.*ezhat-kzm.*eyhat);
rhsbyhatfix=byhat-0.5*dt*(kzm.*exhat-kxm.*ezhat);
rhsbzhatfix=bzhat-0.5*dt*(kxm.*eyhat-kym.*exhat);

rhsexhatfix=exhat+0.5*dt*(kym.*bzhat-kzm.*byhat)-0.5*dt*jxhat;
rhseyhatfix=eyhat+0.5*dt*(kzm.*bxhat-kxm.*bzhat)-0.5*dt*jyhat;
rhsezhatfix=ezhat+0.5*dt*(kxm.*byhat-kym.*bxhat)-0.5*dt*jzhat;

chg=1; t=t+dt;
while (chg>tol)
nonlinu=0.25*( (u+uold).*(ux+uxold)...
+(v+vold).*(uy+uyold)...
+(w+wold).*(uz+uzold)...
-(jy+jyold).*(bz+bzold)...
+(jz+jzold).*(by+byold));
nonlinv=0.25*( (u+uold).*(vx+vxold)...
+(v+vold).*(vy+vyold)...
+(w+wold).*(vz+vzold)...
+(jx+jxold).*(bz+bzold)...
-(jz+jzold).*(bx+bxold));
nonlinw=0.25*( (u+uold).*(wx+wxold)...
+(v+vold).*(wy+wyold)...
+(w+wold).*(wz+wzold)...
-(jx+jxold).*(by+byold)...
+(jy+jyold).*(bx+bxold));
nonlinuhat=fftn(nonlinu);
nonlinvhat=fftn(nonlinv);
nonlinwhat=fftn(nonlinw);
phat=-1.0*(kxm.*nonlinuhat+kym.*nonlinvhat+kzm.*nonlinwhat)...
./(k2xm+k2ym+k2zm+0.1^13);
uhat=(rhsuhatfix-nonlinuhat-kxm.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
vhat=(rhsvhatfix-nonlinvhat-kym.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
what=(rhswhatfix-nonlinwhat-kzm.*phat)...
./(1/dt - (0.5/Re)*(k2xm+k2ym+k2zm));
ux=ifftn(uhat.*kxm);uy=ifftn(uhat.*kym);uz=ifftn(uhat.*kzm);
vx=ifftn(vhat.*kxm);vy=ifftn(vhat.*kym);vz=ifftn(vhat.*kzm);
wx=ifftn(what.*kxm);wy=ifftn(what.*kym);wz=ifftn(what.*kzm);
utemp=u; vtemp=v; wtemp=w;
u=ifftn(uhat); v=ifftn(vhat); w=ifftn(what);

bxhat=rhsbxhatfix-0.5*dt*(kym.*ezhat-kzm.*eyhat);
byhat=rhsbyhatfix-0.5*dt*(kzm.*exhat-kxm.*ezhat);
bzhat=rhsbzhatfix-0.5*dt*(kxm.*eyhat-kym.*exhat);

bxtemp=bx; bytemp=by; bztemp=bz;
bx=ifftn(bxhat); by=ifftn(byhat); bz=ifftn(bzhat);

exhat=rhsexhatfix+0.5*dt*(kym.*bzhat-kzm.*byhat-jxhat);
eyhat=rhseyhatfix+0.5*dt*(kzm.*bxhat-kxm.*bzhat-jyhat);
ezhat=rhsezhatfix+0.5*dt*(kxm.*byhat-kym.*bxhat-jzhat);

extemp=ex; eytemp=ey; eztemp=ez;
ex=ifftn(exhat); ey=ifftn(eyhat); ez=ifftn(ezhat);

jx=ex+uy.*bz-uz.*by;
jy=ey+uz.*bx-ux.*bz;
jz=ez+ux.*by-uy.*bx;

jxhat=fftn(jx);
jyhat=fftn(jy);
jzhat=fftn(jz);

chg=max(abs(utemp-u))   +max(abs(vtemp-v))  +max(abs(wtemp-w))+...
max(abs(bxtemp-bx)) +max(abs(bytemp-by))+max(abs(bztemp-bz))+...
max(abs(extemp-ex)) +max(abs(eytemp-ey))+max(abs(eztemp-ez));
end
% calculate vorticity for plotting
omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy;
omegatot=omegax.^2+omegay.^2+omegaz.^2;

% calculate magnetic vorticity for plotting
omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy;
omegabtot=omegabx.^2+omegaby.^2+omegabz.^2;
figure(1); clf;
subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]);
H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]);
H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z);
xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
axis equal; axis square; view(3);

end

1. A compressible magnetohydrodynamic model is

$\frac{\partial\rho}{\partial t} +\nabla\left(\rho\mathbf{u}\right)=0$

( 18)

$\frac{\partial\left(\rho\mathbf{u}\right)}{\partial t} +\nabla\left(\rho\mathbf{u}\mathbf{u}^T\right)+\nabla p= \nabla\cdot\mathbf{\tau} +\mathbf{J}\times\mathbf{B}$

( 19)

$\frac{\partial\left(\rho E\right)}{\partial t} +\nabla\left(\rho E \mathbf{u}\right)+\nabla\cdot(p \mathbf{u})= -\nabla\cdot\mathbf{q} + \nabla\cdot(\mathbf{\tau}\mathbf{u})$

( 20)

$\frac{\partial \mathbf{H}}{\partial t} - \nabla \times \mathbf{B} = - \mathbf{J}$

( 21)

$\frac{\partial \mathbf{B}}{\partial t} + \nabla \times \mathbf{H} = 0$

( 22)

$\nabla\cdot\mathbf{B}=0$

( 23)

$\mathbf{J}=\sigma\left(\mathbf{H} + \mathbf{U}\times\mathbf{B}\right)$

( 24)

$p=\rho RT$

( 25)

$\mathbf{q}=-K\nabla T$

( 26)

$E=C_vT +\frac{\mathbf{u}\cdot\mathbf{u}}{2}+ \frac{\mathbf{H}\cdot\mathbf{H}}{2}+ \frac{\mathbf{B}\cdot\mathbf{B}}{2}$

( 27)

$\mathbf{\tau}=\mu\left[ \nabla\mathbf{u} +(\nabla\mathbf{u})^T \right] +\lambda\left(\nabla\cdot\mathbf{u}\right)\mathbf{I}$

( 28)

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.

## References

Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. (2006). Spectral Methods: Fundamentals in Single Domains. Springer.

Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. (2007). Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics. Springer.

Carbone, V.; Pouquet, A. (2009). "An Introduction to Fluid and MHD Turbulence for Astrophysical Flows: Theory, Observational and Numerical Data, and Modeling". Lecture Notes in Physics (Springer-Verlag) 778: 71-128. doi:10.1007/978-3-642-00210-6_3.

Chandrasekhar, Subrahmanyan (1961). Hydrodynamic and Hyrdomagnetic Stability. Dover. ISBN 0-486-64071-X.

Davidson, P.A. (2001). An Introduction to Magnetohydrodynamics. Cambridge University Press. doi:10.1017/CBO9780511626333.

Dorch, Soren (2007). "Magnetohydrodynamic". Scholarpedia 4: 2295. doi:10.4249/scholarpedia.2295.

Duvaut, G.; Lions, J.L. (1972). "Inéquations en thermoélasticité et magnétohydrodynamique". Archive for Rational Mechanics and Analysis 46: 241-279. doi:10.1007/BF00250512.

Gerbeau, J.-F.; Le Bris, C.; Lelièvre, T. (2006). Mathematical Methods for the Magnetohydrodynamics of Liquid Metals. Oxford University Press. doi:10.1093/acprof:oso/9780198566656.001.0001.

Johnstone, R. (2012). "Improved Scaling for Direct Numerical Simulations of Turbulence". HECTOR distributed Computational Science and Engineering Report.

Masmoudi, N. (2010). "Global well posedness for the Maxwell-Navier-Stokes system in 2D". J. Math. Pures. Appl. 93: 559-571. doi:10.1017/S002211207900210X.

Orszag, S.; Tang, C.M. (1979). "Small-scale structure of two-dimensional magnetohydrodynamic turbulence". Journal of Fluid Mechanics 90 (1): 129-143. doi:10.1017/S002211207900210X.

Schnack, D.D. (2009). Lectures in Magnetohydrodynamics. Springer. doi:10.1007/978-3-642-00688-3.

Sermange, M.; Temam, R. (1983). "Some mathematical questions related to the MHD equations". Communications on Pure and Applied Mathematics 36 (5): 635-664. doi:10.1002/cpa.3160360506.

Verma, Mahendra (2004). "Statistical theory of magnetohydrodynamic turbulence: recent results". Physics Reports 401: 220-380. doi:10.1016/j.physrep.2004.07.007.

Zheng, Linjin, ed (2009). Topics in Magnetohydrodynamics. InTech. doi:10.5772/2080.