并行谱数值方法/不可压缩磁流体力学

不可压缩磁流体力学

背景

不可压缩磁流体力学方程类似于纳维-斯托克斯方程，但可以显示各种其他现象。本节假设读者熟悉纳维-斯托克斯情况的材料 此处。 不可压缩磁流体力学方程的数学研究也像纳维-斯托克斯方程一样糟糕。 可以在 Carbone 和 Pouquet 、Chandrasekhar 、Davidson 、Gerbeau、Le Bris 和 Lelièvre 、Schnak 和 Verma 中找到对这些方程的物理介绍。 可以在此处找到更一般的概述 。 在流体速度明显小于光速以及流体几乎不可压缩的假设下，无量纲方程表示为

${\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)

在这些方程式中，$\mathbf {u} (x,y,z)=(u,v,w)$ 是流体速度，其分量在 $x$、$y$ 和 $z$ 方向，$\mathbf {b} (x,y,z)=(bx,by,bz)$ 是磁场，其分量在 $x$、$y$ 和 $z$ 方向，$p$ 是压强场。方程式 1 代表动量守恒，方程式 2 对应于麦克斯韦方程式，用于非相对论性流动，其中电场的时变被忽略（参见 Verma），方程式 3 是连续性方程式，代表不可压缩流体的质量守恒，方程式 4 代表没有磁单极子源。

二维情况

我们将首先考虑二维情况。模拟不可压缩磁流体动力学方程的一个困难是数值上满足方程式 3 和 4 中的约束，这些约束有时被称为无散度条件或螺线约束。为了自动满足这些约束，在二维情况下，其中

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

和

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

可以使用势能公式重新编写这些方程。在这种情况下，我们令

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

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

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

和

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

其中 $\psi (x,y)$ 是流函数，而 $\phi (x,y)$ 是磁势。注意，

$\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$,

和

$\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$,

因此，公式 3 和 4 会自动满足。通过这种变量变换，我们可以对动量方程 1 和磁输运方程 2 进行旋度运算，从而得到一个简化的包含两个耦合偏微分方程的系统。我们定义流体涡度 $\omega$，使得

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

以及磁涡度，使得

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

公式 1 和 2 变为

${\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)

和

$\Delta \phi =-\alpha$.

( 8)

其中 $fx$ 和 $fy$ 代表力的 $x$ 和 $y$ 分量 $\mathbf {f}$.

三维情况

在三维情况下，方程的势函数表达并不方便。因此，引入磁压的概念，以便为执行无散度磁场约束提供一个简单的算法。那么方程变为：

${\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)

和

${\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)

这两个方程都可以用隐式中点规则离散化。这里 $Re$ 代表流体的雷诺数，$Rem$ 代表磁雷诺数。流体压力由以下公式给出

$\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}$

磁“压力”由以下公式给出

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

串行程序

我们首先编写 Matlab 程序来演示如何在单处理器上求解这些方程。第一个程序使用隐式中点规则时间步进求解二维磁流体动力学方程，并在列表 A 中。在这个程序中，我们还求解了一个具有扩散常数 $D_{\theta }$ 的被动标量的方程

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

( 11)

( A)

一个 Matlab 程序，它找到二维磁流体动力学方程的数值解 代码下载

% 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)

一个 Python 程序，它找到二维磁流体动力学方程的数值解 代码下载

#!/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

第二个程序使用隐式中点规则对三维不可压缩磁流体动力学方程进行时间步进，它位于列表 B 中。

( B)

一个 Matlab 程序，它找到三维磁流体动力学方程的数值解 代码下载

% 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)

一个 Python 程序，它找到三维磁流体动力学方程的数值解 代码下载

#!/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 程序

本节中的程序基于 https://wikibooks.cn/wiki/Parallel_Spectral_Numerical_Methods/The_Two-_and_Three-Dimensional_Navier-Stokes_Equations 中的 Navier-Stokes 方程。

( C)

一个串行 Fortran 程序，它找到二维磁流体动力学方程的数值解 代码下载

	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)

一个用于编译求解二维磁流体动力学方程的数值解程序的 makefile 代码下载

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

并行 Fortran 程序

