并行谱数值方法/三次非线性薛定谔方程

三次非线性薛定谔方程

背景

三次非线性薛定谔方程出现在许多领域，包括量子力学、非线性光学和水面波。可以在http://en.wikipedia.org/wiki/Schrodinger_equation和http://en.wikipedia.org/wiki/Nonlinear_Schrodinger_equation找到一般介绍。薛定谔方程的数学介绍可以在 Sulem 和 Sulem和杨中找到。在本节中，我们将介绍算子分裂的概念，然后解释如何将其应用于一维、二维和三维的非线性薛定谔方程。在一维中，可以证明三次非线性薛定谔方程是亚临界的，因此存在所有时间的解。在二维中，它是$H^{1}$临界的，因此解可能表现出$H^{1}$范数的爆炸，即解的梯度的平方积分可以在有限时间内变为无穷大。最后，在三维中，非线性薛定谔方程是$L^{2}$超临界的，因此解的平方积分也可以在有限时间内变为无穷大。有关范数和希尔伯特空间的介绍，请参阅偏微分方程或分析的教科书，例如 Evans、Linares 和 Ponce、Lieb 和 Loss或 Renardy 和 Rogers。一个有趣的问题是如何发生这种爆炸，数值模拟经常被用来理解这一点；参见 Sulem 和 Sulem以获取此类示例。三次非线性薛定谔方程由下式给出

$i\psi _{t}+\Delta \psi \pm |\psi |^{2}\psi =0,$

其中 $\psi$ 是波函数，$\Delta$ 是拉普拉斯算子，因此在一维中它是 $\partial _{xx}$，在二维中，$\partial _{xx}+\partial _{yy}$，在三维中它是 $\partial _{xx}+\partial _{yy}+\partial _{zz}$。$+$ 对应于聚焦三次非线性薛定谔方程，而 $-$ 对应于散焦三次非线性薛定谔方程。该方程具有许多守恒量，包括“质量”，

$\int _{\Omega }|\psi |^{2}{\text{d}}^{n}\mathbf {x}$

和“能量”，

$\int _{\Omega }{\frac {1}{2}}|\nabla \psi |^{2}\mp {\frac {1}{4}}|\psi |^{4}{\text{d}}^{n}\mathbf {x}$

其中 $n$ 是维度，$\Omega$ 是解的域。正如 Klein 解释的那样，这两个量可以为数值生成解的精度提供有用的检验。

分裂

我们将考虑一种称为分裂的数值方法来求解此方程。这种方法在多个应用中出现，并且当方程可以分裂成两个单独的方程时，它是一种有用的数值方法，每个方程都可以精确求解，或者每个部分最适合由不同的数值方法求解。 Holden 等人，McLachlan 和 Quispel，Thalhammer，Shen、Tang 和 Wang，Weideman 和 Herbst 以及 Yang 提供了对分裂的介绍，http://en.wikipedia.org/wiki/Split-step_method 也对此进行了介绍。 对于那些对非线性薛定谔方程的时间步进方法比较感兴趣的人，请参阅 Klein。 为了描述该方法的基本原理，我们考虑 Holden 等人 给出的一个例子，即常微分方程：

$u_{t}=u(u-1)$ 其中 $u(t=0)=0.8$.

( 1)

我们可以通过变量分离相对简单地求解此方程，发现

$u(t)={\frac {4}{4+\exp(t)}}.$

现在，一个有趣的观察结果是，我们也可以单独求解方程 $u_{t}=u^{2}$ 和 $u_{t}=-u$。对于第一个，我们得到 $u(t)={\frac {u(0)}{1-tu(0)}}$，对于第二个，我们得到 $u(t)=u(0)\exp(-t)$。分裂背后的原理是交替地求解这两个单独的方程，时间很短。我们将描述 Strang 分裂，虽然还有其他形式的分裂，例如 Godunov 分裂以及加性分裂。我们在这里不介绍这些，而是将您推荐给前面提到的参考文献，特别是 Holden 等人。 为了理解如何使用分裂求解微分方程，请考虑线性常微分方程

$u_{t}=u+2u$ 其中 $u(0)=1$.

We can first solve $p_{t}=p$ for a time $\delta t/2$ and then using $q(0)=p(\delta t/2)$, we solve $q_{t}=2q$ also for a time $\delta t$ to get $q(\delta t)$ and finally solve $r_{t}=r$ for a time $\delta t/2$ with initial data $r(0)=q(\delta t)$. Thus in this case $p(\delta t)=\exp(\delta t/2)$, $q(\delta t)=p(\delta t/2)\exp(2\delta t)=\exp(5\delta t/2)$ and $u(\delta t)\approx r(\delta t/2)=q(\delta t)\exp(\delta t/2)=\exp(3\delta t)$, which in this case is the exact solution. One can perform a similar splitting for matrix differential equations. Consider solving $\mathbf {u} _{t}=(\mathbf {A} +\mathbf {B} )\mathbf {u}$, where $\mathbf {A}$ and $\mathbf {B}$ are $n\times n$ matrices, the exact solution is $\mathbf {u} =\exp \left((\mathbf {A} +\mathbf {B} )t\right)\mathbf {u} (t=0)$, and an approximate solution produced after one time step of splitting is $u(\delta t)\approx u(0)\exp(\mathbf {A} \delta t)\exp(\mathbf {B} \delta t)$, which is not in general equal to $u(t=0)\exp \left((\mathbf {A} +\mathbf {B} )\delta t\right)$ unless the matrices $\mathbf {A}$ and $\mathbf {B}$ commute, and so the error in doing splitting in this case is of the form $(\mathbf {A} \mathbf {B} -\mathbf {B} \mathbf {A} )\delta t$. Listing A uses Matlab to demonstrate how to do splitting for eq. 1 .

( A)

使用 Strang 分裂求解 ODE 的 Matlab 程序 代码下载

% A program to solve the u_t=u(u-1) using a
% Strang Splitting method

clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
    'defaultaxesfontweight','bold')
Nt = 1000; % number of time slices
tmax = 1; % maximum time
dt=tmax/Nt; % increment between times
time=(linspace(1,Nt,Nt)-1)*dt; % time
uexact=4./(4+exp(time)); % exact solution
u(1)=0.8

for i=1:Nt-1
    c=-1/u(i);
    utemp=-1/(c+0.5*dt);
    utemp2=utemp*exp(-dt);
    c=-1/utemp2;
    u(i+1)=-1/(c+0.5*dt);
end
figure(1)
plot(time,u,'r+',time,uexact,'b-');

( Ap)

使用 Strang 分裂求解 ODE 的 Python 程序。 代码下载

"""
A program to solve u_t'=u(u-1) using a Strang
splitting method
"""

import time
import numpy
import matplotlib.pyplot as plt

Nt = 1000	# Number of timesteps
tmax = 1.0	# Maximum time
dt=tmax/Nt # increment between times
u0 = 0.8 # Initial value
t0 = 0	# Starting time
u = [u0]	# Variables holding the values of iterations
t = [t0]	# Times of discrete solutions

T0 = time.clock()
for i in xrange(Nt):
c = -1.0/u[i]
utemp=-1.0/(c+0.5*dt)
utemp2=utemp*numpy.exp(-dt)
c = -1.0/utemp2
unew=-1.0/(c+0.5*dt)
u.append(unew)
t.append(t[i]+dt)

T = time.clock() - T0	
uexact = [4.0/(4.0+numpy.exp(tt)) for tt in t]

print
print "Elapsed time is %f" % (T)

plt.figure()
plt.plot(t,u,'r+',t,uexact,'b-.')
plt.xlabel('Time')
plt.ylabel('Solution')
plt.legend(('Numerical Solution', 'Exact solution'))
plt.title('Numerical solution of du/dt=u(u-1)')
plt.show()

练习

1. 修改 Matlab 代码，以计算时间 1 处的误差，并选择多个不同的时间步长。 数值验证 Strang 分裂的二阶精度。
2. 修改 Matlab 代码以使用 Godunov 分裂，其中求解 $u1_{t}=u1$ 时间为 $\delta t$，然后使用 $u1(\delta t)$ 作为初始数据，求解 $u2_{t}=2u2$ 也为时间 $\delta t$，以获得对 $u(\delta t)$ 的近似值。计算时间为 1 时的误差，对于不同的时间步长选择。数值验证 Godunov 分裂是一阶精度的。

串行

对于非线性薛定谔方程

$i\psi _{t}\pm |\psi |^{2}\psi +\Delta \psi =0$,

( 2)

我们首先求解

$i\psi _{t}+\Delta \psi =0$

( 3)

使用傅里叶变换精确求解，得到

$\psi (\delta {t}/2,\cdot )$.

然后我们求解

$i\psi _{t}\pm |\psi |^{2}\psi =0$

( 4)

其中

$\psi (\delta {t}/2,\cdot )$

作为初始数据，时间步长为 $\delta t$。正如 Klein 和 Thalhammer 所解释的，这可以在实空间中精确求解，因为在等式 4 中，$|\psi |^{2}$ 是空间和时间的每个点的守恒量。为了证明这一点，令 $\psi ^{*}$ 表示 $\psi$ 的复共轭，所以

${\frac {\mathrm {d} |\psi |^{2}}{\mathrm {d} t}}=\psi ^{*}{\frac {\mathrm {d} \psi }{\mathrm {d} t}}+{\frac {\mathrm {d} \psi ^{*}}{\mathrm {d} t}}\psi =\psi ^{*}\left(\pm i|\psi |^{2}\psi \right)+\left(\pm i|\psi |^{2}\psi \right)^{*}\psi =0$.

然后使用方程式 3 计算另一个半步，该半步使用通过求解方程式 4 生成的解来获得时间 $\delta t$ 处的近似解。下面是演示拆分的示例 Matlab 代码。

非线性薛定谔方程的示例 Matlab 和 Python 程序

清单 B 中的程序计算对聚焦非线性薛定谔方程的显式已知精确解的近似值。

( B)

一个使用 Strang 拆分求解一维非线性薛定谔方程的 Matlab 程序 代码下载

% A program to solve the nonlinear Schr\"{o}dinger equation using a
% splitting method

clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
    'defaultaxesfontweight','bold')

Lx = 20; % period 2*pi * L
Nx = 16384; % number of harmonics
Nt = 1000; % number of time slices
dt = 0.25*pi/Nt; % time step
U=zeros(Nx,Nt/10);

Es = -1; % focusing or defocusing parameter

% 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
k2x = kx.^2; % square of wave vector
% initial conditions
t=0; tdata(1)=t;
u=4*exp(1i*t)*(cosh(3*x)+3*exp(8*1i*t)*cosh(x))...
    ./(cosh(4*x)+4*cosh(2*x)+3*cos(8*t));
v=fft(u);
figure(1); clf; plot(x,u);xlim([-2,2]); drawnow;
U(:,1)=u;

% mass
ma = fft(abs(u).^2);
ma0 = ma(1);

