二維和三維 Navier-Stokes 方程

背景

Navier-Stokes 方程描述了流體的運動。為了推匯出 Navier-Stokes 方程，我們假設流體是連續體（不是由單個粒子組成，而是連續的物質），並且質量和動量守恆。在做出一些假設並在不可壓縮流體粒子上應用牛頓第二定律之後，可以完整地推匯出 Navier-Stokes 方程。所有細節都省略了，因為有很多關於此資訊的來源，兩個特別清晰的來源是 Tritton^[1] 和 Doering 和 Gibbon^[2]; Gallavotti^[3] 還應注意引入這些方程的數學和物理方面，以及 Uecker^[4] 包括快速推導和一些示例傅立葉譜 Matlab 程式碼。有關 Navier-Stokes 方程譜方法的更詳細介紹，請參閱 Canuto 等人^[5]^[6]。不可壓縮 Navier-Stokes 方程為

\rho \left({\frac {\partial \mathbf {u} }{\partial t}}+\mathbf {u} \cdot \nabla \mathbf {u} \right)=-\nabla p+\mu \Delta \mathbf {u} +\mathbf {f}

( 1)

\nabla \cdot \mathbf {u} =0

.

( 2)

在這些方程中， $\rho$ 是密度， $\mathbf {u} (x,y,z)=(u,v,w)$ 是速度，其分量分別在 $x$ 、 $y$ 和 $z$ 方向， $p$ 是壓強場， $\mu$ 是動態粘度（在不可壓縮情況下為常數），而 $\mathbf {f}$ 是體力（作用於整個體積的力）。方程 1 表示動量守恆，而方程 2 是連續性方程，表示不可壓縮流體的質量守恆。

二維情況

我們將首先考慮二維情況。模擬不可壓縮 Navier-Stokes 方程的一個困難是方程 2 中不可壓縮約束的數值滿足，這有時被稱為無散度條件或螺線管約束。為了在二維情況下自動滿足此不可壓縮約束，其中

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

可以使用不同的公式重新編寫這些方程，即流函式渦度公式。在這種情況下，令

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

以及

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

其中 $\psi (x,y)$ 是流函式。流函式的等值線代表流線（流線是一條連續曲線，其上瞬時速度是切線，關於流線的更多資訊，請參閱 Tritton^[7]）。需要注意的是

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

因此，公式 2 自動滿足。進行這種變數替換後，透過對動量方程（公式 1 ）進行旋度運算，可以得到一個標量偏微分方程。定義渦度 $\omega$ ，使得

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

而公式 1 變成

${\frac {\partial }{\partial x}}\left[\rho \left({\frac {\partial v}{\partial t}}+u{\frac {\partial v}{\partial x}}+v{\frac {\partial v}{\partial y}}\right)\right]-{\frac {\partial }{\partial y}}\left[\rho \left({\frac {\partial u}{\partial t}}+u{\frac {\partial u}{\partial x}}+v{\frac {\partial u}{\partial y}}\right)\right]$ $={\frac {\partial }{\partial x}}\left[\mu \left({\frac {\partial ^{2}v}{\partial x^{2}}}+{\frac {\partial ^{2}v}{\partial y^{2}}}\right)+fy\right]-{\frac {\partial }{\partial y}}\left[\mu \left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)+fx\right]$

其中 $fx$ 和 $fy$ 代表力 $\mathbf {f}$ 的 $x$ 和 $y$ 分量。由於流動是無散度的，

${\frac {\partial u}{\partial x}}=-{\frac {\partial v}{\partial y}}$

因此可以簡化非線性項得到

${\frac {\partial }{\partial x}}\left[\left(u{\frac {\partial v}{\partial x}}+v{\frac {\partial v}{\partial y}}\right)\right]-{\frac {\partial }{\partial y}}\left[\left(u{\frac {\partial u}{\partial x}}+v{\frac {\partial u}{\partial y}}\right)\right]$

$={\frac {\partial u}{\partial x}}{\frac {\partial v}{\partial x}}+u{\frac {\partial ^{2}v}{\partial x^{2}}}+{\frac {\partial v}{\partial x}}{\frac {\partial v}{\partial y}}+v{\frac {\partial ^{2}v}{\partial x\partial y}}-{\frac {\partial u}{\partial y}}{\frac {\partial u}{\partial x}}-u{\frac {\partial ^{2}u}{\partial x\partial y}}-{\frac {\partial v}{\partial y}}{\frac {\partial u}{\partial y}}-v{\frac {\partial ^{2}u}{\partial y^{2}}}$

$=u\left({\frac {\partial ^{2}v}{\partial x^{2}}}-{\frac {\partial ^{2}u}{\partial x\partial y}}\right)+v\left({\frac {\partial ^{2}v}{\partial x\partial y}}-{\frac {\partial ^{2}u}{\partial ^{2}y}}\right)$ .

最後我們得到

\rho \left({\frac {\partial \omega }{\partial t}}+u{\frac {\partial \omega }{\partial x}}+v{\frac {\partial \omega }{\partial y}}\right)=\mu \Delta \omega +{\frac {\partial fy}{\partial x}}-{\frac {\partial fx}{\partial y}}

( 3)

以及

\Delta \psi =-\omega

.

( 4)

需要注意的是，在這個公式中，Navier-Stokes 方程類似於渦度的強制熱方程，具有非區域性和非線性項。我們可以利用這種結構，透過修改熱方程的近似解的數值程式來找到數值解。

該方程的一個簡單時間離散化方法是 Crank-Nicolson 方法，其中非線性項使用不動點迭代求解。關於 Navier-Stokes 方程時間離散化方案收斂性的教程可以在 Temam^[8] 中找到。時間離散方程變為

\rho \left[{\frac {\omega ^{n+1,k+1}-\omega ^{n}}{\delta t}}\right.

( 5)

$\left.+{\frac {1}{2}}\left(u^{n+1,k}{\frac {\partial \omega ^{n+1,k}}{\partial x}}+v^{n+1,k}{\frac {\partial \omega ^{n+1,k}}{\partial y}}+u^{n}{\frac {\partial \omega ^{n}}{\partial x}}+v^{n}{\frac {\partial \omega ^{n}}{\partial y}}\right)\right]$

$={\frac {\mu }{2}}\Delta \left(\omega ^{n+1,k+1}+\omega ^{n}\right)+\left.\left({\frac {\partial fx}{\partial y}}-{\frac {\partial fy}{\partial x}}\right)\right|_{t=(n+0.5)\delta t}$ ,

以及

\Delta \psi ^{n+1,k+1}=-\omega ^{n+1,k+1}

,

( 6)

$u^{n+1,k+1}={\frac {\partial \psi ^{n+1,k+1}}{\partial y}}$ ,

$v^{n+1,k+1}=-{\frac {\partial \psi ^{n+1,k+1}}{\partial x}}$ .

在這些方程中，上標 $n$ 表示時間步長，上標 $k$ 表示迭代次數。另一種時間離散化的選擇是隱式中點法則，它給出，

\rho \left[{\frac {\omega ^{n+1,k+1}-\omega ^{n}}{\delta t}}\right.

( 7)

$\left.+\left({\frac {u^{n+1,k}+u^{n}}{2}}\right){\frac {\partial }{\partial x}}\left({\frac {\omega ^{n+1,k}+\omega ^{n}}{2}}\right)+\left({\frac {v^{n+1,k}+v^{n}}{2}}\right){\frac {\partial }{\partial y}}\left({\frac {\omega ^{n+1,k}+\omega ^{n}}{2}}\right)\right]$

$={\frac {\mu }{2}}\Delta \left(\omega ^{n+1,k+1}+\omega ^{n}\right)+\left.\left({\frac {\partial fx}{\partial y}}-{\frac {\partial fy}{\partial x}}\right)\right|_{t=(n+0.5)\delta t}$ ,

以及

\Delta \psi ^{n+1,k+1}=-\omega ^{n+1,k+1}

,

( 8)

$u^{n+1,k+1}={\frac {\partial \psi ^{n+1,k+1}}{\partial y}}$ ,