为了编写并行程序，我们将使用 2DECOMP&FFT 库进行域分解和傅里叶变换。2DECOMP&FFT 可从 http://www.2decomp.org/index.html 获取。该网站包含一些示例，说明了如何使用该库，特别是 http://www.2decomp.org/case_study1.html 中的示例代码，它非常清楚地说明了如何将使用 FFTW 的代码转换为使用 MPI 和上述库的代码。为了使程序正常运行，首先编译 2DECOMP&FFT，更新 makefile 中的链接，然后编译并行 Fortran 程序。

( E)

一个并行 Fortran 程序，它找到三维磁流体动力学方程的数值解 代码下载

	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)

一个用于保存并行 MPI Fortran 程序中三维不可压缩磁流体动力学方程解的实数数组数据的子程序 代码下载

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)

一个用于编译程序的 Makefile，该程序可以找到三维磁流体动力学方程的数值解 代码下载

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

其他公开可用的用于求解磁流体动力学方程的傅里叶谱方法程序

另一个公开可用的基于快速傅里叶变换的用于求解不可压缩磁流体动力学方程的程序是

Tarang，可从 http://turbulencehub.org/index.php/codes/tarang/ 获取。

练习

1. 我们编写的三维程序并非最高效的，因为可以使用实数到复数变换将工作量减半。在其中一个 Navier-Stokes 程序中实现一个实数到复数变换。
2. 我们编写的程序也可能引入一些混叠误差。阅读有关谱方法的书籍，例如 Canuto 等人的书，找出什么是混叠误差。解释 Johnstone 中解释的策略如何减少混叠误差。
3. 看看你是否可以找到其他用于求解不可压缩磁流体动力学方程的开源快速傅里叶变换求解器，并将它们添加到上面的列表中。
4. 有一些建议表明，不可压缩磁流体动力学方程可以通过磁场或速度场导数的其中一个范数变为无穷大而表现出奇异行为。进行模拟实验，看看你是否能找到这种现象的证据。在这一领域报告计算实验的第一批论文之一是 Orszag 和 Tang。
5. 将上面的程序之一翻译成其他语言。以下网站可能有所帮助：http://rosettacode.org/wiki/Rosetta_Code
6. 完整的磁流体动力学方程为

${\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)

其中 $\mathbf {U}$ 是速度，$\mathbf {p}$ 是压力，$\mathbf {J}$ 是电流密度，$\mathbf {B}$ 是磁场，$\mathbf {H}$ 是电场。模拟这些完整方程在二维和三维中的情况 - 目前尚不清楚这些方程的解在三维中是否始终存在，参见 Germain、Ibrahim 和 Masmoudi，尽管 Masmoudi 表明某些解在二维中表现良好。简化方程的分析见 Duvaut 和 Lions，也可以在 Gerbeau、Le Bris 和 Lelièvre 中找到。示例程序在列表 H 中。

( H)

一个用于求解三维完整磁流体动力学方程的数值解的 Matlab 程序 代码下载

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

    jxold=jx; jyold=jy; jzold=jz;
    
    rhsuhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*uhat;
    rhsvhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*vhat;
    rhswhatfix=(1/dt+(0.5/Re)*(k2xm+k2ym+k2zm)).*what;
    
    rhsbxhatfix=bxhat-0.5*dt*(kym.*ezhat-kzm.*eyhat);
    rhsbyhatfix=byhat-0.5*dt*(kzm.*exhat-kxm.*ezhat);
    rhsbzhatfix=bzhat-0.5*dt*(kxm.*eyhat-kym.*exhat);

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

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

        exhat=rhsexhatfix+0.5*dt*(kym.*bzhat-kzm.*byhat-jxhat);
        eyhat=rhseyhatfix+0.5*dt*(kzm.*bxhat-kxm.*bzhat-jyhat);
        ezhat=rhsezhatfix+0.5*dt*(kxm.*byhat-kym.*bxhat-jzhat);
        
        extemp=ex; eytemp=ey; eztemp=ez;
        ex=ifftn(exhat); ey=ifftn(eyhat); ez=ifftn(ezhat);
 
        jx=ex+uy.*bz-uz.*by;
        jy=ey+uz.*bx-ux.*bz;
        jz=ez+ux.*by-uy.*bx;

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

        chg=max(abs(utemp-u))   +max(abs(vtemp-v))  +max(abs(wtemp-w))+...
            max(abs(bxtemp-bx)) +max(abs(bytemp-by))+max(abs(bztemp-bz))+...
            max(abs(extemp-ex)) +max(abs(eytemp-ey))+max(abs(eztemp-ez));
    end
    % calculate vorticity for plotting
    omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy; 
    omegatot=omegax.^2+omegay.^2+omegaz.^2;
    
    % calculate magnetic vorticity for plotting
    omegabx=bxy-byz; omegaby=bxz-bzx; omegabz=byx-bxy; 
    omegabtot=omegabx.^2+omegaby.^2+omegabz.^2;
    figure(1); clf; 
    subplot(2,1,1); title(['|Fluid omega|^2 ',num2str(t)]);
    H = vol3d('CData',omegatot,'texture','3D','XData',x,'YData',y,'ZData',z);
    xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
    axis equal; axis square; view(3); 
      
    subplot(2,1,2); title(['|Magnetic omega|^2 ',num2str(t)]);
    H = vol3d('CData',omegabtot,'texture','3D','XData',x,'YData',y,'ZData',z);
    xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
    axis equal; axis square; view(3); 
       