% solve pde and plot results
for n =2:Nt+1
    
    vna=exp(0.5*1i*dt*k2x).*v;
    una=ifft(vna);
    pot=2*(una.*conj(una));
    unb=exp(-1i*Es*dt*pot).*una;
    vnb=fft(unb);
    v=exp(0.5*1i*dt*k2x).*vnb;
    t=(n-1)*dt;
    
    if (mod(n,10)==0)
        tdata(n/10)=t;
        u=ifft(v);
        U(:,n/10)=u;
        uexact=4*exp(1i*t)*(cosh(3*x)+3*exp(8*1i*t)*cosh(x))...
            ./(cosh(4*x)+4*cosh(2*x)+3*cos(8*t));
        figure(1); clf; plot(x,abs(u).^2); ...
            xlim([-0.5,0.5]); title(num2str(t));
        figure(2); clf; plot(x,abs(u-uexact).^2);...
            xlim([-0.5,0.5]); title(num2str(t));
        drawnow;
        ma = fft(abs(u).^2);
        ma = ma(1);
        test = log10(abs(1-ma/ma0))
    end
end
figure(3); clf; mesh(tdata(1:(n-1)/10),x,abs(U(:,1:(n-1)/10)).^2);

( Bp)

一个使用 Strang 拆分求解一维非线性薛定谔方程的 Python 程序。 代码下载

"""
A program to solve the 1D Nonlinear Schrodinger equation using a
second order splitting method. The numerical solution is compared
to an exact soliton solution. The exact equation solved is
iu_t+u_{xx}+2|u|^2u=0

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
import matplotlib.pyplot as plt
import time

plt.ion()

# Grid
Lx=16.0 # Period 2*pi*Lx
Nx=8192 # Number of harmonics
Nt=1000 # Number of time slices
tmax=1.0 # Maximum time
dt=tmax/Nt # time step
plotgap=10 # time steps between plots
Es= -1.0 # focusing (+1) or defocusing (-1) parameter
numplots=Nt/plotgap # number of plots to make

x = [i*2.0*math.pi*(Lx/Nx) for i in xrange(-Nx/2,1+Nx/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)])

k2xm=numpy.zeros((Nx), dtype=float)
xx=numpy.zeros((Nx), dtype=float)

for i in xrange(Nx):
   k2xm[i] = numpy.real(k_x[i]**2)
   xx[i]=x[i]
        

# allocate arrays
usquared=numpy.zeros((Nx), dtype=float)
pot=numpy.zeros((Nx), dtype=float)
u=numpy.zeros((Nx), dtype=complex)
uexact=numpy.zeros((Nx), dtype=complex)
una=numpy.zeros((Nx), dtype=complex)
unb=numpy.zeros((Nx), dtype=complex)
v=numpy.zeros((Nx), dtype=complex)
vna=numpy.zeros((Nx), dtype=complex)
vnb=numpy.zeros((Nx), dtype=complex)
mass=numpy.zeros((Nx), dtype=complex)
test=numpy.zeros((numplots-1),dtype=float)
tdata=numpy.zeros((numplots-1), dtype=float)

t=0.0
u=4.0*numpy.exp(complex(0,1.0)*t)*\
   (numpy.cosh(3.0*xx)+3.0*numpy.exp(8.0*complex(0,1.0)*t)*numpy.cosh(xx))\
   /(numpy.cosh(4*xx)+4.0*numpy.cosh(2.0*xx)+3.0*numpy.cos(8.0*t))
uexact=u
v=numpy.fft.fftn(u)
usquared=abs(u)**2
fig =plt.figure()
ax = fig.add_subplot(311)
ax.plot(xx,numpy.real(u),'b-')
plt.xlabel('x')
plt.ylabel('real u')
ax = fig.add_subplot(312)
ax.plot(xx,numpy.imag(u),'b-')
plt.xlabel('x')
plt.ylabel('imaginary u')
ax = fig.add_subplot(313)
ax.plot(xx,abs(u-uexact),'b-')
plt.xlabel('x')
plt.ylabel('error')
plt.show()


# initial mass
usquared=abs(u)**2
mass=numpy.fft.fftn(usquared)
ma=numpy.real(mass[0])
maO=ma
tdata[0]=t
plotnum=0
#solve pde and plot results
for nt in xrange(numplots-1):
    for n in xrange(plotgap):
        vna=v*numpy.exp(complex(0,0.5)*dt*k2xm)
        una=numpy.fft.ifftn(vna)
        usquared=2.0*abs(una)**2
        pot=Es*usquared
        unb=una*numpy.exp(complex(0,-1)*dt*pot)
        vnb=numpy.fft.fftn(unb)
        v=vnb*numpy.exp(complex(0,0.5)*dt*k2xm)
        u=numpy.fft.ifftn(v)
        t+=dt
    plotnum+=1
    usquared=abs(u)**2
    uexact = 4.0*numpy.exp(complex(0,1.0)*t)*\
      (numpy.cosh(3.0*xx)+3.0*numpy.exp(8.0*complex(0,1.0)*t)*numpy.cosh(xx))\
      /(numpy.cosh(4*xx)+4.0*numpy.cosh(2.0*xx)+3.0*numpy.cos(8.0*t))
    ax = fig.add_subplot(311)
    plt.cla()
    ax.plot(xx,numpy.real(u),'b-')
    plt.title(t)
    plt.xlabel('x')
    plt.ylabel('real u')
    ax = fig.add_subplot(312)
    plt.cla()
    ax.plot(xx,numpy.imag(u),'b-')
    plt.xlabel('x')
    plt.ylabel('imaginary u')
    ax = fig.add_subplot(313)
    plt.cla()
    ax.plot(xx,abs(u-uexact),'b-')
    plt.xlabel('x')
    plt.ylabel('error')
    plt.draw()
    mass=numpy.fft.fftn(usquared)
    ma=numpy.real(mass[0])
    test[plotnum-1]=numpy.log(abs(1-ma/maO))
    print(test[plotnum-1])
    tdata[plotnum-1]=t
       
plt.ioff()
plt.show()

( C)

一个使用 Strang 拆分求解二维非线性薛定谔方程的 Matlab 程序 代码下载

% A program to solve the 2D nonlinear Schr\"{o}dinger equation using a
% splitting method

clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,'defaultaxesfontweight','bold')

% set up grid
tic
Lx = 20; % period 2*pi*L
Ly = 20; % period 2*pi*L
Nx = 2*256; % number of harmonics
Ny = 2*256; % number of harmonics
Nt = 100; % number of time slices
dt = 5.0/Nt; % time step

Es = 1.0;

% 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
[xx,yy]=meshgrid(x,y);
[k2xm,k2ym]=meshgrid(kx.^2,ky.^2);
% initial conditions
u = exp(-(xx.^2+yy.^2));
v=fft2(u);
figure(1); clf; mesh(xx,yy,u); drawnow;
t=0; tdata(1)=t;

% mass
ma = fft2(abs(u).^2);
ma0 = ma(1,1);

% solve pde and plot results
for n =2:Nt+1
    vna=exp(0.5*1i*dt*(k2xm + k2ym)).*v;
    una=ifft2(vna);
    pot=Es*((abs(una)).^2);
    unb=exp(-1i*dt*pot).*una;
    vnb=fft2(unb);
    v=exp(0.5*1i*dt*(k2xm + k2ym)).*vnb;
    u=ifft2(v);
    t=(n-1)*dt;
    tdata(n)=t;
     if (mod(n,10)==0)
         figure(2); clf; mesh(xx,yy,abs(u).^2); title(num2str(t));
         drawnow;
         ma = fft2(abs(u).^2);
         ma = ma(1,1);
         test = log10(abs(1-ma/ma0))
     end
end
figure(4); clf; mesh(xx,yy,abs(u).^2);
toc

( Cp)

一个使用 Strang 拆分求解二维非线性薛定谔方程的 Python 程序。 代码下载

"""
A program to solve the 2D Nonlinear Schrodinger equation using a
second order splitting method

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=4.0 # Period 2*pi*Lx
Ly=4.0 # Period 2*pi*Ly
Nx=64 # Number of harmonics
Ny=64 # Number of harmonics
Nt=100 # Number of time slices
tmax=1.0 # Maximum time
dt=tmax/Nt # time step
plotgap=10 # time steps between plots
Es= 1.0 # focusing (+1) or defocusing (-1) parameter
numplots=Nt/plotgap # number of plots to make

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

k2xm=numpy.zeros((Nx,Ny), dtype=float)
k2ym=numpy.zeros((Nx,Ny), dtype=float)
xx=numpy.zeros((Nx,Ny), dtype=float)
yy=numpy.zeros((Nx,Ny), dtype=float)


for i in xrange(Nx):
    for j in xrange(Ny):
            k2xm[i,j] = numpy.real(k_x[i]**2)
            k2ym[i,j] = numpy.real(k_y[j]**2)
            xx[i,j]=x[i]
            yy[i,j]=y[j]
        

# allocate arrays
usquared=numpy.zeros((Nx,Ny), dtype=float)
pot=numpy.zeros((Nx,Ny), dtype=float)
u=numpy.zeros((Nx,Ny), dtype=complex)
una=numpy.zeros((Nx,Ny), dtype=complex)
unb=numpy.zeros((Nx,Ny), dtype=complex)
v=numpy.zeros((Nx,Ny), dtype=complex)
vna=numpy.zeros((Nx,Ny), dtype=complex)
vnb=numpy.zeros((Nx,Ny), dtype=complex)
mass=numpy.zeros((Nx,Ny), dtype=complex)
test=numpy.zeros((numplots-1),dtype=float)
tdata=numpy.zeros((numplots-1), dtype=float)

u=numpy.exp(-(xx**2 + yy**2 ))
v=numpy.fft.fftn(u)
usquared=abs(u)**2
src = mlab.surf(xx,yy,usquared,colormap='YlGnBu',warp_scale='auto')
mlab.scalarbar()
mlab.xlabel('x',object=src)
mlab.ylabel('y',object=src)
mlab.zlabel('abs(u)^2',object=src)

# initial mass
usquared=abs(u)**2
mass=numpy.fft.fftn(usquared)
ma=numpy.real(mass[0,0])
print(ma)
maO=ma
t=0.0
tdata[0]=t
plotnum=0
#solve pde and plot results
for nt in xrange(numplots-1):
    for n in xrange(plotgap):
        vna=v*numpy.exp(complex(0,0.5)*dt*(k2xm+k2ym))
        una=numpy.fft.ifftn(vna)
        usquared=abs(una)**2
        pot=Es*usquared
        unb=una*numpy.exp(complex(0,-1)*dt*pot)
        vnb=numpy.fft.fftn(unb)
        v=vnb*numpy.exp(complex(0,0.5)*dt*(k2xm+k2ym) )
        u=numpy.fft.ifftn(v)
        t+=dt
    plotnum+=1
    usquared=abs(u)**2
    src.mlab_source.scalars = usquared
    mass=numpy.fft.fftn(usquared)
    ma=numpy.real(mass[0,0])
    test[plotnum-1]=numpy.log(abs(1-ma/maO))
    print(test[plotnum-1])
    tdata[plotnum-1]=t
       
plt.figure()
plt.plot(tdata,test,'r-')
plt.title('Time Dependence of Change in Mass')
plt.show()

( D)

一个使用 Strang 拆分求解三维非线性薛定谔方程的 Matlab 程序 代码下载

% A program to solve the 3D nonlinear Schr\"{o}dinger equation using a
% splitting method
% updated by Abdullah Bargash , AbdulAziz Al-thiban
% vol3d can be downloaded from http://www.mathworks.com/matlabcentral/fileexchange/22940-vol3d-v2
%
clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
    'defaultaxesfontweight','bold')

