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

More information on visualization can be found on the Mayavi
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
	!
	! FURTHER COMMENTS
	! 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
	!
	! FURTHER COMMENTS
	!
	! 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
! see http://www.2decomp.org/ for more information
! 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
!
! FURTHER COMMENTS
!--------------------------------------------------------------------
! 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:

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

Exercises

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.
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.
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.
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].
Translate one of the programs above to another language. The following site may help, http://rosettacode.org/wiki/Rosetta_Code
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

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_{v}T+{\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.

Notes

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". Turbulence in Space Plasmas. Lecture Notes in Physics. Vol. 778. Springer-Verlag. pp. 71–128. doi:10.1007/978-3-642-00210-6_3. ISBN 978-3-642-00209-0. {{cite book}}: Check |authorlink2= value (help); External link in |authorlink2= (help)

Chandrasekhar, Subrahmanyan (1961). Hydrodynamic and Hyrdomagnetic Stability. Dover. ISBN 0-486-64071-X. {{cite book}}: Check |authorlink= value (help); External link in |authorlink= (help)

Davidson, P.A. (2001). An Introduction to Magnetohydrodynamics. Cambridge University Press. doi:10.1017/CBO9780511626333. ISBN 9780521794879. {{cite book}}: Check |authorlink= value (help); External link in |authorlink= (help)

Dorch, Soren (2007). "Magnetohydrodynamic". Scholarpedia. 2 (4): 2295. doi:10.4249/scholarpedia.2295. {{cite journal}}: Check |authorlink= value (help); External link in |authorlink= (help)

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

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. ISBN 978-0-19-856665-6. {{cite book}}: Check |authorlink2= value (help); Check |authorlink3= value (help); Check |authorlink= value (help); External link in |authorlink2=, |authorlink3=, and |authorlink= (help)

Germain, P.; Ibrahim, S.; Masmoudi, N. (2012). "Well Posedness of the Navier-Stokes-Maxwell Equations". arΧiv:1207.6187 [math.AP].

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" (PDF). J. Math. Pures. Appl. 93: 559–571. doi:10.1017/S002211207900210X. S2CID 123354420. {{cite journal}}: Check |authorlink= value (help); External link in |authorlink= (help)

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. S2CID 123354420.

Schnack, D.D. (2009). Lectures in Magnetohydrodynamics. Lecture Notes in Physics. Vol. 780. Springer. doi:10.1007/978-3-642-00688-3. ISBN 978-3-642-00687-6.

Sermange, M.; Temam, R. (1983). "Some mathematical questions related to the MHD equations" (PDF). Communications on Pure and Applied Mathematics. 36 (5): 635–664. doi:10.1002/cpa.3160360506. {{cite journal}}: Check |authorlink2= value (help); External link in |authorlink2= (help)

Verma, Mahendra (2004). "Statistical theory of magnetohydrodynamic turbulence: recent results". Physics Reports. 401 (5–6): 220–380. arXiv:nlin/0404043. doi:10.1016/j.physrep.2004.07.007. S2CID 119352240. {{cite journal}}: Check |authorlink= value (help); External link in |authorlink= (help)

"Magnetohydrodynamics".

Zheng, Linjin, ed. (2009). Topics in Magnetohydrodynamics. InTech. doi:10.5772/2080. ISBN 978-953-51-0211-3.

[1] Carbone and Pouquet (2009)

[2] Chandrasekhar (1961)

[3] Davidson (2001)

[4] Gerbeau, Le Bris and Lelièvre (2006)

[5] Schnak (2009)

[6] Verma (2004)

[7] Magnetohydrodynamics

[8] Verma (2004)

[9] Canuto et al (2006)

[10] Canuto et al (2007)

[11] Johnstone (2012)

[12] Orszag and Tang (1979)

[13] Germain, Ibrahim and Masmoudi (2012)

[14] Masmoudi (2010)

[15] Duvaut and Lions (1972)

[16] Gerbeau, Le Bris, and Lelièvre (2006)

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]