end

1. 一个可压缩磁流体动力学模型是

${\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)

其中 $\mathbf {U}$ 是速度，$\mathbf {p}$ 是压力，$\mathbf {J}$ 是电流密度，$\mathbf {B}$ 是磁场，$\mathbf {H}$ 是电场，$\mathbf {E}$ 是能量密度，$\mathbf {\tau }$ 是粘性流体应力，$T$ 是温度场，$C_{v}$ 是定容比热容，$R$ 是理想气体常数，$\mathbf {q}$ 是热通量。在二维和三维空间中模拟这些完整方程。

参考文献

Canuto, C.；Hussaini, M.Y.；Quarteroni, A.；Zang, T.A. (2006)。谱方法：单域基础。施普林格。

Canuto, C.；Hussaini, M.Y.；Quarteroni, A.；Zang, T.A. (2007)。谱方法：向复杂几何和流体力学应用的演变。施普林格。

Carbone, V.；Pouquet, A. (2009)。“天体物理流体和 MHD 湍流简介：理论、观测和数值数据以及建模”。空间等离子体湍流。物理讲义。第 778 卷。施普林格出版社。第 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)

钱德拉塞卡，苏布拉马尼安 (1961)。流体动力学和磁流体力学稳定性。多佛。 ISBN 0-486-64071-X. {{cite book}}: Check |authorlink= value (help); External link in |authorlink= (help)

Davidson, P.A. (2001). 磁流体力学导论. 劍橋大學出版社. doi:10.1017/CBO9780511626333. ISBN 9780521794879. {{cite book}}: Check |authorlink= value (help); External link in |authorlink= (help)

Dorch, Soren (2007). "磁流体力学". 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). "热弹性与磁流体力学中的不等式". 理性力学与分析档案. 46 (4): 241–279. doi:10.1007/BF00250512. S2CID 120571424.

Gerbeau, J.-F.; Le Bris, C.; Lelièvre, T. (2006). 液态金属磁流体力学的数学方法. 牛津大学出版社. 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). "Navier-Stokes-Maxwell 方程组的适定性". arΧiv:1207.6187 [math.AP]. 

Johnstone, R. (2012). "湍流直接数值模拟的改进缩放". HECTOR 分布式计算科学与工程报告.

Masmoudi, N. (2010). "二维 Maxwell-Navier-Stokes 系统的全局适定性" (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). "二维磁流体力学湍流的小尺度结构". 流体力学杂志. 90 (1): 129–143. doi:10.1017/S002211207900210X. S2CID 123354420.

Schnack, D.D. (2009). 磁流体力学讲义. 物理学讲义. 第 780 卷. 施普林格. doi:10.1007/978-3-642-00688-3. ISBN 978-3-642-00687-6.

Sermange, M.; Temam, R. (1983). "与 MHD 方程相关的某些数学问题" (PDF). 纯数学与应用数学通讯. 36 (5): 635–664. doi:10.1002/cpa.3160360506. {{cite journal}}: Check |authorlink2= value (help); External link in |authorlink2= (help)

Verma, Mahendra (2004). "磁流体动力学湍流的统计理论：最新成果". 物理报告. 401 (5–6): 220–380. arXiv:nlin/0404043. doi:10.1016/j.physrep.2004.07.007. S2CID 119352240. {{cite journal}}: 检查 |authorlink= 值 (帮助); 外部链接在 |authorlink= (帮助)

"磁流体动力学".

郑林进，编 (2009). 磁流体动力学主题. InTech. doi:10.5772/2080. ISBN 978-953-51-0211-3.