% set up grid
tic
Lx = 4; % period 2*pi*L
Ly = 4; % period 2*pi*L
Lz = 4; % period 2*pi*L
Nx = 64; % number of harmonics
Ny = 64; % number of harmonics
Nz = 64; % number of harmonics
Nt = 100; % number of time slices
dt = 1.0/Nt; % time step

Es = 1.0; % focusing or defocusing parameter

% 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);
[k2xm,k2ym,k2zm]=meshgrid(kx.^2,ky.^2,kz.^2);

% initial conditions
u = exp(-(xx.^2+yy.^2+zz.^2));
v=fftn(u);
figure(1); clf; UP = abs(u).^2;
    H = vol3d('CData',UP,'texture','3D','XData',x,'YData',y,'ZData',z);
    xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
    axis equal; axis square; view(3);
 
    xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); drawnow;
t=0; tdata(1)=t;
   
% mass
ma = fftn(abs(u).^2);
ma0 = ma(1,1,1);

% solve pde and plot results

for n =2:Nt+1
    vna=exp(0.5*1i*dt*(k2xm + k2ym + k2zm)).*v;
    una=ifftn(vna);
    pot=Es*((abs(una)).^2);
    unb=exp(-1i*dt*pot).*una;
    vnb=fftn(unb);
    v=exp(0.5*1i*dt*(k2xm + k2ym + k2zm)).*vnb;
    u=ifftn(v);
    t=(n-1)*dt;
    tdata(n)=t;
    if (mod(n,10)==0)
        
    figure(1); clf; UP = abs(u).^2;
   
    H = vol3d('CData',UP,'texture','3D','XData',x,'YData',y,'ZData',z);
    xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
    axis equal; axis square; view(3);
    title(num2str(t));

        
        xlabel('x'); ylabel('y'); zlabel('z');
        axis equal; axis square; view(3); drawnow;
        ma = fftn(abs(u).^2);
        ma = ma(1,1,1); test = log10(abs(1-ma/ma0))
         
    end
     
end

    figure(1); clf; UP = abs(u).^2;
   
    H = vol3d('CData',UP,'texture','3D','XData',x,'YData',y,'ZData',z);
    xlabel('x'); ylabel('y'); zlabel('z'); colorbar;
    axis equal; axis square; view(3);

非线性薛定谔方程的示例一维 Fortran 程序

在考虑并行程序之前，我们需要了解如何为一维非线性薛定谔方程编写 Fortran 代码。下面是一个示例 Fortran 程序，后面是一个 Matlab 绘图脚本，用于可视化结果。在编译 Fortran 程序时，需要标准 Fortran 编译器和 FFTW 库。由于所需命令与热方程 makefile 中的命令相似，因此此处不再赘述。

( E)

一个使用拆分求解一维非线性薛定谔方程的 Fortran 程序 代码下载

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program solves nonlinear Schrodinger equation in 1 dimension
! i*u_t+Es*|u|^2u+u_{xx}=0
! using a second order time spectral splitting scheme
!
! The boundary conditions are u(0)=u(2*L*\pi)
! The initial condition is u=exp(-x^2)
!
! .. Parameters ..
! Nx = number of modes in x - power of 2 for FFT
! Nt = number of timesteps to take
! Tmax = maximum simulation time
! plotgap = number of timesteps between plots
! 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
! L = width of box
! ES = +1 for focusing and -1 for defocusing
! .. Scalars ..
! i = loop counter in x direction
! n = loop counter for timesteps direction
! allocatestatus = error indicator during allocation
! 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
! dt = timestep
! .. Arrays ..
! u = approximate solution
! v = Fourier transform of approximate solution
! .. Vectors ..
! una = temporary field
! unb = temporary field
! vna = temporary field
! pot = potential
! kx = fourier frequencies in x direction
! x = x 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 x Fourier transform
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
! Check that the initial iterate is consistent with the
! boundary conditions for the domain specified
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! FFTW3 -- Fast Fourier Transform in the West Library
! (http://www.fftw.org/)


PROGRAM main

! Declare variables
IMPLICIT NONE
INTEGER(kind=4), PARAMETER :: Nx=8*256	
INTEGER(kind=4), PARAMETER	:: Nt=200	
REAL(kind=8), PARAMETER	&
:: pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8), PARAMETER	:: L=5.0d0	
REAL(kind=8), PARAMETER	:: Es=1.0d0	
REAL(kind=8)	:: dt=2.0d0/Nt	
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx	
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x	
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE	:: u	
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE	:: v
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: una,vn
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: unb,pot
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: time
INTEGER(kind=4)	:: i,j,k,n,modes,AllocateStatus
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)	:: planfx,planbx
   CHARACTER*100	:: name_config

CALL system_clock(start,count_rate)
ALLOCATE(kx(1:Nx),x(1:Nx),u(1:Nx,1:Nt+1),v(1:Nx,1:Nt+1),&
una(1:Nx),vn(1:Nx),unb(1:Nx),pot(1:Nx),time(1:Nt+1),&
fftfx(1:Nx),fftbx(1:Nx),stat=AllocateStatus)
IF (allocatestatus .ne. 0) STOP	
! set up ffts
CALL dfftw_plan_dft_1d_(planfx,Nx,fftfx(1:Nx),fftbx(1:Nx),&
FFTW_FORWARD,FFTW_PATIENT)
CALL dfftw_plan_dft_1d_(planbx,Nx,fftbx(1:Nx),fftfx(1:Nx),&
FFTW_BACKWARD,FFTW_PATIENT)
PRINT *,'Setup FFTs'
! setup fourier frequencies
DO i=1,1+Nx/2
kx(i)= cmplx(0.0d0,1.0d0)*(i-1.0d0)/L
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
x(i)=(-1.0d0 + 2.0d0*REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)))*pi*L
END DO
PRINT *,'Setup grid and fourier frequencies'

DO i=1,Nx
u(i,1)=exp(-1.0d0*(x(i)**2))
END DO
! transform initial data
CALL dfftw_execute_dft_(planfx,u(1:Nx,1),v(1:Nx,1))
PRINT *,'Got initial data, starting timestepping'
time(1)=0.0d0
DO n=1,Nt
time(n+1)=n*dt
DO i=1,Nx
vn(i)=exp(0.5d0*dt*kx(i)*kx(i)*cmplx(0.0d0,1.0d0))*v(i,n)
END DO
CALL dfftw_execute_dft_(planbx,vn(1:Nx),una(1:Nx))
! normalize
DO i=1,Nx
una(i)=una(1:Nx)/REAL(Nx,kind(0d0))	
pot(i)=Es*una(i)*conjg(una(i))
unb(i)=exp(cmplx(0.0d0,-1.0d0)*dt*pot(i))*una(i)
END DO
CALL dfftw_execute_dft_(planfx,unb(1:Nx),vn(1:Nx))
DO i=1,Nx
v(i,n+1)=exp(0.50d0*dt*kx(i)*kx(i)*cmplx(0.0d0,1.0d0))*vn(i)
END DO
CALL dfftw_execute_dft_(planbx,v(1:Nx,n+1),u(1:Nx,n+1))
! normalize
DO i=1,Nx
u(i,n+1)=u(i,n+1)/REAL(Nx,kind(0d0))
END DO
END DO
PRINT *,'Finished time stepping'
CALL system_clock(finish,count_rate)
PRINT*,'Program took ',&
REAL(finish-start,kind(0d0))/REAL(count_rate,kind(0d0)),'for execution'

name_config = 'u.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Nt
DO i=1,Nx
WRITE(11,*) abs(u(i,j))**2
END DO
END DO
CLOSE(11)

name_config = 'tdata.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Nt
WRITE(11,*) time(j)
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)

PRINT *,'Saved data'

CALL dfftw_destroy_plan_(planbx)
CALL dfftw_destroy_plan_(planfx)
CALL dfftw_cleanup_()

DEALLOCATE(kx,x,u,v,una,vn,unb,&
pot,time,fftfx,fftbx,&
stat=AllocateStatus)
IF (allocatestatus .ne. 0) STOP
PRINT *,'deallocated memory'
PRINT *,'Program execution complete'
END PROGRAM main

( F)

一个 Matlab 程序，它绘制清单 E 生成的对一维非线性薛定谔方程的数值解 代码下载

% A program to plot the computed results

clear all; format compact, format short,
set(0,'defaultaxesfontsize',18,'defaultaxeslinewidth',.9,...
    'defaultlinelinewidth',3.5,'defaultpatchlinewidth',5.5);

% Load data
load('./u.dat');
load('./tdata.dat');
load('./xcoord.dat');
Tsteps = length(tdata);

Nx = length(xcoord); Nt = length(tdata);
 
u = reshape(u,Nx,Nt);

% Plot data
figure(3); clf; mesh(tdata,xcoord,u); xlabel t; ylabel x; zlabel('|u|^2');

共享内存并行：OpenMP

我们回顾一下，OpenMP 是一组编译器指令，可以让您轻松地让 Fortran、C 或 C++ 程序在共享内存机器上运行——也就是说，所有计算进程都可以访问同一全局寻址内存空间的计算机。它允许对已用上述语言之一编写的串行程序进行轻松的并行化。

我们将演示二维非线性薛定谔方程的一种并行形式，其中我们将使用 OpenMP 命令对循环进行并行化，但将使用线程化的 FFTW 库为我们并行化变换。示例程序在清单 E 中。在练习中概述了使用 OpenMP 命令显式并行化循环和快速傅里叶变换的第二种方法。

( G)

一个使用拆分和线程化 FFTW 求解二维非线性薛定谔方程的 OpenMP Fortran 程序 代码下载

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program solves nonlinear Schrodinger equation in 2 dimensions
! i*u_t+Es*|u|^2u+u_{xx}+u_{yy}=0
! using a second order time spectral splitting scheme
!
! The boundary conditions are u(x=0,y)=u(2*Lx*\pi,y),
! u(x,y=0)=u(x,y=2*Ly*\pi)
! The initial condition is u=exp(-x^2-y^2)
!
! .. 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
! plotgap = number of timesteps between plots
! 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
! Lx = width of box in x direction
! Ly = width of box in y direction
! ES = +1 for focusing and -1 for defocusing
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
! n = loop counter for timesteps direction
! allocatestatus = error indicator during allocation
! numthreads = number of openmp threads
! ierr = error return code
! 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 = approximate solution
! v = Fourier transform of approximate solution
! unax = temporary field
! vnax = temporary field
! vnbx = temporary field
! vnay = temporary field
! vnby = temporary field
! potx = potential
! .. Vectors ..
! kx = fourier frequencies in x direction
! ky = 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 x Fourier transform
! fftfy = array to setup y Fourier transform
! fftby = array to setup y Fourier transform
!
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
! Check that the initial iterate is consistent with the
! boundary conditions for the domain specified
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! FFTW3 -- Fast Fourier Transform in the West Library
! (http://www.fftw.org/)
! OpenMP library
PROGRAM main
USE omp_lib	
IMPLICIT NONE	
! Declare variables
INTEGER(kind=4), PARAMETER :: Nx=1024
INTEGER(kind=4), PARAMETER :: Ny=1024	
INTEGER(kind=4), PARAMETER	:: Nt=20	
INTEGER(kind=4), PARAMETER	:: plotgap=5	
REAL(kind=8), PARAMETER	:: &
pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8), PARAMETER	:: Lx=2.0d0	
REAL(kind=8), PARAMETER	:: Ly=2.0d0	
REAL(kind=8), PARAMETER	:: Es=1.0d0	
REAL(kind=8)	:: dt=0.10d0/Nt	
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx	
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: ky	
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x	
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: y	
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: unax,vnax,vnbx,potx
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: vnay,vnby
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: time
INTEGER(kind=4)	:: i,j,k,n,allocatestatus,ierr
INTEGER(kind=4)	:: start, finish, count_rate, numthreads
INTEGER(kind=8), PARAMETER :: FFTW_IN_PLACE=8, FFTW_MEASURE=0,&
FFTW_EXHAUSTIVE=8, FFTW_PATIENT=32,&
                    FFTW_ESTIMATE=64
