二维和三维 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://wikibooks.cn/wiki/Parallel_Spectral_Numerical_Methods/The_Cubic_Nonlinear_Schrodinger_Equation#equation_G 中的程序列表类似。
使用 OpenMP 指令将列表 E 中的三维 Navier-Stokes 方程的示例 Fortran 代码并行化。尽可能提高效率。
编写另一个代码，使用线程化的 FFTW 来进行三维 Navier-Stokes 方程的快速傅里叶变换。该代码的结构应该与 https://wikibooks.cn/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]