$v^{n+1,k+1}=-{\frac {\partial \psi ^{n+1,k+1}}{\partial x}}$ .

三維情況

這裡 $\mathbf {u} =(u(x,y,z,t),v(x,y,z,t),w(x,y,z,t))$ - 不幸的是，目前尚不清楚該方程對於合理的邊界條件和初始資料是否具有唯一解。迄今為止的數值方法似乎表明解是唯一的，但在沒有證明的情況下，我們提醒讀者，我們正在 編寫巨大的程式碼，這些程式碼應該為 Navier-Stokes 方程產生解，而我們真正研究的是演算法的輸出，我們希望這能告訴我們一些關於這些方程的資訊（這是對 Gallavoti^[9] 的改寫 - 在實踐中，儘管其數學基礎不確定，但這些程式碼似乎確實提供了關於許多（但並非所有）實際感興趣情況下的近乎不可壓縮流體運動的資訊。有關這些方程的這方面的更多資訊，請參見 Doering 和 Gibbon^[10]。）

我們再次考慮具有周期性邊界條件的模擬，以便於應用傅立葉變換。這也有助於透過使用 Orszag 和 Patterson^[11] 提出的一個想法（在 Canuto 等人^[12]^[13] 中也有解釋）來更容易地執行不可壓縮約束。如果我們對 Navier-Stokes 方程進行散度運算，我們將得到 $\nabla \cdot \left(\mathbf {u} \cdot \nabla \mathbf {u} \right)=-\Delta p$ ，因為 $\nabla \cdot \mathbf {u} =0$ 。因此 ${}p=-\Delta ^{-1}\left[\nabla \cdot \left(\mathbf {u} \cdot \nabla \mathbf {u} \right)\right]$ ，其中 $\Delta ^{-1}$ 是使用傅立葉變換定義的，因此如果 $f(x,y,z)$ 是一個平均值為零的週期性標量場，而 ${\hat {f}}$ 是它的傅立葉變換，那麼 ${\widehat {\Delta ^{-1}f}}=-{\frac {\hat {f}}{k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}$ ，其中 $k_{x}$ 、 $k_{y}$ 和 $k_{z}$ 是波數。然後 Navier-Stokes 方程變為

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

,

( 9)

只要初始資料滿足不可壓縮約束，就可以滿足不可壓縮約束。

為了在時間上離散化 9 ，我們將使用隱式中點規則。這將給出，

{\frac {\mathbf {u} ^{n+1}-\mathbf {u} ^{n}}{\delta t}}={}{\frac {0.5}{\text{Re}}}\Delta \left(\mathbf {u} ^{n+1}+\mathbf {u} ^{n}\right)-0.25\left(\mathbf {u} ^{n+1}+\mathbf {u} ^{n}\right)\cdot \nabla \left(\mathbf {u} ^{n+1}+\mathbf {u} ^{n}\right)

( 10)

$+0.25\nabla \left[\Delta ^{-1}\left(\nabla \cdot \left[\left(\mathbf {u} ^{n+1}+\mathbf {u} ^{n}\right)\cdot \nabla \left(\mathbf {u} ^{n+1}+\mathbf {u} ^{n}\right)\right]\right)\right]$ .

透過將程式與精確解進行比較，可以幫助測試程式的正確性。Shapiro^[14] 找到了以下精確解，它是氣象颶風模擬程式以及具有周期性邊界條件的 Navier-Stokes 求解器的良好測試

$u=-{\frac {A}{k^{2}+l^{2}}}\left[\lambda l\cos(kx)\sin(ly)\sin(mz)+mk\sin(kx)\cos(ly)\cos(mz)\right]\exp \left(-{\frac {\lambda ^{2}t}{\text{Re}}}\right)$

$v={\frac {A}{k^{2}+l^{2}}}\left[\lambda k\sin(kx)\cos(ly)\sin(mz)-ml\cos(kx)\sin(ly)\cos(mz)\right]\exp \left(-{\frac {\lambda ^{2}t}{\text{Re}}}\right)$

$w=A\cos(kx)\cos(ly)\sin(mz)\exp \left(-{\frac {\lambda ^{2}t}{\text{Re}}}\right)$

其中常數 $\lambda ={\sqrt {k^{2}+l^{2}+m^{2}}}$ 和 $l$ , $k$ 和 $m$ 是常數，選擇這些常數的限制條件是解在空間上是週期的。此類解的更多示例可以在 Majda 和 Bertozzi^[15] 中找到。

序列程式

首先，我們編寫 Matlab 程式來演示如何在單處理器上求解這些方程。第一個程式使用 Crank-Nicolson 時間步進法來求解二維 Navier-Stokes 方程，程式碼列在 A 中。為了測試程式，遵循 Laizet 和 Lamballais 的方法^[16]，我們使用在 $(x,y)\in [0,1]\times [0,1]$ 上的精確 Taylor-Green 渦旋解，並使用週期性邊界條件，表示式如下：

${}u(x,y,t)=\sin(2\pi x)\cos(2\pi y)\exp(-8\pi ^{2}\mu t)$

$v(x,y,t)=-\cos(2\pi x)\sin(2\pi y)\exp(-8\pi ^{2}\mu t)$ .

( A)

一個找到二維 Navier Stokes 方程數值解的 Matlab 程式程式碼下載

% Numerical solution of the 2D incompressible Navier-Stokes on a
% Square Domain [0,1]x[0,1] using a Fourier pseudo-spectral method
% and Crank-Nicolson timestepping. The numerical solution is compared to
% the exact Taylor-Green Vortex solution of the Navier-Stokes equations
%
%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)
clear all; format compact; format short; clc; clf;

Re=1;%Reynolds number

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

dt=0.0025; t(1)=0; tmax=.1;
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{\omega}^{n+1,k}
omegahat=fft2(omega);
%nonlinear term
nonlinhat=fft2(u.*ifft2(omegahat.*kx)+v.*ifft2(omegahat.*ky));
for i=1:nplots
    chg=1;
    % save old values
    uold=u; vold=v; omegaold=omega; omegacheck=omega;
    omegahatold=omegahat; nonlinhatold=nonlinhat;
    while chg>tol
        %nonlinear {n+1,k}
        nonlinhat=fft2(u.*ifft2(omegahat.*kx)+v.*ifft2(omegahat.*ky));
       
        %Crank–Nicolson timestepping
        omegahat=((1/dt + 0.5*(1/Re)*(kxx+kyy)).*omegahatold...
            -.5*(nonlinhatold+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}
        
        %\omega^{n+1,k+1}
        omega=ifft2(omegahat);
        % check for convergence
        chg=max(max(abs(omega-omegacheck)))
        % store omega to check for convergence of next iteration
        omegacheck=omega;
    end
     t(i+1)=t(i)+dt;
     uexact_y=-2*pi*sin(2*pi*xx).*sin(2*pi*yy).*exp(-8*pi^2*t(i+1)/Re);
     vexact_x=2*pi*sin(2*pi*xx).*sin(2*pi*yy).*exp(-8*pi^2*t(i+1)/Re);
     omegaexact=vexact_x-uexact_y;
     figure(1); pcolor(omega); xlabel x; ylabel y;
     title Numerical; colorbar; drawnow;
     figure(2); pcolor(omegaexact); xlabel x; ylabel y;
     title Exact; colorbar; drawnow;
     figure(3); pcolor(omega-omegaexact); xlabel x; ylabel y;
     title Error; colorbar; drawnow;
end

( Ap)

一個找到二維 Navier-Stokes 方程數值解的 Python 程式。程式碼下載

#!/usr/bin/env python
"""
Numerical solution of the 2D incompressible Navier-Stokes on a
Square Domain [0,1]x[0,1] using a Fourier pseudo-spectral method
and Crank-Nicolson timesteping. The numerical solution is compared to
the exact Taylor-Green Vortex solution of the Navier-Stokes equations

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*nu*t)
v(x,y,t)=-cos(2*pi*x)sin(2*pi*y)exp(-8*pi^2*nu*t)
"""

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=0.10
#nplots=int(tmax/dt)
Rey=1

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)
# 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]