INTEGER(kind=8),PARAMETER :: FFTW_FORWARD=-1, FFTW_BACKWARD=1	
INTEGER(kind=8)	:: planfxy,planbxy
    CHARACTER*100	:: name_config,number_file

numthreads=omp_get_max_threads()
PRINT *,'There are ',numthreads,' threads.'

ALLOCATE(kx(1:Nx),ky(1:Nx),x(1:Nx),y(1:Nx),unax(1:Nx,1:Ny),&
vnax(1:Nx,1:Ny),potx(1:Nx,1:Ny),time(1:1+Nt/plotgap),&
stat=allocatestatus)	
IF (allocatestatus .ne. 0) stop
PRINT *,'allocated memory'

! set up multithreaded ffts
CALL dfftw_init_threads_(ierr)
PRINT *,'Initiated threaded FFTW'
CALL dfftw_plan_with_nthreads_(numthreads)
PRINT *,'Inidicated number of threads to be used in planning'
CALL dfftw_plan_dft_2d_(planfxy,Nx,Ny,unax(1:Nx,1:Ny),vnax(1:Nx,1:Ny),&
FFTW_FORWARD,FFTW_ESTIMATE)
  CALL dfftw_plan_dft_2d_(planbxy,Nx,Ny,vnax(1:Nx,1:Ny),unax(1:Nx,1:Ny),&
  FFTW_BACKWARD,FFTW_ESTIMATE)
PRINT *,'Setup FFTs'

! setup fourier frequencies
!$OMP PARALLEL PRIVATE(i,j)
!$OMP DO SCHEDULE(static)
DO i=1,1+Nx/2
kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/Lx
END DO
!$OMP END DO
kx(1+Nx/2)=0.0d0
!$OMP DO SCHEDULE(static)
DO i = 1,Nx/2 -1
kx(i+1+Nx/2)=-kx(1-i+Nx/2)
END DO
!$OMP END DO
!$OMP DO SCHEDULE(static)
   DO i=1,Nx
x(i)=(-1.0d0+2.0d0*REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)) )*pi*Lx
END DO
!$OMP END DO
!$OMP DO SCHEDULE(static)
DO j=1,1+Ny/2
ky(j)= cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0))/Ly
END DO
!$OMP END DO
ky(1+Ny/2)=0.0d0
!$OMP DO SCHEDULE(static)
DO j = 1,Ny/2 -1
ky(j+1+Ny/2)=-ky(1-j+Ny/2)
END DO
!$OMP END DO
!$OMP DO SCHEDULE(static)
   DO j=1,Ny
y(j)=(-1.0d0+2.0d0*REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0)) )*pi*Ly
END DO
!$OMP END DO
PRINT *,'Setup grid and fourier frequencies'
!$OMP DO SCHEDULE(static)
DO j=1,Ny
unax(1:Nx,j)=exp(-1.0d0*(x(1:Nx)**2 +y(j)**2))
END DO
!$OMP END DO
!$OMP END PARALLEL
name_config = 'uinitial.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Ny
DO i=1,Nx
WRITE(11,*) abs(unax(i,j))**2
END DO
END DO
CLOSE(11)
! transform initial data and do first half time step
CALL dfftw_execute_dft_(planfxy,unax(1:Nx,1:Ny),vnax(1:Nx,1:Ny))

PRINT *,'Got initial data, starting timestepping'
time(1)=0.0d0
CALL system_clock(start,count_rate)	
DO n=1,Nt	
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
DO i=1,Nx
vnax(i,j)=exp(0.5d0*dt*(kx(i)*kx(i) + ky(j)*ky(j))&
*cmplx(0.0d0,1.0d0))*vnax(i,j)
END DO
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planbxy,vnax(1:Nx,1:Ny),unax(1:Nx,1:Ny))
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
DO i=1,Nx
unax(i,j)=unax(i,j)/REAL(Nx*Ny,kind(0d0))
potx(i,j)=Es*unax(i,j)*conjg(unax(i,j))
unax(i,j)=exp(cmplx(0.0d0,-1.0d0)*dt*potx(i,j))&
*unax(i,j)
END DO
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planfxy,unax(1:Nx,1:Ny),vnax(1:Nx,1:Ny))
!$OMP PARALLEL DO PRIVATE(i) SCHEDULE(static)
DO j=1,Ny
DO i=1,Nx
vnax(i,j)=exp(0.5d0*dt*(kx(i)*kx(i) + ky(j)*ky(j))&	
*cmplx(0.0d0,1.0d0))*vnax(i,j)
END DO
END DO
!$OMP END PARALLEL DO
IF (mod(n,plotgap)==0) then
time(1+n/plotgap)=n*dt
PRINT *,'time',n*dt
CALL dfftw_execute_dft_(planbxy,vnax(1:Nx,1:Ny),unax(1:Nx,1:Ny))
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
DO i=1,Nx
unax(i,j)=unax(i,j)/REAL(Nx*Ny,kind(0d0))
END DO
END DO
!$OMP END PARALLEL DO
name_config='./data/u'
WRITE(number_file,'(i0)') 10000000+1+n/plotgap
ind=index(name_config,' ') -1
name_config=name_config(1:ind)//numberfile
ind=index(name_config,' ') -1
name_config=name_config(1:ind)//'.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Ny
DO i=1,Nx
WRITE(11,*) abs(unax(i,j))**2
END DO
END DO
CLOSE(11)
END IF
END DO
PRINT *,'Finished time stepping'
CALL system_clock(finish,count_rate)
PRINT*,'Program took ',REAL(finish-start)/REAL(count_rate),&
'for Time stepping'


name_config = 'tdata.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt/plotgap
WRITE(11,*) time(j)
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)
PRINT *,'Saved data'

CALL dfftw_destroy_plan_(planbxy)
CALL dfftw_destroy_plan_(planfxy)
CALL dfftw_cleanup_threads_()

DEALLOCATE(unax,vnax,potx,stat=allocatestatus)	
IF (allocatestatus .ne. 0) STOP
PRINT *,'Deallocated memory'

PRINT *,'Program execution complete'
END PROGRAM main

( H)

一个编译清单 G 中 OpenMP 程序的示例 makefile 代码下载

#define the complier
COMPILER = gfortran
# compilation settings, optimization, precision, parallelization
FLAGS = -O3 -fopenmp


# libraries
LIBS = -L/usr/local/lib -lfftw3 -lm
# source list for main program
SOURCES = NLSsplitting.f90

test: $(SOURCES)
${COMPILER} -o NLSsplitting $(FLAGS) $(SOURCES) $(LIBS)

clean:
rm *.o

clobber:
rm NLSsplitting

该示例假设您使用的是 Flux 并且已加载 GCC 编译器以及 GCC 编译的 FFTW 版本的环境。要使用 Intel 编译器使用此代码，需要明确设置 OMP 栈大小以确保足够大。如果您使用的是 PGI 编译器而不是 GCC 编译器，请将 makefile 中的标志 $-fopenmp$ 更改为 $-mp$。

( I)

一个 Matlab 程序，它绘制清单 G 生成的对二维非线性薛定谔方程的数值解 代码下载

% A program to plot the computed results for the 2D NLS equation

clear all; format compact, format short,
set(0,'defaultaxesfontsize',18,'defaultaxeslinewidth',.9,...
    'defaultlinelinewidth',3.5,'defaultpatchlinewidth',5.5);

% Load data
load('./ufinal.dat');
load('./tdata.dat');
load('./ycoord.dat');
load('./xcoord.dat');

Ny = length(ycoord); Nx = length(xcoord); Nt = length(tdata);
 
ufinal = reshape(ufinal,Nx,Ny);

% Plot data
figure(3); clf; mesh(xcoord,ycoord,ufinal); xlabel x; ylabel y; zlabel('|u|^2');

( J)

一个在 Flux 上使用的示例提交脚本。将 {your_username} 替换为您自己的用户名。 代码下载

#!/bin/bash
#PBS -N NLS
#PBS -l nodes=1:ppn=2,walltime=00:03:00
#PBS -q flux
#PBS -l qos=math471f11_flux
#PBS -A math471f11_flux
#PBS -M [email protected]
#PBS -m abe
#PBS -V
#
# Create a local directory to run and copy your files to local.
# Let PBS handle your output
cp ${HOME}/parallelspectralintro/NLSsplitting /nobackup/your_username/NLSsplitting
cd /nobackup/your_username

export OMP_NUM_THREADS=2
./NLSsplitting

练习

1. 下载与 Klein、Muite 和 Roidot 预印本 附带的示例 Matlab 程序。检查这些薛定谔型方程的质量和能量是如何计算的。在非线性薛定谔方程的 Matlab 程序中添加代码以检查质量和能量的守恒性。

2. Gross-Pitaevskii 方程 由 $i\psi _{t}+|\psi |^{2}\psi +V(\mathbf {x} )\psi =0$ 给出，其中我们将取 $V(\mathbf {x} )=\lVert \mathbf {x} \rVert _{l^{2}}^{2}=\sum _{k=1}^{N}x_{k}^{2}$，其中 $N$ 是空间维度。证明该方程可以通过将其拆分为以下两个方程来求解

   $i\psi _{t}+\Delta \psi =0$

   ( 5)

   和

   $i\psi _{t}+|\psi |^{2}\psi +V(\mathbf {x} )\psi =0$.

   ( 6)

   请务必解释如何求解公式 5 和 6 。

3. 修改Matlab代码以求解一维、二维和三维的Gross-Pitaevskii方程。

4. 修改串行Fortran代码以求解一维、二维和三维的Gross-Pitaevskii方程。

5. 清单 J 和 K 给出了一种并行化OpenMP程序的替代方法。使清单 G 中的程序尽可能高效，并尽可能类似于 J 中的程序，但不要更改并行化策略。比较两个不同程序的速度。尝试改变使用的网格点数和内核数。哪段代码在您的系统上更快？您认为这是为什么？

   
   

   ( J)

   一个使用分裂方法求解二维非线性薛定谔方程的OpenMP Fortran程序 代码下载

   !--------------------------------------------------------------------
   !
   !
   ! PURPOSE
   !
   ! This program solves nonlinear Schrodinger equation in 2 dimensions
   ! i*u_t+Es*|u|^2u+u_{xx}+u_{yy}=0
   ! using a second order time spectral splitting scheme
   !
   ! The boundary conditions are u(x=0,y)=u(2*Lx*\pi,y),
   ! u(x,y=0)=u(x,y=2*Ly*\pi)
   ! The initial condition is u=exp(-x^2-y^2)
   !
   ! .. 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
   ! plotgap = number of timesteps between plots
   ! 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
   ! Lx = width of box in x direction
   ! Ly = width of box in y direction
   ! ES = +1 for focusing and -1 for defocusing
   ! .. Scalars ..
   ! i = loop counter in x direction
   ! j = loop counter in y direction
   ! n = loop counter for timesteps direction
   ! allocatestatus = error indicator during allocation
   ! 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 = approximate solution
   ! v = Fourier transform of approximate solution
   ! unax = temporary field
   ! vnax = temporary field
   ! vnbx = temporary field
   ! vnay = temporary field
   ! vnby = temporary field
   ! potx = potential
   ! .. Vectors ..
   ! kx = fourier frequencies in x direction
   ! ky = 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 x Fourier transform
   ! fftfy = array to setup y Fourier transform
   ! fftby = array to setup y Fourier transform
   !
   ! REFERENCES
   !
   ! ACKNOWLEDGEMENTS
   !
   ! ACCURACY
   !
   ! ERROR INDICATORS AND WARNINGS
   !
   ! FURTHER COMMENTS
   ! Check that the initial iterate is consistent with the
   ! boundary conditions for the domain specified
   !--------------------------------------------------------------------
   ! External routines required
   !
   ! External libraries required
   ! FFTW3 -- Fast Fourier Transform in the West Library
   ! (http://www.fftw.org/)
   ! OpenMP library
   
   PROGRAM main
   USE omp_lib	
   IMPLICIT NONE	
   ! Declare variables
   INTEGER(kind=4), PARAMETER :: Nx=2**8	
   INTEGER(kind=4), PARAMETER :: Ny=2**8	
   INTEGER(kind=4), PARAMETER	:: Nt=20	
   INTEGER(kind=4), PARAMETER	:: plotgap=5	
   REAL(kind=8), PARAMETER	:: &
   pi=3.14159265358979323846264338327950288419716939937510d0
   REAL(kind=8), PARAMETER	:: Lx=2.0d0	
   REAL(kind=8), PARAMETER	:: Ly=2.0d0	
   REAL(kind=8), PARAMETER	:: Es=0.0d0	
   REAL(kind=8)	:: dt=0.10d0/Nt	
   COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx,ky
   REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x,y
   COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: unax,vnax,vnbx,potx
   COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: vnay,vnby
   REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: time
   INTEGER(kind=4)	:: i,j,k,n,allocatestatus
   INTEGER(kind=4)	:: start, finish, count_rate
   INTEGER(kind=8), PARAMETER :: FFTW_IN_PLACE=8, FFTW_MEASURE=0,&
   FFTW_EXHAUSTIVE=8, FFTW_PATIENT=32,&
                       FFTW_ESTIMATE=64
   INTEGER(kind=8),PARAMETER :: FFTW_FORWARD=-1, FFTW_BACKWARD=1	
   COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE :: fftfx,fftbx,fftfy,fftby
   INTEGER(kind=8)	:: planfx,planbx,planfy,planby
       CHARACTER*100	:: name_config
   
   ALLOCATE(kx(1:Nx),ky(1:Nx),x(1:Nx),y(1:Nx),unax(1:Nx,1:Ny),&
   vnax(1:Nx,1:Ny),vnbx(1:Nx,1:Ny),potx(1:Nx,1:Ny),fftfx(1:Nx),&
   fftbx(1:Nx),fftfy(1:Nx),fftby(1:Nx),vnay(1:Ny,1:Nx),&
   vnby(1:Ny,1:Nx),time(1:1+Nt/plotgap),stat=allocatestatus)	
   IF (allocatestatus .ne. 0) stop
   PRINT *,'allocated memory'
   ! set up ffts
   CALL dfftw_plan_dft_1d_(planfx,Nx,fftfx(1:Nx),fftbx(1:Nx),&
   FFTW_FORWARD,FFTW_ESTIMATE)
     CALL dfftw_plan_dft_1d_(planbx,Nx,fftbx(1:Nx),fftfx(1:Nx),&
     FFTW_BACKWARD,FFTW_ESTIMATE)
   CALL dfftw_plan_dft_1d_(planfy,Ny,fftfy(1:Ny),fftby(1:Ny),&
   FFTW_FORWARD,FFTW_ESTIMATE)
     CALL dfftw_plan_dft_1d_(planby,Ny,fftby(1:Ny),fftfy(1:Ny),&
     FFTW_BACKWARD,FFTW_ESTIMATE)	
   PRINT *,'Setup FFTs'
   
   ! setup fourier frequencies
   !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(static)
   DO i=1,1+Nx/2
   kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/Lx
   END DO
   !$OMP END PARALLEL DO
   kx(1+Nx/2)=0.0d0
   !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(static)
   DO i = 1,Nx/2 -1
   kx(i+1+Nx/2)=-kx(1-i+Nx/2)
   END DO
   !$OMP END PARALLEL DO
   !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(static)
      DO i=1,Nx
   x(i)=(-1.0d0+2.0d0*REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)) )*pi*Lx
   END DO
   !$OMP END PARALLEL DO
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
   DO j=1,1+Ny/2
   ky(j)= cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0))/Ly
   END DO
   !$OMP END PARALLEL DO
   ky(1+Ny/2)=0.0d0
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
   DO j = 1,Ny/2 -1
   ky(j+1+Ny/2)=-ky(1-j+Ny/2)
   END DO
   !$OMP END PARALLEL DO
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
      DO j=1,Ny
   y(j)=(-1.0d0+2.0d0*REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0)) )*pi*Ly
   END DO
   !$OMP END PARALLEL DO
   PRINT *,'Setup grid and fourier frequencies'
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
   DO j=1,Ny
   DO i=1,Nx
   unax(i,j)=exp(-1.0d0*(x(i)**2 +y(j)**2))
   END DO
   END DO
   !$OMP END PARALLEL DO
   name_config = 'uinitial.dat'
   OPEN(unit=11,FILE=name_config,status="UNKNOWN")
   REWIND(11)
   DO j=1,Ny
   DO i=1,Nx
   WRITE(11,*) abs(unax(i,j))**2
   END DO
   END DO
   CLOSE(11)
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
   DO j=1,Ny
   DO i=1,Nx
   CALL dfftw_execute_dft_(planfx,unax(i,j),vnax(i,j))
   END DO
   END DO
   !$OMP END PARALLEL DO
   vnay(1:Ny,1:Nx)=TRANSPOSE(vnax(1:Nx,1:Ny))
   ! transform initial data and do first half time step
   !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(static)
   DO i=1,Nx
   CALL dfftw_execute_dft_(planfy,vnay(1:Ny,i),vnby(1:Ny,i))
   DO j=1,Ny
   vnby(j,i)=exp(0.5d0*dt*(kx(i)*kx(i) + ky(j)*ky(j))&
   *cmplx(0.0d0,1.0d0))*vnby(j,i)
   END DO
   CALL dfftw_execute_dft_(planby,vnby(j,i),vnay(j,i))
   END DO
   !$OMP END PARALLEL DO
   PRINT *,'Got initial data, starting timestepping'
   time(1)=0.0d0
   CALL system_clock(start,count_rate)	
   DO n=1,Nt	
   vnbx(1:Nx,1:Ny)=TRANSPOSE(vnay(1:Ny,1:Nx))/REAL(Ny,kind(0d0))	
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
   DO j=1,Ny
   CALL dfftw_execute_dft_(planbx,vnbx(1:Nx,j),unax(1:Nx,j))
   DO i=1,Nx
   unax(i,j)=unax(1:Nx,j)/REAL(Nx,kind(0d0))
   potx(i,j)=Es*unax(i,j)*conjg(unax(i,j))
   unax(i,j)=exp(cmplx(0.0d0,-1.0d0)*dt*potx(i,j))&
   *unax(i,j)
   END DO
   CALL dfftw_execute_dft_(planfx,unax(1:Nx,j),vnax(1:Nx,j))
   END DO
   !$OMP END PARALLEL DO
   vnby(1:Ny,1:Nx)=TRANSPOSE(vnax(1:Nx,1:Ny))	
   !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(static)
   DO i=1,Nx
   CALL dfftw_execute_dft_(planfy,vnby(1:Ny,i),vnay(1:Ny,i))
   DO j=1,Ny
   vnby(j,i)=exp(dt*(kx(i)*kx(i) + ky(j)*ky(j))&	
   *cmplx(0.0d0,1.0d0))*vnay(j,i)	
   END DO
   CALL dfftw_execute_dft_(planby,vnby(1:Ny,i),vnay(1:Ny,i))
   END DO
   !$OMP END PARALLEL DO
   IF (mod(n,plotgap)==0) then
   time(1+n/plotgap)=n*dt
   PRINT *,'time',n*dt
   END IF
   END DO
   PRINT *,'Finished time stepping'
   CALL system_clock(finish,count_rate)
   PRINT*,'Program took ',REAL(finish-start)/REAL(count_rate),&
   'for Time stepping'
   
   ! transform back final data and do another half time step
   vnbx(1:Nx,1:Ny)=transpose(vnay(1:Ny,1:Nx))/REAL(Ny,kind(0d0))	
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
   DO j=1,Ny
   CALL dfftw_execute_dft_(planbx,vnbx(1:Nx,j),unax(1:Nx,j))
   unax(1:Nx,j)=unax(1:Nx,j)/REAL(Nx,kind(0d0))
   potx(1:Nx,j)=Es*unax(1:Nx,j)*conjg(unax(1:Nx,j))
   unax(1:Nx,j)=exp(cmplx(0,-1)*dt*potx(1:Nx,j))*unax(1:Nx,j)
   CALL dfftw_execute_dft_(planfx,unax(1:Nx,j),vnax(1:Nx,j))
   END DO
   !$OMP END PARALLEL DO
   vnby(1:Ny,1:Nx)=TRANSPOSE(vnax(1:Nx,1:Ny))	
   !$OMP PARALLEL DO PRIVATE(i) SCHEDULE(static)
   DO i=1,Nx
   CALL dfftw_execute_dft_(planfy,vnby(1:Ny,i),vnay(1:Ny,i))
   vnby(1:Ny,i)=exp(0.5d0*dt*(kx(i)*kx(i) + ky(1:Ny)*ky(1:Ny))&
   *cmplx(0,1))*vnay(1:Ny,i)	
   CALL dfftw_execute_dft_(planby,vnby(1:Ny,i),vnay(1:Ny,i))
   END DO
   !$OMP END PARALLEL DO
   vnbx(1:Nx,1:Ny)=TRANSPOSE(vnay(1:Ny,1:Nx))/REAL(Ny,kind(0d0))	
   !$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
   DO j=1,Ny
   CALL dfftw_execute_dft_(planbx,vnbx(1:Nx,j),unax(1:Nx,j))
   unax(1:Nx,j)=unax(1:Nx,j)/REAL(Nx,kind(0d0))
   END DO	
   !$OMP END PARALLEL DO
   name_config = 'ufinal.dat'
   OPEN(unit=11,FILE=name_config,status="UNKNOWN")
   REWIND(11)
   DO j=1,Ny
   DO i=1,Nx
   WRITE(11,*) abs(unax(i,j))**2
   END DO
   END DO
   CLOSE(11)
   
   name_config = 'tdata.dat'
   OPEN(unit=11,FILE=name_config,status="UNKNOWN")
   REWIND(11)
   DO j=1,1+Nt/plotgap
   WRITE(11,*) time(j)
   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)
   PRINT *,'Saved data'
   
   CALL dfftw_destroy_plan_(planbx)
   CALL dfftw_destroy_plan_(planfx)
   CALL dfftw_destroy_plan_(planby)
   CALL dfftw_destroy_plan_(planfy)
   CALL dfftw_cleanup_()
   
   DEALLOCATE(unax,vnax,vnbx,potx, vnay,vnby,stat=allocatestatus)	
   IF (allocatestatus .ne. 0) STOP
   PRINT *,'Deallocated memory'
   
   PRINT *,'Program execution complete'
   END PROGRAM main
   

   
   

   ( K)

   一个用于编译清单 J 中OpenMP程序的示例makefile 代码下载

   #define the complier
   COMPILER = gfortran
   # compilation settings, optimization, precision, parallelization
   FLAGS = -O0 -fopenmp
   
   
   # libraries
   LIBS = -L/usr/local/lib -lfftw3 -lm
   # source list for main program
   SOURCES = NLSsplitting.f90
   
   test: $(SOURCES)
   ${COMPILER} -o NLSsplitting $(FLAGS) $(SOURCES) $(LIBS)
   
   clean:
   rm *.o
   
   clobber:
   rm NLSsplitting
   

   此示例假设您使用的是Flux，并且已加载了intel编译器以及Intel编译版本的FFTW的环境。如果您使用的是免费的GCC编译器而不是Intel编译器，请将makefile中的标志 $-openmp$ 更改为 $-fopenmp$。