src = mlab.imshow(xx,yy,omega,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)
omegahat=numpy.zeros((N,N), dtype=complex)
omegahatold=numpy.zeros((N,N), dtype=complex)
nlhat=numpy.zeros((N,N), dtype=complex)
nlhatold=numpy.zeros((N,N), dtype=complex)
dpsix=numpy.zeros((N,N), dtype=float)
dpsiy=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)
nlhat=numpy.fft.fft2(u*numpy.fft.ifft2(omegahat*kx)+\
v*numpy.fft.ifft2(omegahat*ky))
while (t<=tmax):
    chg=1.0

    # Save old values
    uold=u
    vold=v
    omegaold=omega
    omegacheck=omega
    omegahatold = omegahat
    nlhatold=nlhat
    
    while(chg>tol):
        # nolinear {n+1,k}
        nlhat=numpy.fft.fft2(u*numpy.fft.ifft2(omegahat*kx)+\
        v*numpy.fft.ifft2(omegahat*ky))

        # Crank-Nicolson timestepmath.ping
        omegahat=((1/dt + 0.5*(1/Rey)*(kxx+kyy))*omegahatold \
        -0.5*(nlhatold+nlhat)) \
        /(1/dt -0.5*(1/Rey)*(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))
        temp=abs(omega-omegacheck)
        chg=numpy.max(temp)
        print(chg)
        omegacheck=omega
    t+=dt
    src.mlab_source.scalars = omega
 
omegaexact=numpy.zeros((N,N), dtype=float)
for i in range(len(x)):
    for j in range(len(y)):
        uexact_y=-2*math.pi*numpy.sin(2*math.pi*x[i])*numpy.sin(2*math.pi*x[j])\
        *numpy.exp(-8*(math.pi**2)*t/Rey)
        vexact_x=2*math.pi*numpy.sin(2*math.pi*x[i])*numpy.sin(2*math.pi*y[j])\
        *numpy.exp(-8*(math.pi**2)*t/Rey)
        omegaexact[i][j]=vexact_x-uexact_y
numpy.savetxt('Error.txt',abs(omegaexact-omega))

第二個程式使用隱式中點規則對三維 Navier-Stokes 方程進行時間步進，程式碼列在 B 中。它也以 Taylor-Green 渦旋作為初始條件，因為該渦旋已被廣泛研究，因此提供了一個基線案例來比較結果。

( B)

一個找到三維 Navier Stokes 方程數值解的 Matlab 程式程式碼下載

% A program to solve the 3D Navier stokes 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 exact solution used to check the numerical method is in
% A. Shapiro "The use of an exact solution of the Navier-Stokes equations
% in a validation test of a three-dimensional nonhydrostatic numerical
% model" Monthly Weather Review vol. 121 pp. 2420-2425 (1993)

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 = 64; % number of harmonics
Ny = 64; % number of harmonics
Nz = 64; % number of harmonics
Nt = 10; % number of time slices
dt = 0.2/Nt; % time step
t=0; % initial time
Re = 1.0; % 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);

% exact solution
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);