6. 修改OpenMP Fortran代码以求解二维和三维的Gross-Pitaevskii方程。

7. 一些等离子体的量子流体动力学模型与非线性薛定谔方程非常相似，也可以使用分裂方法进行数值近似。Eliasson和Shukla 使用的一个等离子体模型是

   $i\Psi _{t}+\Delta \Psi +\phi \Psi -|\Psi |^{4/D}\Psi =0$

   ( 8)

   和

   $\Delta \phi =|\Psi |^{2}-1$

   ( 9)

   其中 $\Psi$ 是有效波函数，$\phi$ 是静电势，$D$ 是维度，通常为 1、2 或 3。此方程可以以与Klein、Muite和Roidot 中的Davey-Stewartson方程类似的方式求解。具体来说，首先求解

   $iP_{t}+\Delta P=0$

   ( 10)

   使用傅里叶变换，使得

   $P(\delta t)=\exp \left(-i\Delta \delta t\right)P(0)$

   其中 $P(0)=\Psi (0)$。然后求解

   $\phi =\Delta ^{-1}\left(|P|^{2}-1\right)$

   ( 11)

   使用傅里叶变换。最后，求解

   $iQ_{t}+\phi Q-|Q|^{4/D}Q=0$

   ( 12)

   利用在每个网格点上$\phi -|Q|^{4/D}$ 是一个常数，所以解是 $Q(\delta t)=\exp \left[i\left(\phi -|\Phi |^{4/D}\right)\delta t\right]Q(0)$ 其中$Q(0)=P(\delta t)$ 以及$\Psi (\delta t)\approx Q(\delta t)$.

8. 算子分裂法可以用于除非线性薛定谔方程以外的其他方程。另一个可以使用算子分裂法的方程是复 Ginzburg-Landau 方程 ${\frac {\partial A}{\partial t}}=A+(1+i\alpha )\Delta A-(1+i\beta )|A|^{2}A,$ 其中 $A$ 是一个复函数，通常是一个、两个或三个变量的函数。基于 Blanes、Casa、Chartier 和 Miura 早期有限差分代码，列表 L 提供了一个一维代码示例，该示例使用了 Blanes 等人 描述的方法。通过使用复系数，Blanes 等人 可以为抛物线方程创建高阶分裂方法。以前的尝试失败了，因为如果只使用实系数，则对于高于二阶的方法，所需的向后步骤会导致数值不稳定。修改示例代码，以解决一维、二维和三维空间中的复 Ginzburg-Landau 方程。线性部分 ${\frac {\partial A}{\partial t}}=A+(1+i\alpha )\Delta A$ 可以使用傅里叶变换显式求解。要解决非线性部分，${\frac {\partial A}{\partial t}}=-(1+i\beta )|A|^{2}A$ 考虑 ${\frac {\partial |A|^{2}}{\partial t}}={\frac {\partial A}{\partial t}}A^{*}+{\frac {\partial A^{*}}{\partial t}}A=2|A|^{4}$ 并精确地求解 $|A|^{2}$。为了恢复相位，观察 ${\frac {\partial \log(A)}{\partial t}}=-(1+i\beta )|A|^{2}$，它也可以显式积分，因为 $|A|^{2}(t)$ 是已知的。

   
   

   ( L)

   一个使用 16 阶分裂法来解决三次非线性薛定谔方程的 Matlab 程序 代码下载

   % A program to solve the nonlinear Schr\"{o}dinger equation using a
   % splitting method. The numerical solution is compared to an exact
   % solution.
   % S. Blanes, F. Casas, P. Chartier and A. Murua
   % "Optimized high-order splitting methods for some classes of parabolic
   % equations"
   % ArXiv pre-print 1102.1622v2
   % Forthcoming Mathematics of Computation
   
   clear all; format compact; format short;
   set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
       'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
       'defaultaxesfontweight','bold')
   
   % set up grid
   Lx = 20; % period 2*pi * L
   Nx = 16384; % number of harmonics
   Nt = 2000; % number of time slices
   dt = 0.25*pi/Nt;% time step
   U=zeros(Nx,Nt/10);
   method=3; % splitting method: 1 Strang, 2 CCDV10, 3 Blanes et al 2012
   
   % 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
   
   % initial conditions
   t=0; tdata(1)=t;
   u=4*exp(1i*t)*(cosh(3*x)+3*exp(8*1i*t)*cosh(x))...
       ./(cosh(4*x)+4*cosh(2*x)+3*cos(8*t));
   v=fft(u);
   figure(1); clf; plot(x,u);xlim([-2,2]); drawnow;
   U(:,1)=u;
   
   % mass
   ma = fft(abs(u).^2);
   ma0 = ma(1);
   
   if method==1,
       %
       % Strang-Splitting
       %
       s=2;
       a=[1;0];
       b=[1/2;1/2];
       %
   elseif method==2,
       %
       % Method of Castella, Chartier, Descombes and Vilmart
       % BIT Numerical Analysis vol 49 pp 487-508, 2009
       %
       s=5;
       a=[1/4;1/4;1/4;1/4;0];
       b=[1/10-1i/30;4/15+2*1i/15;4/15-1i/5;4/15+2*1i/15;1/10-1i/30];
       %
   elseif method==3,
       %
       % Method of Blanes, Casas, Chartier and Murua 2012
       %
       s=17;
       a=1/16*[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0];
       b=[0.028920177910074098791 - 0.005936580835725746103*1i;
          0.056654351383649876160 + 0.020841963949772627119*1i;
          0.067258385822722143569 - 0.039386393748812362460*1i;
          0.070333980553260772061 + 0.058952097930307840316*1i;
          0.077095100838099173580 - 0.038247636602014810025*1i;
          0.042022140317231098258 - 0.033116379859951038579*1i;
          0.050147397749937784280 + 0.061283684958324249562*1i;
          0.047750191909146143447 - 0.032332468814362628262*1i;
          0.119636547031757819706 + 0.015883426044923736862*1i;
          0.047750191909146143447 - 0.032332468814362628262*1i;
          0.050147397749937784280 + 0.061283684958324249562*1i;
          0.042022140317231098258 - 0.033116379859951038579*1i;
          0.077095100838099173580 - 0.038247636602014810025*1i;
          0.070333980553260772061 + 0.058952097930307840316*1i;
          0.067258385822722143569 - 0.039386393748812362460*1i;
          0.056654351383649876160 + 0.020841963949772627119*1i;
          0.028920177910074098791 - 0.005936580835725746103*1i];
   end;
   
   
   % solve pde and plot results
   for n =2:Nt+1
       for m=1:(s-1)
           vna=exp(b(m)*1i*dt*kx.*kx).*v;
           una=ifft(vna);
           pot=(2*una.*conj(una));
           unb=exp(-1i*a(m)*(-1)*dt*pot).*una;
           v=fft(unb);
       end
       v=exp(b(s)*1i*dt*kx.*kx).*v;
       u=ifft(v);
       t=(n-1)*dt;
       if (mod(n,10)==0)
           tdata(n/10)=t;
           u=ifft(v);
           U(:,n/10)=u;
           uexact=...
               4*exp(1i*t)*(cosh(3*x)+3*exp(8*1i*t)*cosh(x))...
               ./(cosh(4*x)+4*cosh(2*x)+3*cos(8*t));
           figure(1); clf; plot(x,abs(u).^2); ...
               xlim([-0.5,0.5]); title(num2str(t));
           figure(2); clf; loglog(abs(v(1:Nx/2))); ...
               title('Fourier Coefficients');
           figure(3); clf; plot(x,abs(u-uexact).^2); ...
               xlim([-0.5,0.5]); title('error');
           drawnow;
           ma = fft(abs(u).^2);
           ma = ma(1);
           test = log10(abs(1-ma/ma0))
       end
   end
   figure(4); clf; mesh(tdata(1:(n-1)/10),x,abs(U(:,1:(n-1)/10)).^2);
   xlim([0,t]);
   

分布式内存并行：MPI

在本节中，我们将使用从 http://www.2decomp.org/index.html 获得的 2DECOMP&FFT 库。该网站包含一些示例，说明了如何使用此库，特别是 http://www.2decomp.org/case_study1.html 中的示例代码，非常有助于说明如何将使用 FFTW 的代码转换为使用 MPI 和上述库的代码。

在使用 2DECOMP&FFT 创建并行 MPI 代码之前，我们将生成一个使用分裂法来解决 3D 非线性薛定谔方程的串行 Fortran 代码。我们不会使用基于循环的并行化来执行一系列一维快速傅里叶变换，而是使用 FFTW 的三维 FFT，这样串行版本和 MPI 并行版本就具有相同的结构。串行版本在列表 E 中。此文件可以通过与 Fortran 程序章节中的列表 A 中相同的方式进行编译。

( M)

一个使用分裂法和 FFTW 解决 3D 非线性薛定谔方程的 Fortran 程序 代码下载

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program solves nonlinear Schrodinger equation in 3 dimensions
! i*u_t+Es*|u|^2u+u_{xx}+u_{yy}+u_{zz}=0
! using a second order time spectral splitting scheme
!
! The boundary conditions are u(x=0,y,z)=u(2*Lx*\pi,y,z),
! u(x,y=0,z)=u(x,y=2*Ly*\pi,z), u(x,y,z=0)=u(x,y,z=2*Lz*\pi)
! The initial condition is u=exp(-x^2-y^2)
!
! .. 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
! plotgap = number of timesteps between plots
! 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
! Lx = width of box in x direction
! Ly = width of box in y direction
! Lz = width of box in z direction
! ES = +1 for focusing and -1 for defocusing
! .. 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
! start = variable to record start time of program
! finish = variable to record end time of program
! count_rate = variable for clock count rate
! count = keep track of information written to disk
! iol = size of array to write to disk
! planfxyz = Forward 3d fft plan
! planbxyz = Backward 3d fft plan
! dt = timestep
! modescalereal = Number to scale after backward FFT
! ierr = error code
! .. Arrays ..
! unax = approximate solution
! vnax = Fourier transform of approximate solution
! potx = potential
! .. Vectors ..
! kx = fourier frequencies in x direction
! ky = 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
! fftfxy = array to setup 2D Fourier transform
! fftbxy = array to setup 2D Fourier transform
!
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
! Check that the initial iterate is consistent with the
! boundary conditions for the domain specified
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! FFTW3 -- Fast Fourier Transform in the West Library
! (http://www.fftw.org/)
PROGRAM main
! Declare variables
IMPLICIT NONE
INTEGER(kind=4), PARAMETER	:: Nx=2**5	
INTEGER(kind=4), PARAMETER :: Ny=2**5	
INTEGER(kind=4), PARAMETER :: Nz=2**5	
INTEGER(kind=4), PARAMETER	:: Nt=50	
INTEGER(kind=4), PARAMETER	:: plotgap=10
REAL(kind=8), PARAMETER	::&
pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8), PARAMETER	:: Lx=2.0d0,Ly=2.0d0,Lz=2.0d0	
REAL(kind=8), PARAMETER	:: Es=1.0d0	
REAL(kind=8)	:: dt=0.10d0/Nt	
REAL(kind=8) :: modescalereal
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx,ky,kz
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x,y,z
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE	:: unax,vnax,potx
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: time
INTEGER(kind=4)	:: i,j,k,n,AllocateStatus,count,iol
! timing
INTEGER(kind=4)	:: start, finish, count_rate
   ! fftw variables
INTEGER(kind=8), PARAMETER	:: FFTW_IN_PLACE = 8, FFTW_MEASURE = 0, &
FFTW_EXHAUSTIVE = 8, FFTW_PATIENT = 32, FFTW_ESTIMATE = 64
INTEGER(kind=8),PARAMETER	:: FFTW_FORWARD = -1, FFTW_BACKWARD=1	
INTEGER(kind=8)	::	planfxyz,planbxyz
CHARACTER*100	::	name_config, number_file

CALL system_clock(start,count_rate)	
ALLOCATE(unax(1:Nx,1:Ny,1:Nz),vnax(1:Nx,1:Ny,1:Nz),potx(1:Nx,1:Ny,1:Nz),&
    kx(1:Nx),ky(1:Ny),kz(1:Nz),x(1:Nx),y(1:Ny),z(1:Nz),&
    time(1:1+Nt/plotgap),stat=AllocateStatus)	
IF (AllocateStatus .ne. 0) STOP
PRINT *,'allocated space'
modescalereal=1.0d0/REAL(Nx,KIND(0d0))
modescalereal=modescalereal/REAL(Ny,KIND(0d0))
modescalereal=modescalereal/REAL(Nz,KIND(0d0))

! set up ffts
CALL dfftw_plan_dft_3d_(planfxyz,Nx,Ny,Nz,unax(1:Nx,1:Ny,1:Nz),&
vnax(1:Nx,1:Ny,1:Nz),FFTW_FORWARD,FFTW_ESTIMATE)
  CALL dfftw_plan_dft_3d_(planbxyz,Nx,Ny,Nz,vnax(1:Nx,1:Ny,1:Nz),&
  unax(1:Nx,1:Ny,1:Nz),FFTW_BACKWARD,FFTW_ESTIMATE)

PRINT *,'Setup FFTs'

! setup fourier frequencies and grid points
DO i=1,1+Nx/2
kx(i)= cmplx(0.0d0,1.0)*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
DO i=1,Nx
x(i)=(-1.0d0+2.0d0*REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)) )*pi*Lx
END DO
DO j=1,1+Ny/2
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
DO j=1,Ny
y(j)=(-1.0d0+2.0d0*REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0)) )*pi*Ly
END DO
DO k=1,1+Nz/2
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
DO k=1,Nz
z(k)=(-1.0d0+2.0d0*REAL(k-1,kind(0d0))/REAL(Nz,kind(0d0)) )*pi*Lz
END DO

PRINT *,'Setup grid and fourier frequencies'

DO k=1,Nz;	DO j=1,Ny; DO i=1,Nx
unax(i,j,k)=exp(-1.0d0*(x(i)**2 +y(j)**2+z(k)**2))
END DO; END DO; END DO
name_config = 'uinitial.dat'
INQUIRE(iolength=iol) unax(1,1,1)
OPEN(unit=11,FILE=name_config,form="unformatted", &
access="direct",recl=iol)
count=1
DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
WRITE(11,rec=count) unax(i,j,k)
        count=count+1
END DO; END DO; END DO
CLOSE(11)

CALL dfftw_execute_dft_(planfxyz,unax(1:Nx,1:Ny,1:Nz),vnax(1:Nx,1:Ny,1:Nz))

PRINT *,'Got initial data, starting timestepping'
time(1)=0
DO n=1,Nt
DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
vnax(i,j,k)=exp(0.50d0*dt*&
(kz(k)*kz(k) + kx(i)*kx(i) + ky(j)*ky(j))&
*cmplx(0.0d0,1.0d0))*vnax(i,j,k)
END DO; END DO; END DO
CALL dfftw_execute_dft_(planbxyz,vnax(1:Nx,1:Ny,1:Nz),&
unax(1:Nx,1:Ny,1:Nz))

DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
unax(i,j,k)=unax(i,j,k)*modescalereal
potx(i,j,k)=Es*unax(i,j,k)*conjg(unax(i,j,k))
unax(i,j,k)=exp(cmplx(0.0d0,-1.0d0)*dt*potx(i,j,k))&
*unax(i,j,k)
END DO; END DO; END DO
CALL dfftw_execute_dft_(planfxyz,unax(1:Nx,1:Ny,1:Nz),&
vnax(1:Nx,1:Ny,1:Nz))

DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
vnax(i,j,k)=exp(0.5d0*dt*&
(kx(i)*kx(i) + ky(j)*ky(j)+ kz(k)*kz(k))&
*cmplx(0.0d0,1.0d0))*vnax(i,j,k)	
END DO; END DO; END DO
IF (mod(n,plotgap)==0) THEN
time(1+n/plotgap)=n*dt
PRINT *,'time',n*dt
CALL dfftw_execute_dft_(planbxyz,vnax(1:Nx,1:Ny,1:Nz),unax(1:Nx,1:Ny,1:Nz))
DO k=1,Nz; DO j=1,Ny; DO i=1,Nx
unax(i,j,k)=unax(i,j,k)*modescalereal
END DO; END DO; END DO
name_config='./data/u'
WRITE(number_file,'(i0)') 10000000+1+n/plotgap
ind=index(name_config,' ') -1
name_config=name_config(1:ind)//numberfile
ind=index(name_config,' ') -1
name_config=name_config(1:ind)//'.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Ny
DO i=1,Nx
WRITE(11,*) abs(unax(i,j))**2
END DO
END DO
CLOSE(11)

END IF
END DO
PRINT *,'Finished time stepping'

! transform back final data and do another half time step
CALL system_clock(finish,count_rate)
PRINT*,'Program took ',REAL(finish-start)/REAL(count_rate),'for execution'

name_config = 'tdata.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt/plotgap
WRITE(11,*) time(j)
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)

name_config = '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'

CALL dfftw_destroy_plan_(planbxyz)
CALL dfftw_destroy_plan_(planfxyz)
CALL dfftw_cleanup_()

DEALLOCATE(unax,vnax,potx,&
    kx,ky,kz,x,y,z,&
    time,stat=AllocateStatus)	
IF (AllocateStatus .ne. 0) STOP
PRINT *,'Program execution complete'
END PROGRAM main

与之前的程序相比，列表 E 中的程序将其最终数据作为二进制文件写入。这通常比写入文本文件快得多，并且生成的二进制文件的大小通常也小得多。当写入多个此类文件时，或者当单个文件很大时，这一点很重要。由于格式发生了变化，二进制文件也需要使用稍微不同的方式读入。列表 N 中的 Matlab 脚本展示了如何做到这一点。

( N)

一个绘制列表 E 或 O 生成的 3D 非线性薛定谔方程数值解的 Matlab 程序 代码下载

% A program to plot the computed results

clear all; format compact, format short,
set(0,'defaultaxesfontsize',18,'defaultaxeslinewidth',.9,...
    'defaultlinelinewidth',3.5,'defaultpatchlinewidth',5.5);

% Load data
tdata=load('./tdata.dat');
x=load('./xcoord.dat');
y=load('./ycoord.dat');
z=load('./zcoord.dat');
Tsteps = length(tdata);

Nx = length(x); Nt = length(tdata);
Ny = length(y); Nz = length(z);
fid=fopen('./ufinal.datbin','r');
[fname,mode,mformat]=fopen(fid);
u=fread(fid,Nx*Ny*Nz,'double',mformat);
u = reshape(u,Nx,Ny,Nz);

% Plot data
figure (1); clf ; UP = abs(u).^2;
p1 = patch(isosurface(x,y,z,UP,.0025) ,...
'FaceColor','yellow','EdgeColor','none');
p2 = patch(isocaps(x,y,z,UP,.0025) ,...
'FaceColor','interp','EdgeColor','none');
isonormals(UP,p1); lighting phong;
xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); drawnow;

现在，我们将修改上述代码，以使用 MPI 和 2DECOMP&FFT 库。2DECOMP&FFT 库隐藏了大多数 MPI 的细节，尽管有一些命令对用户理解非常有用。这些命令是

• USE mpi或INCLUDE 'mpif.h'
• MPI_INIT
• MPI_COMM_SIZE
• MPI_COMM_RANK
• MPI_FINALIZE

程序在清单 O 中列出，请将其与 E 中的串行代码进行比较。库 2DECOMP&FFT 对数组进行域分解，以便将数组的不同部分放到不同的处理器上。即使数组存储在不同的处理器上，该库也可以对数组执行傅里叶变换——该库执行执行傅里叶变换所需的所有必要消息传递和转置。需要注意的是，傅里叶变换后数组中条目顺序不一定与 FFTW 使用的顺序相同。但是，结构将返回条目的正确排序decomp因此，此结构用于获取循环的起始和结束条目。我们假设库 2DECOMP&FFT 已安装在适当的位置。

( O)