% calculate vorticity for plotting
omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy;
omegatot=omegax.^2+omegay.^2+omegaz.^2;
figure(1); clf; n=0;
subplot(2,2,1); title(['omega x ',num2str(n*dt)]);
p1 = patch(isosurface(x,y,z,omegax,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(x,y,z,omegax,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegax,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); colorbar;
subplot(2,2,2); title(['omega y ',num2str(n*dt)]);
p1 = patch(isosurface(x,y,z,omegay,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(x,y,z,omegay,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegay,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); colorbar;
subplot(2,2,3); title(['omega z ',num2str(n*dt)]);
p1 = patch(isosurface(x,y,z,omegaz,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(x,y,z,omegaz,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegaz,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); colorbar;
subplot(2,2,4); title(['|omega|^2 ',num2str(n*dt)]);
p1 = patch(isosurface(x,y,z,omegatot,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
p2 = patch(isocaps(x,y,z,omegatot,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegatot,p1); lighting phong;
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;
    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;
    chg=1; t=t+dt;
    while (chg>tol)
        nonlinu=0.25*((u+uold).*(ux+uxold)...
                      +(v+vold).*(uy+uyold)...
                      +(w+wold).*(uz+uzold));
        nonlinv=0.25*((u+uold).*(vx+vxold)...
                      +(v+vold).*(vy+vyold)...
                      +(w+wold).*(vz+vzold));
        nonlinw=0.25*((u+uold).*(wx+wxold)...
                      +(v+vold).*(wy+wyold)...
                      +(w+wold).*(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);
        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);
        chg=max(abs(utemp-u))+max(abs(vtemp-v))+max(abs(wtemp-w));
    end
    % calculate vorticity for plotting
    omegax=wy-vz; omegay=uz-wx; omegaz=vx-uy;
    omegatot=omegax.^2+omegay.^2+omegaz.^2;
    figure(1); clf;
    subplot(2,2,1); title(['omega x ',num2str(t)]);
    p1 = patch(isosurface(x,y,z,omegax,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(x,y,z,omegax,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegax,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z');
    axis equal; axis square; view(3); colorbar;
    subplot(2,2,2); title(['omega y ',num2str(t)]);
    p1 = patch(isosurface(x,y,z,omegay,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(x,y,z,omegay,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegay,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z');
    axis equal; axis square; view(3); colorbar;
    subplot(2,2,3); title(['omega z ',num2str(t)]);
    p1 = patch(isosurface(x,y,z,omegaz,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(x,y,z,omegaz,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegaz,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z');
    axis equal; axis square; view(3); colorbar;
    subplot(2,2,4); title(['|omega|^2 ',num2str(t)]);
    p1 = patch(isosurface(x,y,z,omegatot,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(x,y,z,omegatot,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegatot,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
    axis equal; axis square; view(3);
end
toc

uexact=-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);
        
vexact=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);
      
wexact=cos(sk*xx).*cos(sl*yy).*sin(sm*zz)*exp(-t*(lamlkm^2)/Re);


error= max(max(max(abs(u-uexact))))+...
        max(max(max(abs(v-vexact))))+...
        max(max(max(abs(w-wexact))))

( Bp)

一個找到三維 Navier-Stokes 方程數值解的 Python 程式。程式碼下載

#!/usr/bin/env python
"""
A program to solve the 3D Navier Stokes 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 exact solution used to check the numerical method is in
A. Shapiro "The use of an exact solution of the Navier-Stokes equations
in a validation test of a three-dimensional nonhydrostatic numerical
model" Monthly Weather Review vol. 121 pp. 2420-2425 (1993)

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
Rey=1.0      	# 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)

uexact=numpy.zeros((Nx,Ny,Nz), dtype=float)
vexact=numpy.zeros((Nx,Ny,Nz), dtype=float)
wexact=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)

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)

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)

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)
phat=numpy.zeros((Nx,Ny,Nz), dtype=complex)
temphat=numpy.zeros((Nx,Ny,Nz), dtype=complex)

rhsuhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)
rhswhatfix=numpy.zeros((Nx,Ny,Nz), dtype=complex)
rhsvhatfix=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)

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 vorticity for plotting

omegax=wy-vz
omegay=uz-wx
omegaz=vx-uy
omegatot=omegax**2.0 + omegay**2.0 + omegaz**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)
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
	rhsuhatfix=(1.0/dt + (0.5/Rey)*(k2xm+k2ym+k2zm))*uhat
	rhsvhatfix=(1.0/dt + (0.5/Rey)*(k2xm+k2ym+k2zm))*vhat
	rhswhatfix=(1.0/dt + (0.5/Rey)*(k2xm+k2ym+k2zm))*what
	chg=1.0
	t=t+dt
	while(chg>tol):
		nonlinu=0.25*((u+uold)*(ux+uxold)+(v+vold)*(uy+uyold)+(w+wold)*(uz+uzold))
		nonlinv=0.25*((u+uold)*(vx+vxold)+(v+vold)*(vy+vyold)+(w+wold)*(vz+vzold))
		nonlinw=0.25*((u+uold)*(wx+wxold)+(v+vold)*(wy+wyold)+(w+wold)*(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)
		uhat=(rhsuhatfix-nonlinuhat-kxm*phat)/(1.0/dt - (0.5/Rey)*(k2xm+k2ym+k2zm))
		vhat=(rhsvhatfix-nonlinvhat-kym*phat)/(1.0/dt - (0.5/Rey)*(k2xm+k2ym+k2zm))
		what=(rhswhatfix-nonlinwhat-kzm*phat)/(1.0/dt - (0.5/Rey)*(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))

		chg=numpy.max(abs(u-utemp))+numpy.max(abs(v-vtemp))+numpy.max(abs(w-wtemp))
		
	# 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

       
uexact=-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)
vexact= 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)
wexact= numpy.cos(sk*xx)*numpy.cos(sl*yy)*numpy.sin(sm*zz)*numpy.exp(-t*(lamlkm**2.0)/Rey)

err=numpy.max(abs(u-uexact))+numpy.max(abs(v-vexact))+numpy.max(abs(w-wexact))
print(err)

練習

證明對於 Taylor-Green 渦旋解，二維 Navier-Stokes 方程中的非線性項完全抵消。
編寫一個 Matlab 程式，使用隱式中點規則而不是 Crank-Nicolson 方法來獲得二維 Navier-Stokes 方程的解。將您的數值解與 Taylor-Green 渦旋解進行比較。
編寫一個 Fortran 程式，使用隱式中點規則而不是 Crank-Nicolson 方法來獲得二維 Navier-Stokes 方程的解。將您的數值解與 Taylor-Green 渦旋解進行比較。
編寫一個 Matlab 程式，使用 Crank-Nicolson 方法而不是隱式中點規則來獲得三維 Navier-Stokes 方程的解。
編寫一個 Fortran 程式，使用 Crank-Nicolson 方法而不是隱式中點規則來獲得三維 Navier-Stokes 方程的解。
如方程 3 和 4 所示的 Navier-Stokes 方程還滿足進一步的積分性質。特別地，證明
1. ${\frac {\rho }{2}}{\frac {\mathrm {d} }{\mathrm {d} t}}\lVert \omega \rVert _{l^{2}}^{2}=-\mu \lVert \nabla \omega \rVert _{l^{2}}^{2},$ 其中 $\lVert \omega \rVert _{l^{2}}^{2}=\int \int (\omega )^{2}\mathrm {d} x\mathrm {d} y$ 和 $\lVert \nabla \omega \rVert _{l^{2}}^{2}=\int \int (\nabla \omega )\cdot (\nabla \omega )\mathrm {d} x\mathrm {d} y.$ 提示：將方程 3 乘以 $\omega$ ，然後用分部積分。
2. 證明 (a) 意味著 $\lVert \omega (t=T)\rVert _{l^{2}}^{2}-\lVert \omega (t=0)\rVert _{l^{2}}^{2}=-\mu \int _{0}^{T}\lVert \nabla \omega \rVert _{l^{2}}^{2}\mathrm {d} t$
3. 部分 (b) 給出了一個可以在積分二維 Navier-Stokes 方程時檢查的性質。現在我們證明隱式中點規則滿足一個類似的性質。將方程 7 乘以 $0.5(\omega ^{n+1}+\omega ^{n})$ , 在空間上進行分部積分，然後對時間求和以推斷 $\lVert \omega ^{N}\rVert _{l^{2}}^{2}-\lVert \omega ^{0}\rVert _{l^{2}}^{2}=-{\frac {\mu }{4}}\sum _{n=0}^{N-1}\lVert \nabla (\omega ^{n}+\omega ^{n+1})\rVert _{l^{2}}^{2}\delta t$ .
4. 推斷這意味著隱式中點規則時間步進方法是無條件穩定的，前提是非線性項可以求解^[17]。

並行程式：OpenMP

我們不會提供完全並行的示例程式，而是使用 Fortran 提供一個簡單的實現，用於二維和三維 Navier-Stokes 方程的 Crank-Nicolson 和隱式中點規則演算法，這些演算法在 Matlab 中已經介紹過。二維方程的程式在列表 C 中，用於繪製所得渦量場的示例 Matlab 指令碼在列表 D 中。該程式在列表 E 中，用於繪製所得渦量場的示例 Matlab 指令碼在列表 F 中。

( C)

一個用於求解二維 Navier-Stokes 方程的 Fortran 程式程式碼下載

PROGRAM main
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program numerically solves the 2D incompressible Navier-Stokes
! on a Square Domain [0,1]x[0,1] using pseudo-spectral methods and
! Crank-Nicolson timestepping. The numerical solution is compared to
! the exact Taylor-Green Vortex Solution.
!
! 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*nu*t)
! v(x,y,t)=-cos(2*pi*x)sin(2*pi*y)exp(-8*pi^2*nu*t)
!
! .. 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
! 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
! .. 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=256	
INTEGER(kind=4), PARAMETER :: Ny=256	
REAL(kind=8), PARAMETER	:: dt=0.00125	
REAL(kind=8), PARAMETER	&
:: pi=3.14159265358979323846264338327950288419716939937510
REAL(kind=8), PARAMETER	:: rho=1.0d0
REAL(kind=8), PARAMETER	:: mu=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, omegoldhat_x, omegold_x,&	
omegoldhat_y, omegold_y, nlold, nloldhat,&
omeghat, omeghat_x, omeghat_y, omeg_x, omeg_y,&
nl, nlhat, psihat, psihat_x, psi_x, psihat_y, psi_y
REAL(kind=8),DIMENSION(:,:), ALLOCATABLE	::	uexact_y,vexact_x,omegexact
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),omegoldhat_x(1:Nx,1:Ny), omegold_x(1:Nx,1:Ny), &
omegoldhat_y(1:Nx,1:Ny),omegold_y(1:Nx,1:Ny), nlold(1:Nx,1:Ny), nloldhat(1:Nx,1:Ny),&
omeghat(1:Nx,1:Ny), omeghat_x(1:Nx,1:Ny), omeghat_y(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), &
psihat_x(1:Nx,1:Ny), psi_x(1:Nx,1:Ny), psihat_y(1:Nx,1:Ny), psi_y(1:Nx,1:Ny),&
uexact_y(1:Nx,1:Ny), vexact_x(1:Nx,1:Ny), omegexact(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)
END DO
END DO

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

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!Initial nonlinear term !!!!!!!!!!!!!!!!!!!!!!!!
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! obtain \hat{\omega}_x^{n,k}
DO j=1,Ny
omeghat_x(1:Nx,j)=omeghat(1:Nx,j)*kx(1:Nx)
END DO
! obtain \hat{\omega}_y^{n,k}
DO i=1,Nx
omeghat_y(i,1:Ny)=omeghat(i,1:Ny)*ky(1:Ny)
END DO
! convert to real space
CALL dfftw_execute_dft_(planbxy,omeghat_x(1:Nx,1:Ny),omeg_x(1:Nx,1:Ny))
CALL dfftw_execute_dft_(planbxy,omeghat_y(1:Nx,1:Ny),omeg_y(1:Nx,1:Ny))
! compute nonlinear term in real space
DO j=1,Ny
nl(1:Nx,j)=u(1:Nx,j)*omeg_x(1:Nx,j)/REAL(Nx*Ny,kind(0d0))+&
v(1:Nx,j)*omeg_y(1:Nx,j)/REAL(Nx*Ny,kind(0d0))	
END DO
CALL dfftw_execute_dft_(planfxy,nl(1:Nx,1:Ny),nlhat(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)
nloldhat(1:Nx,1:Ny)=nlhat(1:Nx,1:Ny)
DO WHILE (chg>tol)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!nonlinear fixed (n,k+1)!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! obtain \hat{\omega}_x^{n+1,k}
DO j=1,Ny
omeghat_x(1:Nx,j)=omeghat(1:Nx,j)*kx(1:Nx)
END DO
! obtain \hat{\omega}_y^{n+1,k}
DO i=1,Nx
omeghat_y(i,1:Ny)=omeghat(i,1:Ny)*ky(1:Ny)
END DO
! convert back to real space
CALL dfftw_execute_dft_(planbxy,omeghat_x(1:Nx,1:Ny),omeg_x(1:Nx,1:Ny))
CALL dfftw_execute_dft_(planbxy,omeghat_y(1:Nx,1:Ny),omeg_y(1:Nx,1:Ny))
! calculate nonlinear term in real space
DO j=1,Ny
nl(1:Nx,j)=u(1:Nx,j)*omeg_x(1:Nx,j)/REAL(Nx*Ny,kind(0d0))+&
v(1:Nx,j)*omeg_y(1:Nx,j)/REAL(Nx*Ny,kind(0d0))
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 Crank Nicolson timestepping
DO j=1,Ny
omeghat(1:Nx,j)=( (1.0d0/dt+0.5d0*(mu/rho)*(kxx(1:Nx)+kyy(j)))&
*omegoldhat(1:Nx,j) - 0.5d0*(nloldhat(1:Nx,j)+nlhat(1:Nx,j)))/&
(1.0d0/dt-0.5d0*(mu/rho)*(kxx(1:Nx)+kyy(j)))
END DO

! calculate \hat{\psi}^{n+1,k+1}
DO j=1,Ny
psihat(1:Nx,j)=-omeghat(1:Nx,j)/(kxx(1:Nx)+kyy(j))
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 j=1,Ny
psihat_x(1:Nx,j)=psihat(1:Nx,j)*kx(1:Nx)
END DO
CALL dfftw_execute_dft_(planbxy,psihat_x(1:Nx,1:Ny),psi_x(1:Nx,1:Ny))
DO i=1,Nx
psihat_y(i,1:Ny)=psihat(i,1:Ny)*ky(1:Ny)
END DO
CALL dfftw_execute_dft_(planbxy,psihat_y(1:Ny,1:Ny),psi_y(1:Ny,1:Ny))
DO j=1,Ny
psi_x(1:Nx,j)=psi_x(1:Nx,j)/REAL(Nx*Ny,kind(0d0))
psi_y(1:Nx,j)=psi_y(1:Nx,j)/REAL(Nx*Ny,kind(0d0))
END DO	

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

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

! check for convergence
chg=maxval(abs(omeg-omegcheck))
! saves {n+1,k+1} to {n,k} for next iteration
omegcheck=omeg
END DO
time(n+1)=time(n)+dt
PRINT *,'TIME ',time(n+1)
END DO
DO j=1,Ny
DO i=1,Nx
uexact_y(i,j)=-2.0d0*pi*sin(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j))*&
exp(-8.0d0*mu*(pi**2)*nplots*dt)
vexact_x(i,j)=2.0d0*pi*sin(2.0d0*pi*x(i))*sin(2.0d0*pi*y(j))*&
exp(-8.0d0*mu*(pi**2)*nplots*dt)
omegexact(i,j)=vexact_x(i,j)-uexact_y(i,j)
END DO
END DO
name_config = 'omegafinal.datbin'
INQUIRE(iolength=iol) omegexact(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 = 'omegaexactfinal.datbin'
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) omegexact(i,j)
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, omegoldhat_x, omegold_x,&	
omegoldhat_y, omegold_y, nlold, nloldhat,&
omeghat, omeghat_x, omeghat_y, omeg_x, omeg_y,&
nl, nlhat, psihat, psihat_x, psi_x, psihat_y, psi_y,&
uexact_y,vexact_x,omegexact, &
fftfx,fftbx,stat=AllocateStatus)	
IF (AllocateStatus .ne. 0) STOP
PRINT *,'Program execution complete'
END PROGRAM main

( D)

一個 Matlab 程式，用於繪製列表 C 生成的渦量場和誤差程式碼下載

% A program to create a plot of the computed results
% from the 2D Matlab Navier-Stokes solver

clear all; format compact, format short,
set(0,'defaultaxesfontsize',14,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',2,'defaultpatchlinewidth',3.5);

% Load data
% Get coordinates
X=load('xcoord.dat');
Y=load('ycoord.dat');
% find number of grid points
Nx=length(X);
Ny=length(Y);

% reshape coordinates to allow easy plotting
[xx,yy]=ndgrid(X,Y);

%
% Open file and dataset using the default properties.
%
FILENUM=['omegafinal.datbin'];
FILEEXA=['omegaexactfinal.datbin'];
fidnum=fopen(FILENUM,'r');
[fnamenum,modenum,mformatnum]=fopen(fidnum);
fidexa=fopen(FILEEXA,'r');
[fnameexa,modeexa,mformatexa]=fopen(fidexa);
Num=fread(fidnum,Nx*Ny,'double',mformatnum);
Exa=fread(fidexa,Nx*Ny,'double',mformatexa);
Num=reshape(Num,Nx,Ny);
Exa=reshape(Exa,Nx,Ny);
% close files
fclose(fidnum);
fclose(fidexa);
%
% Plot data on the screen.
%
figure(2);clf;
subplot(3,1,1); contourf(xx,yy,Num);
title(['Numerical Solution ']);
colorbar; axis square;
subplot(3,1,2); contourf(xx,yy,Exa);
title(['Exact Solution ']);
colorbar; axis square;
subplot(3,1,3); contourf(xx,yy,Exa-Num);
title(['Error']);
colorbar; axis square;
drawnow;

( E)

一個用於求解三維 Navier-Stokes 方程的 Fortran 程式程式碼下載

PROGRAM main	
!-----------------------------------------------------------------------------------
!
!
! PURPOSE
!
! This program numerically solves the 3D incompressible Navier-Stokes
! on a Cubic Domain [0,2pi]x[0,2pi]x[0,2pi] using pseudo-spectral methods and
! Implicit Midpoint rule timestepping. The numerical solution is compared to
! an exact solution reported by Shapiro
!
! 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
! 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
! .. 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
!
! A. Shapiro " The use of an exact solution of the Navier-Stokes equations
! in a validation test of a three-dimensional nonhydrostatic numerical model"
! Monthly Weather Review vol. 121, 2420-2425, (1993).
!
! 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/)
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=20
REAL(kind=8), PARAMETER	:: dt=0.2d0/Nt
REAL(kind=8), PARAMETER	:: Re=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	:: dtInv=1.0d0/REAL(dt,kind(0d0))
REAL(kind=8)	:: scalemodes,chg,factor
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x, y, z, time
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

COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx, ky, kz	
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE	:: uhat, vhat, what,&
rhsuhatfix, rhsvhatfix,&
rhswhatfix, nonlinuhat,&
nonlinvhat, nonlinwhat,&
phat,temp_c
REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: realtemp
!FFTW variables
INTEGER(kind=4)	:: ierr
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	
INTEGER(kind=8)	:: planfxyz,planbxyz

!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

    PRINT *,'Grid:',Nx,'X',Ny,'Y',Nz,'Z'
PRINT *,'dt:',dt	
ALLOCATE(x(1:Nx),y(1:Ny),z(1:Nz),time(1:Nt+1),u(1:Nx,1:Ny,1:Nz),&
  v(1:Nx,1:Ny,1:Nz), w(1:Nx,1:Ny,1:Nz), ux(1:Nx,1:Ny,1:Nz),&
uy(1:Nx,1:Ny,1:Nz), uz(1:Nx,1:Ny,1:Nz), vx(1:Nx,1:Ny,1:Nz),&
vy(1:Nx,1:Ny,1:Nz), vz(1:Nx,1:Ny,1:Nz), wx(1:Nx,1:Ny,1:Nz),&
wy(1:Nx,1:Ny,1:Nz), wz(1:Nx,1:Ny,1:Nz), uold(1:Nx,1:Ny,1:Nz),&
uxold(1:Nx,1:Ny,1:Nz), uyold(1:Nx,1:Ny,1:Nz), uzold(1:Nx,1:Ny,1:Nz),&
vold(1:Nx,1:Ny,1:Nz), vxold(1:Nx,1:Ny,1:Nz), vyold(1:Nx,1:Ny,1:Nz),&
vzold(1:Nx,1:Ny,1:Nz), wold(1:Nx,1:Ny,1:Nz), wxold(1:Nx,1:Ny,1:Nz),&
wyold(1:Nx,1:Ny,1:Nz), wzold(1:Nx,1:Ny,1:Nz), utemp(1:Nx,1:Ny,1:Nz),&
  vtemp(1:Nx,1:Ny,1:Nz), wtemp(1:Nx,1:Ny,1:Nz), temp_r(1:Nx,1:Ny,1:Nz),&
kx(1:Nx),ky(1:Ny),kz(1:Nz),uhat(1:Nx,1:Ny,1:Nz), vhat(1:Nx,1:Ny,1:Nz),&
what(1:Nx,1:Ny,1:Nz), rhsuhatfix(1:Nx,1:Ny,1:Nz),&
  rhsvhatfix(1:Nx,1:Ny,1:Nz), rhswhatfix(1:Nx,1:Ny,1:Nz),&
  nonlinuhat(1:Nx,1:Ny,1:Nz), nonlinvhat(1:Nx,1:Ny,1:Nz),&
  nonlinwhat(1:Nx,1:Ny,1:Nz), phat(1:Nx,1:Ny,1:Nz),temp_c(1:Nx,1:Ny,1:Nz),&
realtemp(1:Nx,1:Ny,1:Nz), stat=AllocateStatus)	
IF (AllocateStatus .ne. 0) STOP
PRINT *,'allocated space'

CALL dfftw_plan_dft_3d_(planfxyz,Nx,Ny,Nz,temp_r(1:Nx,1:Ny,1:Nz),&
temp_c(1:Nx,1:Ny,1:Nz),FFTW_FORWARD,FFTW_ESTIMATE)
CALL dfftw_plan_dft_3d_(planbxyz,Nx,Ny,Nz,temp_c(1:Nx,1:Ny,1:Nz),&
temp_r(1:Nx,1:Ny,1:Nz),FFTW_BACKWARD,FFTW_ESTIMATE)
PRINT *,'Setup 3D FFTs'

! 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))
PRINT *,'Setup grid and fourier frequencies'

!initial conditions for Taylor-Green vortex
! factor=2.0d0/sqrt(3.0d0)
! DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
! 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=1,Nz; DO j=1,Ny; DO i=1,Nx
! 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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
! w(i,j,k)=factor*sin(theta)*cos(x(i))*cos(y(j))*sin(z(k))
! END DO ; END DO ; END DO

! Initial conditions for exact solution
time(1)=0.0d0
factor=sqrt(3.0d0)
DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
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=1,Nz; DO j=1,Ny; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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 dfftw_execute_dft_(planfxyz,u(1:Nx,1:Ny,1:Nz),uhat(1:Nx,1:Ny,1:Nz))
CALL dfftw_execute_dft_(planfxyz,v(1:Nx,1:Ny,1:Nz),vhat(1:Nx,1:Ny,1:Nz))
CALL dfftw_execute_dft_(planfxyz,w(1:Nx,1:Ny,1:Nz),what(1:Nx,1:Ny,1:Nz))

! derivative of u with respect to x, y, and z
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=uhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),ux(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=uhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),uy(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),uz(1:Nx,1:Ny,1:Nz))

! derivative of v with respect to x, y, and z
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=vhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),vx(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=vhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),vy(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=vhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),vz(1:Nx,1:Ny,1:Nz))

! derivative of w with respect to x, y, and z
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=what(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),wx(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=what(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),wy(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=what(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),wz(1:Nx,1:Ny,1:Nz))
! save initial data
time(1)=0.0
n=0
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)
!omegay
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)
!omegaz
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)

DO n=1,Nt
!fixed point
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
rhsuhatfix(i,j,k) = (dtInv+(0.5d0*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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
rhsvhatfix(i,j,k) = (dtInv+(0.5d0*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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
rhswhatfix(i,j,k) = (dtInv+(0.5d0*ReInv)*&
(kx(i)*kx(i)+ky(j)*ky(j)+kz(k)*kz(k)))*what(i,j,k)
END DO ; END DO ; END DO
chg=1
DO WHILE (chg .gt. tol)
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)))
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planfxyz,temp_r(1:Nx,1:Ny,1:Nz),nonlinuhat(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)))
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planfxyz,temp_r(1:Nx,1:Ny,1:Nz),nonlinvhat(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)))
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planfxyz,temp_r(1:Nx,1:Ny,1:Nz),nonlinwhat(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=uhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),ux(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=uhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),uy(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=uhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),uz(1:Nx,1:Ny,1:Nz))

! derivative of v with respect to x, y, and z
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=vhat(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),vx(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=vhat(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),vy(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=vhat(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),vz(1:Nx,1:Ny,1:Nz))

! derivative of w with respect to x, y, and z
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=what(i,j,k)*kx(i)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),wx(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=what(i,j,k)*ky(j)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),wy(1:Nx,1:Ny,1:Nz))
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
temp_c(i,j,k)=what(i,j,k)*kz(k)*scalemodes
END DO ; END DO ; END DO
CALL dfftw_execute_dft_(planbxyz,temp_c(1:Nx,1:Ny,1:Nz),wz(1:Nx,1:Ny,1:Nz))

DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
utemp(i,j,k)=u(i,j,k)
END DO ; END DO ; END DO
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
vtemp(i,j,k)=v(i,j,k)
END DO ; END DO ; END DO
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
wtemp(i,j,k)=w(i,j,k)
END DO ; END DO ; END DO

CALL dfftw_execute_dft_(planbxyz,uhat(1:Nx,1:Ny,1:Nz),u(1:Nx,1:Ny,1:Nz))
CALL dfftw_execute_dft_(planbxyz,vhat(1:Nx,1:Ny,1:Nz),v(1:Nx,1:Ny,1:Nz))
CALL dfftw_execute_dft_(planbxyz,what(1:Nx,1:Ny,1:Nz),w(1:Nx,1:Ny,1:Nz))

DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
u(i,j,k)=u(i,j,k)*scalemodes
END DO ; END DO ; END DO
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
v(i,j,k)=v(i,j,k)*scalemodes
END DO ; END DO ; END DO
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
w(i,j,k)=w(i,j,k)*scalemodes
END DO ; END DO ; END DO
chg =maxval(abs(utemp-u))+maxval(abs(vtemp-v))+maxval(abs(wtemp-w))
PRINT *,'chg:',chg
END DO
time(n+1)=n*dt
PRINT *,'time',n*dt
!NOTE: utemp, vtemp, and wtemp are just temporary space that can be used
! instead of creating new arrays.
!omegax
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)
!omegay
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)
!omegaz
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
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)
END DO

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'

! Calculate error in final numerical solution
DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
utemp(i,j,k)=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(Nt+1)/Re))
END DO; END DO; END DO
DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
vtemp(i,j,k)=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(Nt+1)/Re))
END DO ; END DO ; END DO
DO k=1,Nz ; DO j=1,Ny ; DO i=1,Nx
wtemp(i,j,k)=w(i,j,k)-&
(cos(x(i))*cos(y(j))*sin(z(k))*exp(-(factor**2)*time(Nt+1)/Re))
END DO ; END DO ; END DO
chg=maxval(abs(utemp))+maxval(abs(vtemp))+maxval(abs(wtemp))
PRINT*,'The error at the final timestep is',chg

CALL dfftw_destroy_plan_(planfxyz)
CALL dfftw_destroy_plan_(planbxyz)
DEALLOCATE(x,y,z,time,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,kx,ky,kz,uhat,vhat,what,rhsuhatfix,rhsvhatfix,&
  rhswhatfix,phat,nonlinuhat,nonlinvhat,nonlinwhat,temp_c,&
realtemp,stat=AllocateStatus)	
IF (AllocateStatus .ne. 0) STOP
PRINT *,'Program execution complete'
END PROGRAM main

( F)

一個 Matlab 程式，用於繪製列表 E 生成的渦量場程式碼下載

% A program to create a plot of the computed results
% from the 3D Fortran Navier-Stokes solver
clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
    'defaultaxesfontweight','bold')
% Load data
% Get coordinates
tdata=load('./data/tdata.dat');
x=load('./data/xcoord.dat');
y=load('./data/ycoord.dat');
z=load('./data/zcoord.dat');
nplots = length(tdata);

Nx = length(x); Nt = length(tdata);
Ny = length(y); Nz = length(z);

% reshape coordinates to allow easy plotting
[yy,xx,zz]=meshgrid(x,y,z);

for i =1:nplots
    %
    % Open file and dataset using the default properties.
    %
    FILEX=['./data/omegax',num2str(9999999+i),'.datbin'];
    FILEY=['./data/omegay',num2str(9999999+i),'.datbin'];
    FILEZ=['./data/omegaz',num2str(9999999+i),'.datbin'];
    FILEPIC=['./data/pic',num2str(9999999+i),'.jpg'];
    fid=fopen(FILEX,'r');
    [fname,mode,mformat]=fopen(fid);
    omegax=fread(fid,Nx*Ny*Nz,'real*8');
    omegax=reshape(omegax,Nx,Ny,Nz);
    fclose(fid);
    fid=fopen(FILEY,'r');
    [fname,mode,mformat]=fopen(fid);
    omegay=fread(fid,Nx*Ny*Nz,'real*8');
    omegay=reshape(omegay,Nx,Ny,Nz);
    fclose(fid);
    fid=fopen(FILEZ,'r');
    [fname,mode,mformat]=fopen(fid);
    omegaz=fread(fid,Nx*Ny*Nz,'real*8');
    omegaz=reshape(omegaz,Nx,Ny,Nz);
    fclose(fid);
    %
    % Plot data on the screen.
    %
    omegatot=omegax.^2+omegay.^2+omegaz.^2;
    figure(100); clf;
    subplot(2,2,1); title(['omega x ',num2str(tdata(i))]);
    p1 = patch(isosurface(xx,yy,zz,omegax,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(xx,yy,zz,omegax,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegax,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z');
    axis equal; axis square; view(3); colorbar;
    subplot(2,2,2); title(['omega y ',num2str(tdata(i))]);
    p1 = patch(isosurface(xx,yy,zz,omegay,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(xx,yy,zz,omegay,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegay,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z');
    axis equal; axis square; view(3); colorbar;
    subplot(2,2,3); title(['omega z ',num2str(tdata(i))]);
    p1 = patch(isosurface(xx,yy,zz,omegaz,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(xx,yy,zz,omegaz,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegaz,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z');
    axis equal; axis square; view(3); colorbar;
    subplot(2,2,4); title(['|omega|^2 ',num2str(tdata(i))]);
    p1 = patch(isosurface(xx,yy,zz,omegatot,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.3);
    p2 = patch(isocaps(xx,yy,zz,omegatot,.0025),...
            'FaceColor','interp','EdgeColor','none','FaceAlpha',0.1);
        isonormals(omegatot,p1); lighting phong;
    xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
    axis equal; axis square; view(3);
saveas(100,FILEPIC);

end

練習

驗證列表 C 中的程式在時間上是二階精度的。
使用 OpenMP 指令將二維 Navier Stokes 方程的示例 Fortran 程式碼並行化。儘可能提高效率。
編寫另一個程式碼，使用執行緒化的 FFTW 來進行快速傅立葉變換。該程式碼的結構應該與 https://wikibook.tw/wiki/Parallel_Spectral_Numerical_Methods/The_Cubic_Nonlinear_Schrodinger_Equation#equation_G 中的程式列表類似。
使用 OpenMP 指令將列表 E 中的三維 Navier-Stokes 方程的示例 Fortran 程式碼並行化。儘可能提高效率。
編寫另一個程式碼，使用執行緒化的 FFTW 來進行三維 Navier-Stokes 方程的快速傅立葉變換。該程式碼的結構應該與 https://wikibook.tw/wiki/Parallel_Spectral_Numerical_Methods/The_Cubic_Nonlinear_Schrodinger_Equation#equation_J 中的程式類似。

並行程式：MPI

在本節中，我們將使用從 http://www.2decomp.org/index.html 獲取的 2DECOMP&FFT 庫。該網站包含一些示例，說明了如何使用此庫，特別是 http://www.2decomp.org/case_study1.html 中的示例程式碼，這是一個非常有用的說明，說明如何將使用 FFTW 的程式碼轉換為使用 MPI 和上述庫的程式碼。

此程式碼與列表 C 中的序列程式碼非常相似。為了完整性和為了讓您瞭解如何並行化其他程式，我們將其包含在內。該程式使用 2DECOMP&FFT 庫。此程式和序列程式之間的一個區別是，包含一個子程式來寫入資料。由於計算的這部分被重複多次，因此將重複程式碼放在子程式中會使程式更具可讀性。該子程式也足夠通用，可以重複用於其他程式，從而節省了程式開發人員的時間。

( G)

一個用於求解三維 Navier-Stokes 方程的並行 MPI Fortran 程式程式碼下載

PROGRAM main	
!-----------------------------------------------------------------------------------
!
!
! PURPOSE
!
! This program numerically solves the 3D incompressible Navier-Stokes
! on a Cubic Domain [0,2pi]x[0,2pi]x[0,2pi] using pseudo-spectral methods and
! Implicit Midpoint rule timestepping. The numerical solution is compared to
! an exact solution reported by Shapiro
!
! 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 u to check convergence
! wtemp = temporary storage of u 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
! 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
!
! A. Shapiro " The use of an exact solution of the Navier-Stokes equations
! in a validation test of a three-dimensional nonhydrostatic numerical model"
! Monthly Weather Review vol. 121, 2420-2425, (1993).
!
! 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=256
INTEGER(kind=4), PARAMETER :: Ny=256
INTEGER(kind=4), PARAMETER :: Nz=256
    INTEGER(kind=4), PARAMETER :: Lx=1
INTEGER(kind=4), PARAMETER :: Ly=1
INTEGER(kind=4), PARAMETER :: Lz=1
INTEGER(kind=4), PARAMETER	:: Nt=20
REAL(kind=8), PARAMETER	:: dt=0.05d0/Nt
REAL(kind=8), PARAMETER	:: Re=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	:: 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

COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx, ky, kz	
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE	:: uhat, vhat, what,&
rhsuhatfix, rhsvhatfix,&
rhswhatfix, 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:3),allchg(1:3),&
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)),&
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)),&
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)),&
  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 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)
        
!start timer
CALL system_clock(start,count_rate)
DO n=1,Nt
!fixed point
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
chg=1
DO WHILE (chg .gt. tol)
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)))
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)))
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)))
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
mychg(1) =maxval(abs(utemp-u))
mychg(2) =maxval(abs(vtemp-v))
mychg(3) =maxval(abs(wtemp-w))
CALL MPI_ALLREDUCE(mychg,allchg,3,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,ierr)
chg=allchg(1)+allchg(2)+allchg(3)
IF (myid.eq.0) THEN
PRINT *,'chg:',chg
END IF
END DO
time(n+1)=n*dt

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

                !save omegax, omegay, and omegaz
!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)
                !5100 continue
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

! Calculate error in final numerical solution
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) -&
(-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(Nt+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)
vtemp(i,j,k)=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(Nt+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)
wtemp(i,j,k)=w(i,j,k)-&
(cos(x(i))*cos(y(j))*sin(z(k))*exp(-(factor**2)*time(Nt+1)/Re))
END DO ; END DO ; END DO
mychg(1) = maxval(abs(utemp))
mychg(2) = maxval(abs(vtemp))
mychg(3) = maxval(abs(wtemp))
CALL MPI_ALLREDUCE(mychg,allchg,3,MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,ierr)
chg=allchg(1)+allchg(2)+allchg(3)
IF (myid.eq.0) THEN
PRINT*,'The error at the final timestep is',chg
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,kx,ky,kz,uhat,vhat,what,rhsuhatfix,rhsvhatfix,&
  rhswhatfix,phat,nonlinuhat,nonlinvhat,nonlinwhat,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

( H)

一個用於儲存用於求解三維 Navier-Stokes 方程的並行 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

( I)

一個用於編譯用於求解三維 Navier-Stokes 方程的並行 MPI Fortran 程式的 makefile 程式碼下載

COMPILER = mpif90
decompdir= ../2decomp_fft
FLAGS = -O0

DECOMPLIB = -I${decompdir}/include -L${decompdir}/lib -l2decomp_fft
LIBS = #-L${FFTW_LINK} -lfftw3 -lm
SOURCES = NavierStokes3DfftIMR.f90 savedata.f90

ns3d: $(SOURCES)
${COMPILER} -o ns3d $(FLAGS) $(SOURCES) $(LIBS) $(DECOMPLIB)

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

其他公開可用的用於求解 Navier-Stokes 方程的程式

還有其他幾個公開可用的基於快速傅立葉變換的程式可用於求解 Navier-Stokes 方程。其中包括

HIT 3d，可從 http://code.google.com/p/hit3d/ 獲取

Tarang，可從 http://turbulencehub.org/index.php/codes/tarang/ 獲取

SpectralDNS，可從 https://github.com/spectralDNS 獲取

以及

Turbo，可從 http://aqua.ulb.ac.be/home/?page_id=50 獲取

練習

使用 2DECOMP&FFT 來編寫一個二維 Navier-Stokes 求解器。該庫被構建為進行三維 FFT，但是透過選擇其中一個數組只有一個條目，該庫就可以在分散式記憶體機器上進行二維 FFT。
Uecker^[18] 描述了二維各向同性湍流中渦量^[19] 功率譜的預期冪律縮放。查詢 Uecker^[20]，然後嘗試生成數值資料，以驗證儘可能多的波數空間數量級的功率縮放定律，這些數量級在您可訪問的計算資源上是可行的。Boffetta 和 Ecke^[21] 提供了該領域研究工作的最新概述。Fornberg^[22] 討論瞭如何計算功率譜。
如果我們設定 $\mu =0$ ，則 Navier-Stokes 方程變為尤拉方程。嘗試使用隱式中點規則或 Crank-Nicolson 方法來模擬二維或三維尤拉方程。看看你是否能找到好的迭代方案來做到這一點，你可能需要使用牛頓迭代。Majda 和 Bertozzi^[23] 介紹了尤拉方程。
Taylor-Green 渦流初始條件已被研究作為一種可能的流動，它可能在尤拉方程和 Navier-Stokes 方程中導致速度梯度絕對值最大值的爆破。在許多這樣的模擬中，對稱性已被用來獲得更高的有效解析度，例如，參見 Cichowlas 和 Brachet^[24]。考慮使用 Kida-Pelz 渦流或 Taylor-Green 渦流作為尤拉方程的初始條件，並新增非對稱擾動。如果您無法使隱式時間步進方案執行，請考慮使用顯式方案，例如龍格-庫塔方法。與文獻中的先前研究相比，該流動是如何演化的？Majda 和 Bertozzi^[25] 介紹了尤拉方程的爆破。
我們編寫的三維程式不是最高效的，因為可以使用實數到複數變換來將工作量減少一半。在其中一個 Navier-Stokes 程式中實現實數到複數變換。
我們編寫的程式也會引入一些混疊誤差。透過閱讀關於譜方法的書籍，例如 Canuto 等人^[26]^[27]，找出什麼是混疊誤差。解釋 Johnstone^[28] 中解釋的策略為何可以減少混疊誤差。
看看你是否可以找到其他針對 Naiver-Stokes 方程的開源快速傅立葉變換求解器，並將它們新增到上面的列表中。

備註

↑ Tritton (1988)
↑ Doering 和 Gibbon (1995)
↑ Gallavotti (2002)
↑ Uecker
↑ Canuto 等人 (2006)
↑ Canuto 等人 (2007)
↑ Tritton (1988)
↑ Temam (2001)
↑ Gallavotti (2002)
↑ Doering 和 Gibbon (1995)
↑ Orszag 和 Patterson (1972)
↑ Canuto 等人 (2006)
↑ Canuto 等人 (2007)
↑ Shapiro (1993)
↑ Majda 和 Bertozzi (2002)
↑ Laizet 和 Lamballais (2009)
↑ 我們沒有證明空間離散化的收斂性，因此該結果假設空間離散化尚未完成。
↑ Uecker (2009)
↑ 渦量是渦旋的平方。
↑ Uecker (2009)
↑ Boffetta 和 Ecke (2012)
↑ Fornberg (1977)
↑ Majda 和 Bertozzi (2002)
↑ Cichowlas 和 Brachet (2005)
↑ Majda 和 Bertozzi (2012)
↑ Canuto 等人 (2006)
↑ Canuto 等人 (2007)
↑ Johnstone (2012)

參考文獻

Boffetta, G.; Ecke, R.E. (2012). "二維湍流". 流體力學年鑑. 44: 427–451.

Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. (2006). 譜方法：單域基礎. Springer.

Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. (2007). 譜方法：向複雜幾何的演變及流體力學應用. Springer.

Cichowlas, C.; Brachet, M.-E. (2005). "複雜奇點在 Kida-Pelz 和 Taylor-Green 無粘流中的演化". 流體動力學研究. 36: 239–248.

Doering, C.R.; Gibbon, J.D. (1995). Navier-Stokes 方程的應用分析. 劍橋大學出版社.

Fornberg, B. (1977). "二維湍流的數值研究". 計算物理學報. 25: 1–31.

Gallavotti, G. (2002). 流體動力學基礎. 施普林格.

Gottlieb, D.; Orszag, S.A. (1977). 譜方法的數值分析：理論與應用. SIAM.

Johnstone, R. (2012). "湍流直接數值模擬的改進縮放". HECTOR 分散式計算科學與工程報告. {{cite web}}: 引用具有空的未知引數：|coauthors= (幫助)

Laizet, S.; Lamballais, E. (2009). "不可壓縮流的高階緊緻格式：一種簡單高效的準譜精度方法". 計算物理學報. 238: 5989–6015.

Bertozzi, A.L. (2002). 渦量與不可壓縮流. 劍橋大學出版社.

Orszag, S.A.; Patterson Jr., G.S. (1979). 物理評論快報. 28 (2): 76–79. {{cite journal}}: 缺少或空的|title= (幫助)

Peyret, R. (2002). 不可壓縮粘性流的譜方法. 施普林格.

Shapiro, A (1993). "Navier-Stokes 方程精確解在三維非靜力數值模型驗證測試中的應用". 每月天氣評論. 121: 2420–2425. {{cite journal}}: 引用具有空的未知引數：|coauthors= (幫助)

Temam, R. (2001). Navier-Stokes 方程 (第三版). 美國數學學會.

Tritton, D.J. (1988). 物理流體動力學. 牛津大學出版社.

Uecker, H. 關於拋物線型偏微分方程和 Navier-Stokes 方程的譜方法簡短的特別介紹. {{cite book}}: 引用具有空的未知引數：|coauthors= (幫助)

[1] Tritton (1988)

[2] Doering 和 Gibbon (1995)

[3] Gallavotti (2002)

[4] Uecker

[5] Canuto 等人 (2006)

[6] Canuto 等人 (2007)

[7] Tritton (1988)

[8] Temam (2001)

[9] Gallavotti (2002)

[10] Doering 和 Gibbon (1995)

[11] Orszag 和 Patterson (1972)

[12] Canuto 等人 (2006)

[13] Canuto 等人 (2007)

[14] Shapiro (1993)

[15] Majda 和 Bertozzi (2002)

[16] Laizet 和 Lamballais (2009)

[17] 我們沒有證明空間離散化的收斂性，因此該結果假設空間離散化尚未完成。

[18] Uecker (2009)

[19] 渦量是渦旋的平方。

[20] Uecker (2009)

[21] Boffetta 和 Ecke (2012)

[22] Fornberg (1977)

[23] Majda 和 Bertozzi (2002)

[24] Cichowlas 和 Brachet (2005)

[25] Majda 和 Bertozzi (2012)

[26] Canuto 等人 (2006)

[27] Canuto 等人 (2007)

[28] Johnstone (2012)

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]