使用分裂和 2DECOMP&FFT 求解 3D 非线性薛定谔方程的 Fortran 程序 代码下载

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program solves nonlinear Schrodinger equation in 3 dimensions
! i*u_t+Es*|u|^2u+u_{xx}+u_{yy}+u_{zz}=0
! using a second order time spectral splitting scheme
!
! The boundary conditions are u(x=0,y,z)=u(2*Lx*\pi,y,z),
! u(x,y=0,z)=u(x,y=2*Ly*\pi,z), u(x,y,z=0)=u(x,y,z=2*Lz*\pi)
! The initial condition is u=exp(-x^2-y^2)
!
! .. 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
! plotgap = number of timesteps between plots
! pi = 3.14159265358979323846264338327950288419716939937510d0
! Lx = width of box in x direction
! Ly = width of box in y direction
! Lz = width of box in z direction
! ES = +1 for focusing and -1 for defocusing
! .. 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
! start = variable to record start time of program
! finish = variable to record end time of program
! count_rate = variable for clock count rate
! dt = timestep
! modescalereal = Number to scale after backward FFT
! myid = Process id
! ierr = error code
! p_row = number of rows for domain decomposition
! p_col = number of columns for domain decomposition
! filesize = total filesize
! disp = displacement to start writing data from
! ind = index in array to write
! plotnum = number of plot to save
! numberfile = number of the file to be saved to disk
! stat = error indicator when reading inputfile
! .. Arrays ..
! u = approximate solution
! v = Fourier transform of approximate solution
! pot = potential
! .. 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 = z locations
! time = times at which save data
! nameconfig = array to store filename for data to be saved
! InputFileName = name of the Input File
! .. Special Structures ..
! decomp = contains information on domain decomposition
! see http://www.2decomp.org/ for more information
!
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
! Check that the initial iterate is consistent with the
! boundary conditions for the domain specified
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
! 2DECOMP&FFT -- Domain decomposition and Fast Fourier Library
! (http://www.2decomp.org/index.html)
! MPI library

PROGRAM main
USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
USE MPI
! Declare variables
IMPLICIT NONE
INTEGER(kind=4)	:: Nx=2**5
INTEGER(kind=4) :: Ny=2**5
INTEGER(kind=4) :: Nz=2**5
INTEGER(kind=4)	:: Nt=50	
INTEGER(kind=4)	:: plotgap=10
REAL(kind=8), PARAMETER	::&
pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8)	:: Lx=2.0d0,Ly=2.0d0,Lz=2.0d0	
REAL(kind=8)	:: Es=1.0d0	
REAL(kind=8)	:: dt=0.0010d0	
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx,ky,kz
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x,y,z
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE	:: u,v,pot
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: time
INTEGER(KIND=4), DIMENSION(1:5)	:: intcomm
REAL(KIND=8), DIMENSION(1:5)	:: dpcomm
REAL(kind=8) :: modescalereal
INTEGER(kind=4)	:: i,j,k,n,AllocateStatus,stat
INTEGER(kind=4)	:: myid,numprocs,ierr
TYPE(DECOMP_INFO)	:: decomp
INTEGER(kind=MPI_OFFSET_KIND) :: filesize, disp
INTEGER(kind=4)	:: p_row=0, p_col=0	
  INTEGER(kind=4)	:: start, finish, count_rate, ind, plotnum
CHARACTER*50	::	nameconfig
CHARACTER*20	::	numberfile
! initialisation of MPI
CALL MPI_INIT(ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)

! initialisation of 2decomp
! 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
ALLOCATE(u(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
            v(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
            pot(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),&
    x(1:Nx),y(1:Ny),z(1:Nz),&
    time(1:1+Nt/plotgap),stat=AllocateStatus)	
IF (AllocateStatus .ne. 0) STOP

IF (myid.eq.0) THEN
PRINT *,'allocated space'
END IF
modescalereal=1.0d0/REAL(Nx,KIND(0d0))
modescalereal=modescalereal/REAL(Ny,KIND(0d0))
modescalereal=modescalereal/REAL(Nz,KIND(0d0))

! setup fourier frequencies and grid points
DO i=1,1+Nx/2
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
DO i=1,Nx
x(i)=(-1.0d0 + 2.0d0*REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)))*pi*Lx
END DO
DO j=1,1+Ny/2
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
DO j=1,Ny
y(j)=(-1.0d0 + 2.0d0*REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0)))*pi*Ly
END DO
DO k=1,1+Nz/2
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
DO k=1,Nz
z(k)=(-1.0d0 + 2.0d0*REAL(k-1,kind(0d0))/REAL(Nz,kind(0d0)))*pi*Lz
END DO

IF (myid.eq.0) THEN
PRINT *,'Setup grid and fourier frequencies'
END IF
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)=exp(-1.0d0*(x(i)**2 +y(j)**2+z(k)**2))
END DO
END DO
END DO

! write out using 2DECOMP&FFT MPI-IO routines
nameconfig='./data/u'
plotnum=0
WRITE(numberfile,'(i0)') 10000000+plotnum
ind=index(nameconfig,' ') -1
nameconfig=nameconfig(1:ind)//numberfile
ind=index(nameconfig,' ') -1
nameconfig=nameconfig(1:ind)//'.datbin'
CALL decomp_2d_write_one(1,u,nameconfig)
  
CALL decomp_2d_fft_3d(u,v,DECOMP_2D_FFT_FORWARD)
IF (myid.eq.0) THEN
PRINT *,'Got initial data, starting timestepping'
END IF
CALL system_clock(start,count_rate)	
time(1)=0
DO n=1,Nt
! Use Strang splitting
DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
v(i,j,k)=exp(0.50d0*dt*&
(kz(k)*kz(k) + kx(i)*kx(i) + ky(j)*ky(j))&
*cmplx(0.0d0,1.0d0))*v(i,j,k)
END DO
END DO
END DO

CALL decomp_2d_fft_3d(v,u,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)*modescalereal
pot(i,j,k)=Es*u(i,j,k)*conjg(u(i,j,k))
u(i,j,k)=exp(cmplx(0.0d0,-1.0d0)*dt*pot(i,j,k))*u(i,j,k)
END DO
END DO
END DO
CALL decomp_2d_fft_3d(u,v,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)
v(i,j,k)=exp(dt*0.5d0*&
(kx(i)*kx(i) +ky(j)*ky(j) +kz(k)*kz(k))&
*cmplx(0.0d0,1.0d0))*v(i,j,k)	
END DO
END DO
END DO
IF (mod(n,plotgap)==0) THEN
time(1+n/plotgap)=n*dt
IF (myid.eq.0) THEN
PRINT *,'time',n*dt
END IF
CALL decomp_2d_fft_3d(v,u,DECOMP_2D_FFT_BACKWARD)
u=u*modescalereal
nameconfig='./data/u'
plotnum=plotnum+1
WRITE(numberfile,'(i0)') 10000000+plotnum
ind=index(nameconfig,' ') -1
nameconfig=nameconfig(1:ind)//numberfile
ind=index(nameconfig,' ') -1
nameconfig=nameconfig(1:ind)//'.datbin'
! write out using 2DECOMP&FFT MPI-IO routines
CALL decomp_2d_write_one(1,u,nameconfig)
END IF
END DO
IF (myid.eq.0) THEN
PRINT *,'Finished time stepping'
END IF

CALL system_clock(finish,count_rate)

IF (myid.eq.0) THEN
PRINT*,'Program took ',REAL(finish-start)/REAL(count_rate),'for execution'
END IF
IF (myid.eq.0) THEN	
! Save times at which output was made in text format
nameconfig = './data/tdata.dat'
OPEN(unit=11,FILE=nameconfig,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt/plotgap
WRITE(11,*) time(j)
END DO
CLOSE(11)
! Save x grid points in text format
nameconfig = './data/xcoord.dat'
OPEN(unit=11,FILE=nameconfig,status="UNKNOWN")
REWIND(11)
DO i=1,Nx
WRITE(11,*) x(i)
END DO
CLOSE(11)
! Save y grid points in text format
nameconfig = './data/ycoord.dat'
OPEN(unit=11,FILE=nameconfig,status="UNKNOWN")
REWIND(11)
DO j=1,Ny
WRITE(11,*) y(j)
END DO
CLOSE(11)
! Save z grid points in text format
nameconfig = './data/zcoord.dat'
OPEN(unit=11,FILE=nameconfig,status="UNKNOWN")
REWIND(11)
DO k=1,Nz
WRITE(11,*) z(k)
END DO
CLOSE(11)
PRINT *,'Saved data'
END IF

! clean up
   CALL decomp_2d_fft_finalize
   CALL decomp_2d_finalize
  DEALLOCATE(u,v,pot,&
    kx,ky,kz,x,y,z,&
    time,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

要编译代码，请根据您的系统在清单 P 中相应地修改 makefile。键入 make，然后在编译后您应该就可以运行它了。假设 FFTW 是 2DECOMP&FFT 调用的底层一维快速傅里叶变换库，因此如果使用其他库，则应更改它。请注意，如果您要更改问题大小或运行时间，则需要重新构建程序。

( P)

使用分裂和 2DECOMP&FFT 求解 3D 非线性薛定谔方程的并行 Fortran 程序的示例 makefile 代码下载

COMPILER =  mpif90 
decompdir=../2decomp_fft
FLAGS = -O0 
	
DECOMPLIB = -I${decompdir}/include -L${decompdir}/lib -l2decomp_fft
LIBS = -L${FFTW_LINK} -lfftw3_threads -lfftw3 -lm 
SOURCES = NLSsplitting.f90  

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


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

练习

1. ASCII 字符集需要每个字符 7 位，因此至少需要 7 位来存储 0 到 9 之间的数字。IEEE 算术中的双精度数需要 64 位来存储一个双精度数，大约有 15 位十进制数和大约 3 位十进制数指数。存储 IEEE 双精度数需要多少位？假设一个文件有 $10^{6}$ 个双精度数。如果数字存储为 IEEE 双精度数，则文件的最小大小是多少？如果数字存储为字符，则文件的最小大小是多少？
2. 使用 2DECOMP&FFT 编写一个 MPI 代码来求解三维 Gross-Pitaevskii 方程。
3. 学习使用 VisIt (https://wci.llnl.gov/codes/visit/) 或 Paraview (http://www.paraview.org/) 并编写一个脚本以类似于 Matlab 代码的方式可视化二维和三维输出。

参考

Blanes, S.; Casas, F.; Chartier, P.; Murua, A. "Splitting methods with complex coefficients for some classes of evolution equations". Mathematics of Computation.

Shukla, P.K. (2009). "Nonlinear aspects of quantum plasma physics: Nanoplasmonics and nanostructures in dense plasmas". Plasma and Fusion Research: Review Articles. 4: 32.

Evans, L.C. (2010). Partial Differential Equations. American Mathematical Society.

Holden, H.; Karlsen, K.H.; Risebro, N.H.; Tao, T. (2011). "Operator splitting for the KdV equation". Mathematics of Computation. 80: 821–846.

Klein, C. (2008). "Fourth order time-stepping for low dispersion Korteweg-De Vries and nonlinear Schrödinger equations". Electronic Transactions on Numerical Analysis. 29: 116–135.

Klein, C.; Muite, B.K.; Roidot, K. (2011). "Numerical Study of Blowup in the Davey-Stewartson System". {{cite journal}}: Cite journal requires |journal= (help)

Klein, C.; Roidot, K. (2011). "Fourth order time-stepping for Kadomstev-Petviashvili and Davey-Stewartson Equations". SIAM Journal on Scientific Computation. 33: 3333–3356.

Loss, M. (2003). Analysis. American Mathematical Society.

Ponce, G. (2009). Introduction to Nonlinear Dispersive Equations. Springer.

McLachlan, R.I.; Quispel, G.R.W. (2002). "Splitting Methods". Acta Numerica. 11: 341–434.

Rogers, R.C. (2004). An Introduction to Partial Differential Equations. Springer.

Wang, L.-L. (2011). Spectral Methods: Algorithms, Analysis and Applications. Springer.

Sulem, P.L. (1999). The Nonlinear Schrödinger equation: Self-Focusing and Wave Collapse. Springer.

Thalhammer, M. Time-Splitting Spectral Methods for Nonlinear Schrödinger Equations (PDF). {{cite book}}: Cite has empty unknown parameter: |coauthors= (help)

Weideman, J.A.C.; Herbst, B.M. (1986). "Split-step methods for the solution of the nonlinear Schrödinger equation". SIAM Journal on Numerical Analysis. SIAM. 23 (3): 485–507.

Yang, J. (2010). Nonlinear Waves in Integrable and Nonintegrable Systems. SIAM.