克莱因-戈登方程

背景

^[1]聚焦/散焦非线性克莱因-戈登方程描述了可能复标量场 $u$ 的演化，根据：

{\frac {\partial ^{2}u}{\partial t^{2}}}-\Delta u+u=\pm |u|^{2}u

,

( 1)

其中 $+$ 是聚焦情况， $-$ 是散焦情况，类似于非线性薛定谔方程。 Donninger 和 Schlag^[2]最近对该方程的三维径向对称实解的爆破进行了数值研究。克莱因-戈登方程的二维模拟可以在 Yang^[3] 中找到。线性克莱因-戈登方程作为线性薛定谔方程的修正出现，与狭义相对论一致，例如参见 Landau^[4] 或 Grenier^[5]。目前，还没有对该方程解的爆破进行数值研究，而没有假设径向对称性。该方程产生了大量的数学文献，但仍然知之甚少。这些数学文献的大部分集中在分析在无限三维空间上的方程，其初始数据要么随着趋于无穷大而呈指数衰减，要么在一个有限的域集中不为零。在这里，我们将对方程在一个周期性环境中进行模拟。由于该方程是一个波动方程，它具有有限的信息传播速度，就像空气中的声波需要时间从一点移动到另一点一样。因此，对于短时间模拟，对仅在域的有限部分上不为零的解的模拟类似于对无限域上的模拟。然而，在很长一段时间内，解可以在周期性域上扩散并与自身相互作用，而在无限域上，长时间的相互作用会大大减少，解主要扩散出去。理解周期性环境中的相互作用是一个有趣的数学问题。克莱因-戈登方程有一个守恒的能量，由下式给出：

$\int {\frac {1}{2}}\left({\frac {\partial u}{\partial t}}\right)^{2}+{\frac {u^{2}}{2}}+{\frac {1}{2}}\left|\nabla u\right|^{2}\mp {\frac {\left|u\right|^{4}}{4}}\mathrm {d} \mathbf {x}$ .

该方程也是时间可逆的。对于长时间模拟，人们希望构建近似守恒该能量并且也是时间可逆的数值方法。当使用傅里叶谱方法时，我们主要需要确保时间离散化保留了这些属性，因为谱空间离散化通常会自动满足这些属性。遵循 Donninger 和 Schlag^[6]，我们使用两种方案。首先，一个隐式-显式时间步进方案，该方案是时间可逆的，但仅近似守恒能量，并由下式给出：

{\frac {u^{n+1}-2u^{n}+u^{n-1}}{(\delta t)^{2}}}-\Delta {\frac {u^{n+1}+2u^{n}+u^{n-1}}{4}}+{\frac {u^{n+1}+2u^{n}+u^{n-1}}{4}}=\pm \left|u^{n}\right|^{2}u^{n}

( 2)

其次，一个带有固定点迭代的完全隐式时间步进方案

{\frac {u^{n+1,k+1}-2u^{n}+u^{n-1}}{(\delta t)^{2}}}-\Delta {\frac {u^{n+1,k+1}+2u^{n}+u^{n-1}}{4}}+{\frac {u^{n+1,k+1}+2u^{n}+u^{n-1}}{4}}

=\pm {\frac {\left|u^{n+1,k}\right|^{4}-\left|u^{n-1}\right|^{4}}{4(u^{n+1,k}-u^{n-1})}}

( 3)

它准确地保留了离散能量

\int {\frac {1}{2}}\left({\frac {u^{n+1}-u^{n}}{\delta t}}\right)^{2}+{\frac {1}{2}}\left({\frac {u^{n+1}+u^{n}}{2}}\right)^{2}+{\frac {1}{2}}\left|\nabla {\frac {u^{n+1}+u^{n}}{2}}\right|^{2}\mp {\frac {\left|{u}^{n+1}\right|^{4}+\left|{u}^{n}\right|^{4}}{8}}

.

( 4)

如前所述，上标 $n$ 表示时间步长， $k$ 表示不动点迭代方案中的迭代次数。当两个连续迭代之间的差异小于某个容差时，迭代停止。

Matlab 和 Python 程序

列表 A 、 B 、 C 和 D 展示了这些时间步长方案的 Matlab 实现。在一维中，克莱因-戈登方程具有易于计算的精确解（例如，参见 Nakanishi 和 Schlag^[7]），这些解可用于测试数值方案的准确性。这些方程似乎显示了解决方案行为的三个可能性，这些可能性取决于初始条件

解可能会 *分散* 或 *热化*，也就是说，给定的局部初始条件在整个空间中扩散。
解会爆炸或变成无穷大。
解的一部分以局部粒子的形式传播，而解的其余部分则分散。

由于方程是可逆的，因此也可能存在初始分布在空间域上的解，它会自身局部化。

( A)

一个使用 eq. 中的时间离散化来求解一维克莱因-戈登方程 1 的 Matlab 程序。 2 代码下载

% A program to solve the 1D cubic Klein Gordon equation using a
% second order semi-explicit method
% u_{tt}-u_{xx}+u=u^3
clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
    'defaultaxesfontweight','bold')

% set up grid
tic
Lx = 64; % period 2*pi*L
Nx = 4096; % number of harmonics
Nt = 500; % number of time slices
plotgap=10; % time steps to take between plots
c=0.5; % wave speed
dt = 5.00/Nt; % time step

Es = 1.0; % focusing (+1) or defocusing (-1) parameter
t=0; tdata(1)=t;

% 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
u = sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
uexact= sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
uold=sqrt(2)*sech((x+c*dt)/sqrt(1-c^2));
v=fft(u,[],1);
vold=fft(uold,[],1);
figure(1); clf;
% Plot data on
plot(x,u,'r+',x,uexact,'b-'); legend('numerical','exact');
title(num2str(t)); xlabel x; ylabel u; drawnow;


% initial energy
vx=0.5*kx.*(v+vold);
ux=ifft(vx,[],1);
Kineticenergy=0.5*abs( (u-uold)/dt).^2;
Strainenergy=0.5*abs(ux).^2;
Potentialenergy=0.5*abs(0.5*(u+uold)).^2 ...
                    -Es*0.25*((u+uold)*0.5).^4;
Kineticenergy=fft(Kineticenergy,[],1);
Potentialenergy=fft(Potentialenergy,[],1);
Strainenergy=fft(Strainenergy,[],1);
EnKin(1)=Kineticenergy(1);
EnPot(1)=Potentialenergy(1);
EnStr(1)=Strainenergy(1);
En(1)=EnStr(1)+EnKin(1)+EnPot(1);
En0=En(1)

plotnum=1;
% solve pde and plot results

for n =1:Nt+1
    nonlin=u.^3;
    nonlinhat=fft(nonlin,[],1);
    vnew=(0.25*(kx.*kx -1).*(2*v+vold)...
                +(2*v-vold)/(dt*dt) +Es*nonlinhat)./...
            (1/(dt*dt) - (kx.*kx-1)*0.25 );
    unew=ifft(vnew,[],1);
    t=n*dt;
    if (mod(n,plotgap)==0)
        uexact=sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
        figure(1); clf;
        plot(x,u,'r+',x,uexact,'b-'); legend('numerical','exact');
        title(num2str(t)); xlim([-6,6]); xlabel x; ylabel u; drawnow;
        tdata(plotnum+1)=t;
        vx=0.5*kx.*(v+vold);
        ux=ifft(vx,[],1);
        Kineticenergy=0.5*abs( (u-uold)/dt).^2;
        Strainenergy=0.5*abs(ux).^2;
        Potentialenergy=0.5*abs(0.5*(u+uold)).^2 ...
                    -Es*0.25*((u+uold)*0.5).^4;
        Kineticenergy=fft(Kineticenergy,[],1);
        Potentialenergy=fft(Potentialenergy,[],1);
        Strainenergy=fft(Strainenergy,[],1);
        EnKin(plotnum+1)=Kineticenergy(1);
        EnPot(plotnum+1)=Potentialenergy(1);
        EnStr(plotnum+1)=Strainenergy(1);
        En(plotnum+1)=EnStr(plotnum+1)+EnKin(plotnum+1)+EnPot(plotnum+1);
        Enchange(plotnum)=log(abs(1-En(1+plotnum)/En0));
        plotnum=plotnum+1;
    end
    % update old terms
    vold=v;
    v=vnew;
    uold=u;
    u=unew;
end
figure(4); clf;
uexact=sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
plot(x,u,'r+',x,uexact,'b-'); legend('numerical','exact');
title(num2str(t)); xlabel x; ylabel u; drawnow;
max(abs(u-uexact))
figure(5); clf; plot(tdata,En,'r-',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--');
xlabel time; ylabel Energy; legend('Total','Kinetic','Potential','Strain');
figure(6); clf; plot(tdata(2:end),Enchange,'r-'); xlabel time; ylabel('Energy change');

toc

( Ap1)

一个使用 eq. 中的时间离散化来求解一维克莱因-戈登方程 1 的 Python 程序。 2 。代码下载

"""
A program to solve the 1D Klein Gordon equation using a
second order semi-explicit method. The numerical solution is
compared to an exact solution

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=64.0 # Period 2*pi*Lx
Nx=4096 # Number of harmonics
Nt=500 # Number of time slices
tmax=5.0 # Maximum time
c=0.5	# Wave speed
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)])

kxm=numpy.zeros((Nx), dtype=complex)
xx=numpy.zeros((Nx), dtype=float)

for i in xrange(Nx):
        kxm[i] = k_x[i]
        xx[i] = x[i]
        

# allocate arrays
unew=numpy.zeros((Nx), dtype=float)
u=numpy.zeros((Nx), dtype=float)
uexact=numpy.zeros((Nx), dtype=float)
uold=numpy.zeros((Nx), dtype=float)
vnew=numpy.zeros((Nx), dtype=complex)
v=numpy.zeros((Nx), dtype=complex)
vold=numpy.zeros((Nx), dtype=complex)
ux=numpy.zeros((Nx), dtype=float)
vx=numpy.zeros((Nx), dtype=complex)
Kineticenergy=numpy.zeros((Nx), dtype=complex)
Potentialenergy=numpy.zeros((Nx), dtype=complex)
Strainenergy=numpy.zeros((Nx), dtype=complex)
EnKin=numpy.zeros((numplots), dtype=float)
EnPot=numpy.zeros((numplots), dtype=float)
EnStr=numpy.zeros((numplots), dtype=float)
En=numpy.zeros((numplots), dtype=float)
Enchange=numpy.zeros((numplots-1),dtype=float)
tdata=numpy.zeros((numplots), dtype=float)
nonlin=numpy.zeros((Nx), dtype=float)
nonlinhat=numpy.zeros((Nx), dtype=complex)

t=0.0
u=numpy.sqrt(2)/(numpy.cosh((xx-c*t)/numpy.sqrt(1.0-c**2)))
uexact=numpy.sqrt(2)/(numpy.cosh((xx-c*t)/numpy.sqrt(1.0-c**2)))
uold=numpy.sqrt(2)/(numpy.cosh((xx+c*dt)/numpy.sqrt(1.0-c**2)))
v=numpy.fft.fftn(u)
vold=numpy.fft.fftn(uold)
fig=plt.figure()
ax=fig.add_subplot(211)
ax.plot(xx,u,'b-')
plt.xlabel('x')
plt.ylabel('u')
ax=fig.add_subplot(212)
ax.plot(xx,abs(u-uexact),'b-')
plt.xlabel('x')
plt.ylabel('error')
plt.show()
# initial energy
vx=0.5*kxm*(v+vold)
ux=numpy.real(numpy.fft.ifftn(vx))
Kineticenergy=0.5*((u-uold)/dt)**2
Strainenergy=0.5*(ux)**2
Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
Kineticenergy=numpy.fft.fftn(Kineticenergy)
Strainenergy=numpy.fft.fftn(Strainenergy)
Potentialenergy=numpy.fft.fftn(Potentialenergy)
EnKin[0]=numpy.real(Kineticenergy[0])
EnPot[0]=numpy.real(Potentialenergy[0])
EnStr[0]=numpy.real(Strainenergy[0])
En[0]=EnStr[0]+EnPot[0]+EnKin[0]
EnO=En[0]
tdata[0]=t
plotnum=0
#solve pde and plot results
for nt in xrange(numplots-1):
    for n in xrange(plotgap):
        nonlin=u**3
        nonlinhat=numpy.fft.fftn(nonlin)
        vnew=( (0.25*(kxm**2 - 1)*(2*v+vold)
          +(2*v-vold)/(dt*dt) +Es*nonlinhat)/
          (1/(dt*dt) - (kxm**2 -1)*0.25 ) )
        unew=numpy.real(numpy.fft.ifftn(vnew))
        t+=dt
        # update old terms
        vold=v
        v=vnew
        uold=u
        u=unew
    plotnum+=1
    uexact=numpy.sqrt(2)/(numpy.cosh((xx-c*t)/numpy.sqrt(1.0-c**2)))
    ax = fig.add_subplot(211)
    plt.cla()
    ax.plot(xx,u,'b-')
    plt.title(t)
    plt.xlabel('x')
    plt.ylabel('u')
    ax = fig.add_subplot(212)
    plt.cla()
    ax.plot(xx,abs(u-uexact),'b-')
    plt.xlabel('x')
    plt.ylabel('error')
    plt.draw()
    vx=0.5*kxm*(v+vold)
    ux=numpy.real(numpy.fft.ifftn(vx))
    Kineticenergy=0.5*((u-uold)/dt)**2
    Strainenergy=0.5*(ux)**2
    Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
    Kineticenergy=numpy.fft.fftn(Kineticenergy)
    Strainenergy=numpy.fft.fftn(Strainenergy)
    Potentialenergy=numpy.fft.fftn(Potentialenergy)
    EnKin[plotnum]=numpy.real(Kineticenergy[0])
    EnPot[plotnum]=numpy.real(Potentialenergy[0])
    EnStr[plotnum]=numpy.real(Strainenergy[0])
    En[plotnum]=EnStr[plotnum]+EnPot[plotnum]+EnKin[plotnum]
    Enchange[plotnum-1]=numpy.log(abs(1-En[plotnum]/EnO))
    tdata[plotnum]=t

plt.ioff()

plt.figure()
plt.plot(tdata,En,'r+',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--')
plt.xlabel('Time')
plt.ylabel('Energy')
plt.legend(('Total', 'Kinetic','Potential','Strain'))
plt.title('Time Dependence of Energy Components')
plt.show()

plt.figure()
plt.plot(Enchange,'r-')
plt.title('Time Dependence of Change in Total Energy')
plt.show()

( Ap2)

一个使用 eq. 中的时间离散化来求解一维克莱因-戈登方程 1 的 Python 程序。 3 。代码下载

"""
A program to solve the 1D Klein Gordon equation using a
second order semi-explicit method. The numerical solution is
compared to an exact solution

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=64.0 # Period 2*pi*Lx
Nx=4096 # Number of harmonics
Nt=500 # Number of time slices
tmax=5.0 # Maximum time
c=0.5	# Wave speed
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
tol=0.1**12 # tolerance for fixed point iterations

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

kxm=numpy.zeros((Nx), dtype=complex)
xx=numpy.zeros((Nx), dtype=float)

for i in xrange(Nx):
        kxm[i] = k_x[i]
        xx[i] = x[i]
        
# allocate arrays
unew=numpy.zeros((Nx), dtype=float)
u=numpy.zeros((Nx), dtype=float)
utemp=numpy.zeros((Nx), dtype=float)
uexact=numpy.zeros((Nx), dtype=float)
uold=numpy.zeros((Nx), dtype=float)
vnew=numpy.zeros((Nx), dtype=complex)
v=numpy.zeros((Nx), dtype=complex)
vold=numpy.zeros((Nx), dtype=complex)
ux=numpy.zeros((Nx), dtype=float)
vx=numpy.zeros((Nx), dtype=complex)
Kineticenergy=numpy.zeros((Nx), dtype=complex)
Potentialenergy=numpy.zeros((Nx), dtype=complex)
Strainenergy=numpy.zeros((Nx), dtype=complex)
EnKin=numpy.zeros((numplots), dtype=float)
EnPot=numpy.zeros((numplots), dtype=float)
EnStr=numpy.zeros((numplots), dtype=float)
En=numpy.zeros((numplots), dtype=float)
Enchange=numpy.zeros((numplots-1),dtype=float)
tdata=numpy.zeros((numplots), dtype=float)
nonlin=numpy.zeros((Nx), dtype=float)
nonlinhat=numpy.zeros((Nx), dtype=complex)

t=0.0
u=numpy.sqrt(2)/(numpy.cosh((xx-c*t)/numpy.sqrt(1.0-c**2)))
uexact=numpy.sqrt(2)/(numpy.cosh((xx-c*t)/numpy.sqrt(1.0-c**2)))
uold=numpy.sqrt(2)/(numpy.cosh((xx+c*dt)/numpy.sqrt(1.0-c**2)))
v=numpy.fft.fftn(u)
vold=numpy.fft.fftn(uold)
fig=plt.figure()
ax=fig.add_subplot(211)
ax.plot(xx,u,'b-')
plt.xlabel('x')
plt.ylabel('u')
ax=fig.add_subplot(212)
ax.plot(xx,abs(u-uexact),'b-')
plt.xlabel('x')
plt.ylabel('error')
plt.show()
# initial energy
vx=0.5*kxm*(v+vold)
ux=numpy.real(numpy.fft.ifftn(vx))
Kineticenergy=0.5*((u-uold)/dt)**2
Strainenergy=0.5*(ux)**2
Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
Kineticenergy=numpy.fft.fftn(Kineticenergy)
Strainenergy=numpy.fft.fftn(Strainenergy)
Potentialenergy=numpy.fft.fftn(Potentialenergy)
EnKin[0]=numpy.real(Kineticenergy[0])
EnPot[0]=numpy.real(Potentialenergy[0])
EnStr[0]=numpy.real(Strainenergy[0])
En[0]=EnStr[0]+EnPot[0]+EnKin[0]
EnO=En[0]
tdata[0]=t
plotnum=0
#solve pde and plot results
for nt in xrange(numplots-1):
    for n in xrange(plotgap):
        nonlin=(u**2+uold**2)*(u+uold)/4.0
        nonlinhat=numpy.fft.fftn(nonlin)
        chg=1
        unew=u
        while (chg>tol):
            utemp=unew
            vnew=( (0.25*(kxm**2 - 1)*(2*v+vold)\
                 +(2*v-vold)/(dt*dt) +Es*nonlinhat)\
                 /(1/(dt*dt) - (kxm**2 -1)*0.25 ) )
            unew=numpy.real(numpy.fft.ifftn(vnew))
            nonlin=(unew**2+uold**2)*(unew+uold)/4.0
            nonlinhat=numpy.fft.fftn(nonlin)
            chg=numpy.max(abs(unew-utemp))
        t+=dt
        # update old terms
        vold=v
        v=vnew
        uold=u
        u=unew
    plotnum+=1
    uexact=numpy.sqrt(2)/(numpy.cosh((xx-c*t)/numpy.sqrt(1.0-c**2)))
    ax = fig.add_subplot(211)
    plt.cla()
    ax.plot(xx,u,'b-')
    plt.title(t)
    plt.xlabel('x')
    plt.ylabel('u')
    ax = fig.add_subplot(212)
    plt.cla()
    ax.plot(xx,abs(u-uexact),'b-')
    plt.xlabel('x')
    plt.ylabel('error')
    plt.draw()
    vx=0.5*kxm*(v+vold)
    ux=numpy.real(numpy.fft.ifftn(vx))
    Kineticenergy=0.5*((u-uold)/dt)**2
    Strainenergy=0.5*(ux)**2
    Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
    Kineticenergy=numpy.fft.fftn(Kineticenergy)
    Strainenergy=numpy.fft.fftn(Strainenergy)
    Potentialenergy=numpy.fft.fftn(Potentialenergy)
    EnKin[plotnum]=numpy.real(Kineticenergy[0])
    EnPot[plotnum]=numpy.real(Potentialenergy[0])
    EnStr[plotnum]=numpy.real(Strainenergy[0])
    En[plotnum]=EnStr[plotnum]+EnPot[plotnum]+EnKin[plotnum]
    Enchange[plotnum-1]=numpy.log(abs(1-En[plotnum]/EnO))
    tdata[plotnum]=t

plt.ioff()

plt.figure()
plt.plot(tdata,En,'r+',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--')
plt.xlabel('Time')
plt.ylabel('Energy')
plt.legend(('Total', 'Kinetic','Potential','Strain'))
plt.title('Time Dependence of Energy Components')
plt.show()

plt.figure()
plt.plot(Enchange,'r-')
plt.title('Time Dependence of Change in Total Energy')
plt.show()

( C)

一个使用 eq. 中的时间离散化来求解二维克莱因-戈登方程 1 的 Matlab 程序。 3 代码下载

% A program to solve the 1D cubic Klein Gordon equation using a
% second order implicit method
% u_{tt}-u_{xx}+u=u^3
clear all; format compact; format short;
set(0,'defaultaxesfontsize',30,'defaultaxeslinewidth',.7,...
    'defaultlinelinewidth',6,'defaultpatchlinewidth',3.7,...
    'defaultaxesfontweight','bold')

% set up grid
tic
Lx = 64; % period 2*pi*L
Nx = 4096; % number of harmonics
Nt = 400; % number of time slices
plotgap=10; % timesteps between plots
tol=0.1^(15); % tolerance for fixed point iterations
dt = 0.500/Nt; % time step
c=0.5; % wave speed

Es = 1.0; % focusing (+1) or defocusing (-1) parameter
t=0; tdata(1)=t;

% 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
u = sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
uexact= sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
uold=sqrt(2)*sech((x+c*dt)/sqrt(1-c^2));
v=fft(u,[],1);
vold=fft(uold,[],1);
figure(1); clf;
% Plot data on
plot(x,u,'r+',x,uexact,'b-'); legend('numerical','exact');
title(num2str(0)); xlim([-6,6]); xlabel x; ylabel u; drawnow;


% initial energy
vx=0.5*kx.*(v+vold);
ux=ifft(vx,[],1);
Kineticenergy=0.5*abs( (u-uold)/dt).^2;
Strainenergy=0.5*abs(ux).^2;
Potentialenergy=0.5*abs(0.5*(u+uold)).^2 ...
                    -Es*0.25*((u+uold)*0.5).^4;
Kineticenergy=fft(Kineticenergy,[],1);
Potentialenergy=fft(Potentialenergy,[],1);
Strainenergy=fft(Strainenergy,[],1);
EnKin(1)=Kineticenergy(1);
EnPot(1)=Potentialenergy(1);
EnStr(1)=Strainenergy(1);
En(1)=EnStr(1)+EnKin(1)+EnPot(1);
En0=En(1)

plotnum=1;
% solve pde and plot results

for n =1:Nt+1
   nonlin=(u.^2 +uold.^2).*(u+uold)/4;
    nonlinhat=fft(nonlin,[],1);
    chg=1;
    unew=u;
    while (chg>tol)
        utemp=unew;
        vnew=(0.25*(kx.*kx -1).*(2*v+vold)...
                +(2*v-vold)/(dt*dt) +Es*nonlinhat)./...
            (1/(dt*dt) - (kx.*kx -1)*0.25 );
        unew=ifft(vnew,[],1);
        nonlin=(unew.^2 +uold.^2).*(unew+uold)/4;
        nonlinhat=fft(nonlin,[],1);
        chg=max(abs(unew-utemp));
    end
    t=n*dt;
    if (mod(n,plotgap)==0)
        uexact=sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
        figure(1); clf;
        plot(x,u,'r+',x,uexact,'b-'); legend('numerical','exact');
        title(num2str(t)); xlim([-6,6]); xlabel x; ylabel u; drawnow;
        tdata(plotnum+1)=t;
        vx=0.5*kx.*(v+vold);
        ux=ifft(vx,[],1);
        Kineticenergy=0.5*abs( (u-uold)/dt).^2;
        Strainenergy=0.5*abs(ux).^2;
        Potentialenergy=0.5*abs(0.5*(u+uold)).^2 ...
                    -Es*0.25*((u+uold)*0.5).^4;
        Kineticenergy=fft(Kineticenergy,[],1);
        Potentialenergy=fft(Potentialenergy,[],1);
        Strainenergy=fft(Strainenergy,[],1);
        EnKin(plotnum+1)=Kineticenergy(1);
        EnPot(plotnum+1)=Potentialenergy(1);
        EnStr(plotnum+1)=Strainenergy(1);
        En(plotnum+1)=EnStr(plotnum+1)+EnKin(plotnum+1)+EnPot(plotnum+1);
        Enchange(plotnum)=log(abs(1-En(1+plotnum)/En0));
        plotnum=plotnum+1;
    end
    % update old terms
    vold=v;
    v=vnew;
    uold=u;
    u=unew;
end
figure(4); clf;
uexact=sqrt(2)*sech((x-c*t)/sqrt(1-c^2));
plot(x,u,'r+',x,uexact,'b-'); legend('numerical','exact');
title(num2str(t)); xlim([-6,6]); xlabel x; ylabel u; drawnow;
max(abs(u-uexact))
figure(5); clf; plot(tdata,En,'r-',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--');
xlabel time; ylabel Energy; legend('Total','Kinetic','Potential','Strain');
figure(6); clf; plot(tdata(2:end),Enchange,'r-'); xlabel time; ylabel('Energy change');

toc

( Cp)

一个使用 eq. 中的时间离散化来求解二维克莱因-戈登方程 1 的 Python 程序。 3 。代码下载

#!/usr/bin/env python
"""
A program to solve the 2D Klein Gordon equation using a
second order semi-explicit 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=3.0 # Period 2*pi*Lx
Ly=3.0 # Period 2*pi*Ly
Nx=512 # Number of harmonics
Ny=512	# Number of harmonics
Nt=200 # Number of time slices
tmax=5.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)])

kxm=numpy.zeros((Nx,Ny), dtype=complex)
kym=numpy.zeros((Nx,Ny), dtype=complex)
xx=numpy.zeros((Nx,Ny), dtype=float)
yy=numpy.zeros((Nx,Ny), dtype=float)


for i in xrange(Nx):
    for j in xrange(Ny):
        kxm[i,j] = k_x[i]
        kym[i,j] = k_y[j]
        xx[i,j] = x[i]
        yy[i,j] = y[j]
        

# allocate arrays
unew=numpy.zeros((Nx,Ny), dtype=float)
u=numpy.zeros((Nx,Ny), dtype=float)
uold=numpy.zeros((Nx,Ny), dtype=float)
vnew=numpy.zeros((Nx,Ny), dtype=complex)
v=numpy.zeros((Nx,Ny), dtype=complex)
vold=numpy.zeros((Nx,Ny), dtype=complex)
ux=numpy.zeros((Nx,Ny), dtype=float)
uy=numpy.zeros((Nx,Ny), dtype=float)
vx=numpy.zeros((Nx,Ny), dtype=complex)
vy=numpy.zeros((Nx,Ny), dtype=complex)
Kineticenergy=numpy.zeros((Nx,Ny), dtype=complex)
Potentialenergy=numpy.zeros((Nx,Ny), dtype=complex)
Strainenergy=numpy.zeros((Nx,Ny), dtype=complex)
EnKin=numpy.zeros((numplots), dtype=float)
EnPot=numpy.zeros((numplots), dtype=float)
EnStr=numpy.zeros((numplots), dtype=float)
En=numpy.zeros((numplots), dtype=float)
Enchange=numpy.zeros((numplots-1),dtype=float)
tdata=numpy.zeros((numplots), dtype=float)
nonlin=numpy.zeros((Nx,Ny), dtype=float)
nonlinhat=numpy.zeros((Nx,Ny), dtype=complex)

u=0.1*numpy.exp(-(xx**2 + yy**2))*numpy.sin(10*xx+12*yy)
uold=u
v=numpy.fft.fft2(u)
vold=numpy.fft.fft2(uold)
src = mlab.surf(xx,yy,u,colormap='YlGnBu',warp_scale='auto')
mlab.scalarbar(object=src)
mlab.xlabel('x',object=src)
mlab.ylabel('y',object=src)
mlab.zlabel('u',object=src)
# initial energy
vx=0.5*kxm*(v+vold)
vy=0.5*kym*(v+vold)
ux=numpy.fft.ifft2(vx)
uy=numpy.fft.ifft2(vy)
Kineticenergy=0.5*((u-uold)/dt)**2
Strainenergy=0.5*(ux)**2 + 0.5*(uy)**2
Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
Kineticenergy=numpy.fft.fft2(Kineticenergy)
Strainenergy=numpy.fft.fft2(Strainenergy)
Potentialenergy=numpy.fft.fft2(Potentialenergy)
EnKin[0]=numpy.real(Kineticenergy[0,0])
EnPot[0]=numpy.real(Potentialenergy[0,0])
EnStr[0]=numpy.real(Strainenergy[0,0])
En[0]=EnStr[0]+EnPot[0]+EnKin[0]
EnO=En[0]
t=0.0
tdata[0]=t
plotnum=0
#solve pde and plot results
for nt in xrange(numplots-1):
    for n in xrange(plotgap):
        nonlin=u**3
        nonlinhat=numpy.fft.fft2(nonlin)
        vnew=( (0.25*(kxm**2 + kym**2 - 1)*(2*v+vold)
          +(2*v-vold)/(dt*dt) +Es*nonlinhat)/
          (1/(dt*dt) - (kxm**2 + kym**2 -1)*0.25 ) )
        unew=numpy.real(numpy.fft.ifft2(vnew))
        t+=dt
        # update old terms
        vold=v
        v=vnew
        uold=u
        u=unew
    plotnum+=1
    src.mlab_source.scalars = unew
    vx=0.5*kxm*(v+vold)
    vy=0.5*kym*(v+vold)
    ux=numpy.fft.ifft2(vx)
    uy=numpy.fft.ifft2(vy)
    Kineticenergy=0.5*((u-uold)/dt)**2
    Strainenergy=0.5*(ux)**2 + 0.5*(uy)**2
    Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
    Kineticenergy=numpy.fft.fft2(Kineticenergy)
    Strainenergy=numpy.fft.fft2(Strainenergy)
    Potentialenergy=numpy.fft.fft2(Potentialenergy)
    EnKin[plotnum]=numpy.real(Kineticenergy[0,0])
    EnPot[plotnum]=numpy.real(Potentialenergy[0,0])
    EnStr[plotnum]=numpy.real(Strainenergy[0,0])
    En[plotnum]=EnStr[plotnum]+EnPot[plotnum]+EnKin[plotnum]
    Enchange[plotnum-1]=numpy.log(abs(1-En[plotnum]/EnO))
    tdata[plotnum]=t

       
plt.figure()
plt.plot(tdata,En,'r+',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--')
plt.xlabel('Time')
plt.ylabel('Energy')
plt.legend(('Total', 'Kinetic','Potential','Strain'))
plt.title('Time Dependence of Energy Components')
plt.show()

plt.figure()
plt.plot(Enchange,'r-')
plt.title('Time Dependence of Change in Total Energy')
plt.show()

( Cp2)

使用等式 2 中的时间离散化，编写一个 Python 程序来求解二维克莱因-戈登方程 1 。代码下载

"""
A program to solve the 2D Klein Gordon equation using a
second order semi-explicit 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=3.0 # Period 2*pi*Lx
Ly=3.0 # Period 2*pi*Ly
Nx=512 # Number of harmonics
Ny=512	# Number of harmonics
Nt=200 # Number of time slices
tmax=5.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)])

kxm=numpy.zeros((Nx,Ny), dtype=complex)
kym=numpy.zeros((Nx,Ny), dtype=complex)
xx=numpy.zeros((Nx,Ny), dtype=float)
yy=numpy.zeros((Nx,Ny), dtype=float)


for i in xrange(Nx):
    for j in xrange(Ny):
        kxm[i,j] = k_x[i]
        kym[i,j] = k_y[j]
        xx[i,j] = x[i]
        yy[i,j] = y[j]
        

# allocate arrays
unew=numpy.zeros((Nx,Ny), dtype=float)
u=numpy.zeros((Nx,Ny), dtype=float)
uold=numpy.zeros((Nx,Ny), dtype=float)
vnew=numpy.zeros((Nx,Ny), dtype=complex)
v=numpy.zeros((Nx,Ny), dtype=complex)
vold=numpy.zeros((Nx,Ny), dtype=complex)
ux=numpy.zeros((Nx,Ny), dtype=float)
uy=numpy.zeros((Nx,Ny), dtype=float)
vx=numpy.zeros((Nx,Ny), dtype=complex)
vy=numpy.zeros((Nx,Ny), dtype=complex)
Kineticenergy=numpy.zeros((Nx,Ny), dtype=complex)
Potentialenergy=numpy.zeros((Nx,Ny), dtype=complex)
Strainenergy=numpy.zeros((Nx,Ny), dtype=complex)
EnKin=numpy.zeros((numplots), dtype=float)
EnPot=numpy.zeros((numplots), dtype=float)
EnStr=numpy.zeros((numplots), dtype=float)
En=numpy.zeros((numplots), dtype=float)
Enchange=numpy.zeros((numplots-1),dtype=float)
tdata=numpy.zeros((numplots), dtype=float)
nonlin=numpy.zeros((Nx,Ny), dtype=float)
nonlinhat=numpy.zeros((Nx,Ny), dtype=complex)

u=0.1*numpy.exp(-(xx**2 + yy**2))*numpy.sin(10*xx+12*yy)
uold=u
v=numpy.fft.fft2(u)
vold=numpy.fft.fft2(uold)
src = mlab.surf(xx,yy,u,colormap='YlGnBu',warp_scale='auto')
mlab.scalarbar(object=src)
mlab.xlabel('x',object=src)
mlab.ylabel('y',object=src)
mlab.zlabel('u',object=src)
# initial energy
vx=0.5*kxm*(v+vold)
vy=0.5*kym*(v+vold)
ux=numpy.fft.ifft2(vx)
uy=numpy.fft.ifft2(vy)
Kineticenergy=0.5*((u-uold)/dt)**2
Strainenergy=0.5*(ux)**2 + 0.5*(uy)**2
Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
Kineticenergy=numpy.fft.fft2(Kineticenergy)
Strainenergy=numpy.fft.fft2(Strainenergy)
Potentialenergy=numpy.fft.fft2(Potentialenergy)
EnKin[0]=numpy.real(Kineticenergy[0,0])
EnPot[0]=numpy.real(Potentialenergy[0,0])
EnStr[0]=numpy.real(Strainenergy[0,0])
En[0]=EnStr[0]+EnPot[0]+EnKin[0]
EnO=En[0]
t=0.0
tdata[0]=t
plotnum=0
#solve pde and plot results
for nt in xrange(numplots-1):
    for n in xrange(plotgap):
        nonlin=u**3
        nonlinhat=numpy.fft.fft2(nonlin)
        vnew=( (0.25*(kxm**2 + kym**2 - 1)*(2*v+vold)
          +(2*v-vold)/(dt*dt) +Es*nonlinhat)/
          (1/(dt*dt) - (kxm**2 + kym**2 -1)*0.25 ) )
        unew=numpy.real(numpy.fft.ifft2(vnew))
        t+=dt
        # update old terms
        vold=v
        v=vnew
        uold=u
        u=unew
    plotnum+=1
    src.mlab_source.scalars = unew
    vx=0.5*kxm*(v+vold)
    vy=0.5*kym*(v+vold)
    ux=numpy.fft.ifft2(vx)
    uy=numpy.fft.ifft2(vy)
    Kineticenergy=0.5*((u-uold)/dt)**2
    Strainenergy=0.5*(ux)**2 + 0.5*(uy)**2
    Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
    Kineticenergy=numpy.fft.fft2(Kineticenergy)
    Strainenergy=numpy.fft.fft2(Strainenergy)
    Potentialenergy=numpy.fft.fft2(Potentialenergy)
    EnKin[plotnum]=numpy.real(Kineticenergy[0,0])
    EnPot[plotnum]=numpy.real(Potentialenergy[0,0])
    EnStr[plotnum]=numpy.real(Strainenergy[0,0])
    En[plotnum]=EnStr[plotnum]+EnPot[plotnum]+EnKin[plotnum]
    Enchange[plotnum-1]=numpy.log(abs(1-En[plotnum]/EnO))
    tdata[plotnum]=t

       
plt.figure()
plt.plot(tdata,En,'r+',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--')
plt.xlabel('Time')
plt.ylabel('Energy')
plt.legend(('Total', 'Kinetic','Potential','Strain'))
plt.title('Time Dependence of Energy Components')
plt.show()

plt.figure()
plt.plot(Enchange,'r-')
plt.title('Time Dependence of Change in Total Energy')
plt.show()

( D)

使用等式 2 中的时间离散化，编写一个 Matlab 程序来求解三维克莱因-戈登方程 1 。代码下载

% A program to solve the 3D Klein Gordon equation using a
% second order semi-explicit 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 = 2; % period 2*pi*L
Ly = 2; % period 2*pi*L
Lz = 2; % period 2*pi*L
Nx = 64; % number of harmonics
Ny = 64; % number of harmonics
Nz = 64; % number of harmonics
Nt = 2000; % number of time slices
plotgap=10;
dt = 10.0/Nt; % time step

Es = -1.0; % focusing (+1) or defocusing (-1) 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);
[kxm,kym,kzm]=meshgrid(kx,ky,kz);

% initial conditions
u = 0.1*exp(-(xx.^2+(yy).^2+zz.^2));
uold=u;
v=fftn(u);
vold=v;
figure(1); clf;
% coordinate slice to show plots on
sx=[0]; sy=[0]; sz=[-Lx*2*pi];
slice(xx,yy,zz,u,sx,sy,sz); colormap jet;
title(num2str(0)); colorbar('location','EastOutside'); drawnow;

xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); drawnow;
t=0; tdata(1)=t;

% initial energy
vx=0.5*kxm.*(v+vold);
vy=0.5*kym.*(v+vold);
vz=0.5*kzm.*(v+vold);
ux=ifftn(vx);
uy=ifftn(vy);
uz=ifftn(vz);
Kineticenergy=0.5*abs( (u-uold)/dt).^2;
Strainenergy=0.5*abs(ux).^2 +0.5*abs(uy).^2+0.5*abs(uz).^2;
Potentialenergy=0.5*abs(0.5*(u+uold)).^2 ...
                    -Es*0.25*((u+uold)*0.5).^4;
Kineticenergy=fftn(Kineticenergy);
Potentialenergy=fftn(Potentialenergy);
Strainenergy=fftn(Strainenergy);
EnKin(1)=Kineticenergy(1,1);
EnPot(1)=Potentialenergy(1,1);
EnStr(1)=Strainenergy(1,1);
En(1)=EnStr(1)+EnKin(1)+EnPot(1);
En0=En(1)

plotnum=1;
% solve pde and plot results

for n =1:Nt+1
    nonlin=u.^3;
    nonlinhat=fftn(nonlin);
    vnew=(0.25*(kxm.^2 + kym.^2 + kzm.^2 -1).*(2*v+vold)...
                +(2*v-vold)/(dt*dt) +Es*nonlinhat)./...
            (1/(dt*dt) - (kxm.^2 + kzm.^2 + kym.^2 - 1)*0.25 );
    unew=ifftn(vnew);
    t=n*dt;
    if (mod(n,plotgap)==0)
        figure(1); clf; sx=[0]; sy=[0]; sz=[0];
        slice(xx,yy,zz,u,sx,sy,sz); colormap jet;
        title(num2str(t)); colorbar('location','EastOutside'); drawnow;
        xlabel('x'); ylabel('y'); zlabel('z');
        axis equal; axis square; view(3); drawnow;
        tdata(plotnum+1)=t;
        t
        vx=0.5*kxm.*(v+vold);
        vy=0.5*kym.*(v+vold);
        vz=0.5*kzm.*(v+vold);
        ux=ifftn(vx);
        uy=ifftn(vy);
        uz=ifftn(vz);
        Kineticenergy=0.5*abs( (u-uold)/dt).^2;
        Strainenergy=0.5*abs(ux).^2 +0.5*abs(uy).^2+0.5*abs(uz).^2;
        Potentialenergy=0.5*abs(0.5*(u+uold)).^2 ...
                    -Es*0.25*((u+uold)*0.5).^4;
        Kineticenergy=fftn(Kineticenergy);
        Potentialenergy=fftn(Potentialenergy);
        Strainenergy=fftn(Strainenergy);
        EnKin(plotnum+1)=Kineticenergy(1,1,1);
        EnPot(plotnum+1)=Potentialenergy(1,1,1);
        EnStr(plotnum+1)=Strainenergy(1,1,1);
        En(plotnum+1)=EnStr(plotnum+1)+EnKin(plotnum+1)+EnPot(plotnum+1);
        Enchange(plotnum)=log(abs(1-En(1+plotnum)/En0));
        plotnum=plotnum+1;
    end
    % update old terms
    vold=v;
    v=vnew;
    uold=u;
    u=unew;
end
figure(4); clf;
% coordinate slice to show plots on
sx=[0]; sy=[0]; sz=[0];
slice(xx,yy,zz,u,sx,sy,sz); colormap jet;
title(num2str(t)); colorbar('location','EastOutside'); drawnow;

xlabel('x'); ylabel('y'); zlabel('z');
axis equal; axis square; view(3); drawnow;

figure(5); clf; plot(tdata,En,'r-',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--');
xlabel time; ylabel Energy; legend('Total','Kinetic','Potential','Strain');
figure(6); clf; plot(tdata(2:end),Enchange,'r-'); xlabel time; ylabel('Energy change');

toc

( Dp)

使用等式 2 中的时间离散化，编写一个 Python 程序来求解三维克莱因-戈登方程，即等式 1 。代码下载

#!/usr/bin/env python
"""
A program to solve the 3D Klein Gordon equation using a
second order semi-explicit 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=2.0 # Period 2*pi*Lx
Ly=2.0 # Period 2*pi*Ly
Lz=2.0 # Period 2*pi*Lz
Nx=64 # Number of harmonics
Ny=64 # Number of harmonics
Nz=64 # Number of harmonics
Nt=2000 # Number of time slices
tmax=10.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)]
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)
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]
            xx[i,j,k]=x[i]
            yy[i,j,k]=y[j]
            zz[i,j,k]=z[k]
        

# allocate arrays
unew=numpy.zeros((Nx,Ny,Nz), dtype=float)
u=numpy.zeros((Nx,Ny,Nz), dtype=float)
uold=numpy.zeros((Nx,Ny,Nz), dtype=float)
vnew=numpy.zeros((Nx,Ny,Nz), dtype=complex)
v=numpy.zeros((Nx,Ny,Nz), dtype=complex)
vold=numpy.zeros((Nx,Ny,Nz), dtype=complex)
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=complex)
vy=numpy.zeros((Nx,Ny,Nz), dtype=complex)
vz=numpy.zeros((Nx,Ny,Nz), dtype=complex)
Kineticenergy=numpy.zeros((Nx,Ny,Nz), dtype=complex)
Potentialenergy=numpy.zeros((Nx,Ny,Nz), dtype=complex)
Strainenergy=numpy.zeros((Nx,Ny,Nz), dtype=complex)
EnKin=numpy.zeros((numplots), dtype=float)
EnPot=numpy.zeros((numplots), dtype=float)
EnStr=numpy.zeros((numplots), dtype=float)
En=numpy.zeros((numplots), dtype=float)
Enchange=numpy.zeros((numplots-1),dtype=float)
tdata=numpy.zeros((numplots), dtype=float)
nonlin=numpy.zeros((Nx,Ny,Nz), dtype=float)
nonlinhat=numpy.zeros((Nx,Ny,Nz), dtype=complex)

u=0.1*numpy.exp(-(xx**2 + yy**2 + zz**2))
uold=u
v=numpy.fft.fftn(u)
vold=numpy.fft.fftn(uold)
#src=mlab.contour3d(xx,yy,zz,u,colormap='jet',opacity=0.1,contours=4)
src = mlab.pipeline.scalar_field(xx,yy,zz,u,colormap='YlGnBu')
mlab.pipeline.iso_surface(src, contours=[u.min()+0.1*u.ptp(), ],
   colormap='YlGnBu',opacity=0.85)
mlab.pipeline.iso_surface(src, contours=[u.max()-0.1*u.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)

# initial energy
vx=0.5*kxm*(v+vold)
vy=0.5*kym*(v+vold)
vz=0.5*kzm*(v+vold)
ux=numpy.fft.ifftn(vx)
uy=numpy.fft.ifftn(vy)
uz=numpy.fft.ifftn(vz)
Kineticenergy=0.5*((u-uold)/dt)**2
Strainenergy=0.5*(ux)**2 + 0.5*(uy)**2 + 0.5*(uz)**2
Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
Kineticenergy=numpy.fft.fftn(Kineticenergy)
Strainenergy=numpy.fft.fftn(Strainenergy)
Potentialenergy=numpy.fft.fftn(Potentialenergy)
EnKin[0]=numpy.real(Kineticenergy[0,0,0])
EnPot[0]=numpy.real(Potentialenergy[0,0,0])
EnStr[0]=numpy.real(Strainenergy[0,0,0])
En[0]=EnStr[0]+EnPot[0]+EnKin[0]
EnO=En[0]
t=0.0
tdata[1]=t
plotnum=0
#solve pde and plot results
for nt in xrange(numplots-1):
    for n in xrange(plotgap):
        nonlin=u**3
        nonlinhat=numpy.fft.fftn(nonlin)
        vnew=( (0.25*(kxm**2 + kym**2 + kzm**2 - 1)*(2*v+vold)
          +(2*v-vold)/(dt*dt) +Es*nonlinhat)/
          (1/(dt*dt) - (kxm**2 + kym**2 + kzm**2 -1)*0.25 ) )
        unew=numpy.real(numpy.fft.ifftn(vnew))
        t+=dt
        # update old terms
        vold=v
        v=vnew
        uold=u
        u=unew
    plotnum+=1
    src.mlab_source.scalars = unew
    vx=0.5*kxm*(v+vold)
    vy=0.5*kym*(v+vold)
    vz=0.5*kzm*(v+vold)
    ux=numpy.fft.ifftn(vx)
    uy=numpy.fft.ifftn(vy)
    uz=numpy.fft.ifftn(vz)
    Kineticenergy=0.5*((u-uold)/dt)**2
    Strainenergy=0.5*(ux)**2 + 0.5*(uy)**2 + 0.5*(uz)**2
    Potentialenergy=0.5*(0.5*(u+uold))**2 - Es*0.25*(0.5*(u+uold))**4
    Kineticenergy=numpy.fft.fftn(Kineticenergy)
    Strainenergy=numpy.fft.fftn(Strainenergy)
    Potentialenergy=numpy.fft.fftn(Potentialenergy)
    EnKin[plotnum]=numpy.real(Kineticenergy[0,0,0])
    EnPot[plotnum]=numpy.real(Potentialenergy[0,0,0])
    EnStr[plotnum]=numpy.real(Strainenergy[0,0,0])
    En[plotnum]=EnStr[plotnum]+EnPot[plotnum]+EnKin[plotnum]
    Enchange[plotnum-1]=numpy.log(abs(1-En[plotnum]/EnO))
    tdata[plotnum]=t

       
plt.figure()
plt.plot(tdata,En,'r+',tdata,EnKin,'b:',tdata,EnPot,'g-.',tdata,EnStr,'y--')
plt.xlabel('Time')
plt.ylabel('Energy')
plt.legend(('Total', 'Kinetic','Potential','Strain'))
plt.title('Time Dependence of Energy Components')
plt.show()

plt.figure()
plt.plot(Enchange,'r-')
plt.title('Time Dependence of Change in Total Energy')
plt.show()

二维 OpenMP Fortran 程序

我们用 Fortran 开发的程序已经相当长了。在这里，我们添加子程序以使程序更短，更易于维护。清单 E 是使用 OpenMP 求解二维克莱因-戈登方程的主 Fortran 程序。注意，通过使用子程序，我们已经使主程序明显缩短，更易于阅读。它仍然不如 Matlab 程序那样简单，但比以前的一些 Fortran 程序好得多。它也更容易维护，并且一旦子程序被编写和调试，就可以在其他程序中重复使用。使用过多子程序的唯一缺点是，由于调用子程序和向其传递数据的开销，可能会导致性能略微下降。子程序在清单 F 、 G 、 H 、 I 、 J 、 K 中，并且示例 makefile 在清单 L 中。最后，清单 M 包含一个 Matlab 程序，该程序从已计算的二进制文件中生成图片。然后可以使用另一个程序来获取图像并创建视频^[8]。

( E)

一个 Fortran 程序来求解二维克莱因-戈登方程代码下载

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program solves nonlinear Klein-Gordon equation in 2 dimensions
! u_{tt}-u_{xx}+u_{yy}+u=Es*|u|^2u
! using a second order implicit-explicit time stepping 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=0.5*exp(-x^2-y^2)*sin(10*x+12*y)
!
! .. Parameters ..
! Nx = number of modes in x - power of 2 for FFT
! Ny = number of modes in y - power of 2 for FFT
! Nt = number of timesteps to take
! Tmax = maximum simulation time
! 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
! planfxy = Forward 2d fft plan
! planbxy = Backward 2d fft plan
! dt = timestep
! ierr = error code
! plotnum = number of plot
! .. Arrays ..
! unew = approximate solution
! vnew = Fourier transform of approximate solution
! u = approximate solution
! v = Fourier transform of approximate solution
! uold = approximate solution
! vold = Fourier transform of approximate solution
! nonlin = nonlinear term, u^3
! nonlinhat = Fourier transform of nonlinear term, u^3
! .. 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
! en = total energy
! enstr = strain energy
! enpot = potential energy
! enkin = kinetic energy
! 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
! getgrid.f90 -- Get initial grid of points
! initialdata.f90 -- Get initial data
! enercalc.f90 -- Subroutine to calculate the energy
! savedata.f90 -- Save initial data
! storeold.f90 -- Store old data
! External libraries required
! FFTW3 -- Fast Fourier Transform in the West Library
! (http://www.fftw.org/)
! OpenMP library

PROGRAM Kg
USE omp_lib	
! Declare variables
IMPLICIT NONE
INTEGER(kind=4), PARAMETER	:: Nx=128
INTEGER(kind=4), PARAMETER :: Ny=128
INTEGER(kind=4), PARAMETER	:: Nt=20
INTEGER(kind=4), PARAMETER	:: plotgap=5	
REAL(kind=8), PARAMETER	:: &
pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8), PARAMETER	:: Lx=3.0d0
REAL(kind=8), PARAMETER	:: Ly=3.0d0
REAL(kind=8), PARAMETER	:: Es=1.0d0	
REAL(kind=8)	:: dt=0.10d0/REAL(Nt,kind(0d0))	
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx,ky
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x,y	
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: u,nonlin
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: v,nonlinhat
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: uold
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: vold
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: unew
COMPLEX(kind=8), DIMENSION(:,:), ALLOCATABLE:: vnew
REAL(kind=8), DIMENSION(:,:), ALLOCATABLE :: savearray
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: time,enkin,enstr,enpot,en
INTEGER(kind=4)	:: ierr,i,j,n,allocatestatus,numthreads
INTEGER(kind=4)	:: start, finish, count_rate, plotnum
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)	:: planfxy,planbxy
CHARACTER*100	:: name_config
! Start short parallel region to count threads
numthreads=omp_get_max_threads()
PRINT *,'There are ',numthreads,' threads.'
ALLOCATE(kx(1:Nx),ky(1:Ny),x(1:Nx),y(1:Ny),u(1:Nx,1:Ny),&
v(1:Nx,1:Ny),nonlin(1:Nx,1:Ny),nonlinhat(1:Nx,1:Ny),&
uold(1:Nx,1:Ny),vold(1:Nx,1:Ny),&
unew(1:Nx,1:Ny),vnew(1:Nx,1:Ny),savearray(1:Nx,1:Ny),&
time(1:1+Nt/plotgap),enkin(1:1+Nt/plotgap),&
enstr(1:1+Nt/plotgap),enpot(1:1+Nt/plotgap),&
en(1:1+Nt/plotgap),stat=allocatestatus)	
IF (allocatestatus .ne. 0) stop
PRINT *,'allocated arrays'

! 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,u,v,&
FFTW_FORWARD,FFTW_ESTIMATE)
  CALL dfftw_plan_dft_2d_(planbxy,Nx,Ny,v,u,&
  FFTW_BACKWARD,FFTW_ESTIMATE)
PRINT *,'Setup FFTs'
! setup fourier frequencies
CALL getgrid(Nx,Ny,Lx,Ly,pi,name_config,x,y,kx,ky)
PRINT *,'Setup grid and fourier frequencies'
CALL initialdata(Nx,Ny,x,y,u,uold)
plotnum=1	
name_config = 'data/u'
savearray=REAL(u)
CALL savedata(Nx,Ny,plotnum,name_config,savearray)

CALL dfftw_execute_dft_(planfxy,u,v)
CALL dfftw_execute_dft_(planfxy,uold,vold)

CALL enercalc(Nx,Ny,planfxy,planbxy,dt,Es,&
enkin(plotnum),enstr(plotnum),&
enpot(plotnum),en(plotnum),&
kx,ky,nonlin,nonlinhat,&
v,vold,u,uold)

PRINT *,'Got initial data, starting timestepping'
time(plotnum)=0.0d0
   CALL system_clock(start,count_rate)	
DO n=1,Nt	
!$OMP PARALLEL DO PRIVATE(i,j) SCHEDULE(static)
DO j=1,Ny
DO i=1,Nx
nonlin(i,j)=(abs(u(i,j))*2)*u(i,j)
END DO
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planfxy,nonlin,nonlinhat)
!$OMP PARALLEL DO PRIVATE(i,j) SCHEDULE(static)
DO j=1,Ny
DO i=1,Nx
vnew(i,j)=( 0.25*(kx(i)*kx(i) + ky(j)*ky(j)-1.0d0)&
*(2.0d0*v(i,j)+vold(i,j))&
+(2.0d0*v(i,j)-vold(i,j))/(dt*dt)&
+Es*nonlinhat(i,j) )&
/(1/(dt*dt)-0.25*(kx(i)*kx(i) + ky(j)*ky(j)-1.0d0))
END DO
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planbxy,vnew,unew)
! normalize result
!$OMP PARALLEL DO PRIVATE(i,j) SCHEDULE(static)
DO j=1,Ny
DO i=1,Nx
unew(i,j)=unew(i,j)/REAL(Nx*Ny,kind(0d0))
END DO
END DO
!$OMP END PARALLEL DO
IF (mod(n,plotgap)==0) then
plotnum=plotnum+1
time(plotnum)=n*dt
PRINT *,'time',n*dt
CALL enercalc(Nx,Ny,planfxy,planbxy,dt,Es,&
enkin(plotnum),enstr(plotnum),&
enpot(plotnum),en(plotnum),&
kx,ky,&
nonlin,nonlinhat,&
vnew,v,unew,u)
savearray=REAL(unew,kind(0d0))
CALL savedata(Nx,Ny,plotnum,name_config,savearray)
END IF
! .. Update old values ..
CALL storeold(Nx,Ny,&
unew,u,uold,&
vnew,v,vold)
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 Time stepping'
CALL saveresults(Nt,plotgap,time(1:1+n/plotgap),en(1:1+n/plotgap),&
enstr(1:1+n/plotgap),enkin(1:1+n/plotgap),enpot(1:1+n/plotgap))	

! Save times at which output was made in text format
PRINT *,'Saved data'

CALL dfftw_destroy_plan_(planbxy)
CALL dfftw_destroy_plan_(planfxy)
CALL dfftw_cleanup_threads_()
DEALLOCATE(kx,ky,x,y,u,v,nonlin,nonlinhat,savearray,&
uold,vold,unew,vnew,time,enkin,enstr,enpot,en,&
stat=allocatestatus)	
IF (allocatestatus .ne. 0) STOP
PRINT *,'Deallocated arrays'
PRINT *,'Program execution complete'
END PROGRAM Kg

( F)

一个 Fortran 子程序来获取在二维克莱因-戈登方程上求解的网格代码下载

SUBROUTINE getgrid(Nx,Ny,Lx,Ly,pi,name_config,x,y,kx,ky)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine gets grid points and fourier frequencies for a
! pseudospectral simulation of the 2D nonlinear Klein-Gordon equation
!
! u_{tt}-u_{xx}+u_{yy}+u=Es*u^3
!
! The boundary conditions are u(x=0,y)=u(2*Lx*\pi,y),
! u(x,y=0)=u(x,y=2*Ly*\pi)
!
! INPUT
!
! .. Scalars ..
! Nx = number of modes in x - power of 2 for FFT
! Ny = number of modes in y - power of 2 for FFT
! pi = 3.142....
! Lx = width of box in x direction
! Ly = width of box in y direction
! OUTPUT
!
! .. Vectors ..
! kx = fourier frequencies in x direction
! ky = fourier frequencies in y direction
! x = x locations
! y = y locations
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
!
! 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
! OpenMP library
IMPLICIT NONE
USE omp_lib	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny
REAL(kind=8), INTENT(IN)	:: Lx,Ly,pi
REAL(KIND=8), DIMENSION(1:NX), INTENT(OUT) :: x
REAL(KIND=8), DIMENSION(1:NY), INTENT(OUT) :: y
COMPLEX(KIND=8), DIMENSION(1:NX), INTENT(OUT)	:: kx
COMPLEX(KIND=8), DIMENSION(1:NY), INTENT(OUT)	:: ky
CHARACTER*100, INTENT(OUT)	:: name_config
INTEGER(kind=4)	:: i,j

!$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
! Save x grid points in text format
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)
! Save y grid points in text format
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)


END SUBROUTINE getgrid
</div>
</div>


<div class="toccolours mw-collapsible mw-collapsed" style="width:850px">
{{NumBlk|||{{EquationRef| G}}}} A Fortran subroutine to get the initial data to solve the 2D Klein-Gordon equation for [https://github.com/openmichigan/PSNM/blob/master/KleinGordon/Programs/KleinGordon2dThreadFFT/initialdata.f90 Code download]
<div class="mw-collapsible-content">
<syntaxhighlight lang="matlab">
SUBROUTINE initialdata(Nx,Ny,x,y,u,uold)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine gets initial data for nonlinear Klein-Gordon equation
! in 2 dimensions
! u_{tt}-u_{xx}+u_{yy}+u=Es*u^3+
!
! The boundary conditions are u(x=-Lx*\pi,y)=u(x=Lx*\pi,y),
! u(x,y=-Ly*\pi)=u(x,y=Ly*\pi)
! The initial condition is u=0.5*exp(-x^2-y^2)*sin(10*x+12*y)
!
! INPUT
!
! .. Parameters ..
! Nx = number of modes in x - power of 2 for FFT
! Ny = number of modes in y - power of 2 for FFT
! .. Vectors ..
! x = x locations
! y = y locations
!
! OUTPUT
!
! .. Arrays ..
! u = initial solution
! uold = approximate solution based on time derivative of
! initial solution
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
!
! 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
! OpenMP library
USE omp_lib	
IMPLICIT NONE	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny
REAL(KIND=8), DIMENSION(1:NX), INTENT(IN) :: x
REAL(KIND=8), DIMENSION(1:NY), INTENT(IN) :: y
COMPLEX(KIND=8), DIMENSION(1:NX,1:NY), INTENT(OUT)	:: u,uold
INTEGER(kind=4)	:: i,j
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
u(1:Nx,j)=0.5d0*exp(-1.0d0*(x(1:Nx)**2 +y(j)**2))*&
sin(10.0d0*x(1:Nx)+12.0d0*y(j))
END DO
!$OMP END PARALLEL DO
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
uold(1:Nx,j)=0.5d0*exp(-1.0d0*(x(1:Nx)**2 +y(j)**2))*&
sin(10.0d0*x(1:Nx)+12.0d0*y(j))
END DO
!$OMP END PARALLEL DO

END SUBROUTINE initialdata

( H)

一个 Fortran 程序来保存二维克莱因-戈登方程的解中的一个场代码下载

SUBROUTINE savedata(Nx,Ny,plotnum,name_config,field)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine saves a two 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
! 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
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y 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
IMPLICIT NONE	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny
INTEGER(KIND=4), INTENT(IN)	:: plotnum
REAL(KIND=8), DIMENSION(1:NX,1:NY), INTENT(IN)	:: field
CHARACTER*100, INTENT(IN)	:: name_config
INTEGER(kind=4)	:: i,j,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'	
INQUIRE( iolength=iol )	field(1,1)
OPEN(unit=11,FILE=number_file,form="unformatted",&
access="direct",recl=iol)
count=1
DO j=1,Ny
DO i=1,Nx
WRITE(11,rec=count) field(i,j)
count=count+1
END DO
END DO
CLOSE(11)

END SUBROUTINE savedata

( I)

一个 Fortran 子程序来更新求解二维克莱因-戈登方程时的数组代码下载

SUBROUTINE savedata(Nx,Ny,plotnum,name_config,field)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine saves a two 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
! 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
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y 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
IMPLICIT NONE	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny
INTEGER(KIND=4), INTENT(IN)	:: plotnum
REAL(KIND=8), DIMENSION(1:NX,1:NY), INTENT(IN)	:: field
CHARACTER*100, INTENT(IN)	:: name_config
INTEGER(kind=4)	:: i,j,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'	
INQUIRE( iolength=iol )	field(1,1)
OPEN(unit=11,FILE=number_file,form="unformatted",&
access="direct",recl=iol)
count=1
DO j=1,Ny
DO i=1,Nx
WRITE(11,rec=count) field(i,j)
count=count+1
END DO
END DO
CLOSE(11)

END SUBROUTINE savedata

( J)

一个 Fortran 子程序来计算求解二维克莱因-戈登方程时的能量代码下载

SUBROUTINE enercalc(Nx,Ny,planfxy,planbxy,dt,Es,enkin,enstr,&
enpot,en,kx,ky,temp1,temp2,v,vold,u,uold)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine program calculates the energy for the nonlinear
! Klein-Gordon equation in 2 dimensions
! u_{tt}-u_{xx}+u_{yy}+u=Es*|u|^2u
!
! The energy density is given by
! 0.5u_t^2+0.5u_x^2+0.5u_y^2+0.5u^2+Es*0.25u^4
!
! INPUT
!
! .. Scalars ..
! Nx = number of modes in x - power of 2 for FFT
! Ny = number of modes in y - power of 2 for FFT
! planfxy = Forward 2d fft plan
! planbxy = Backward 2d fft plan
! dt = timestep
! Es = +1 for focusing, -1 for defocusing
! .. Arrays ..
! u = approximate solution
! v = Fourier transform of approximate solution
! uold = approximate solution
! vold = Fourier transform of approximate solution
! temp1 = array to hold temporary values
! temp2 = array to hold temporary values
! .. Vectors ..
! kx = fourier frequencies in x direction
! ky = fourier frequencies in y direction
!
! OUTPUT
!
! .. Scalars ..
! enkin = Kinetic energy
! enstr = Strain energy
! enpot = Potential energy
! en = Total energy
!
! LOCAL VARIABLES
!
! .. Scalars ..
! j = loop counter in y direction
!
! 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
USE omp_lib	
IMPLICIT NONE	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny
REAL(KIND=8), INTENT(IN)	:: dt,Es
INTEGER(KIND=8), INTENT(IN)	:: planfxy	
INTEGER(KIND=8), INTENT(IN)	:: planbxy	
COMPLEX(KIND=8), DIMENSION(1:Nx),INTENT(IN)	:: kx
COMPLEX(KIND=8), DIMENSION(1:Ny),INTENT(IN)	:: ky
COMPLEX(KIND=8), DIMENSION(1:Nx,1:Ny),INTENT(IN)	:: u,v,uold,vold
COMPLEX(KIND=8), DIMENSION(1:Nx,1:Ny),INTENT(INOUT)	:: temp1,temp2
REAL(KIND=8), INTENT(OUT)	:: enkin,enstr
REAL(KIND=8), INTENT(OUT)	:: enpot,en	
INTEGER(KIND=4)	:: j

!.. Strain energy ..
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
temp1(1:Nx,j)=0.5d0*kx(1:Nx)*(vold(1:Nx,j)+v(1:Nx,j))
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planbxy,temp1(1:Nx,1:Ny),temp2(1:Nx,1:Ny))
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
temp1(1:Nx,j)=abs(temp2(1:Nx,j)/REAL(Nx*Ny,kind(0d0)))**2
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planfxy,temp1(1:Nx,1:Ny),temp2(1:Nx,1:Ny))
enstr=0.5d0*REAL(abs(temp2(1,1)),kind(0d0))/REAL(Nx*Ny,kind(0d0))
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
temp1(1:Nx,j)=0.5d0*ky(j)*(vold(1:Nx,j)+v(1:Nx,j))
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planbxy,temp1(1:Nx,1:Ny),temp2(1:Nx,1:Ny))
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
temp1(1:Nx,j)=abs(temp2(1:Nx,j)/REAL(Nx*Ny,kind(0d0)))**2
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planfxy,temp1(1:Nx,1:Ny),temp2(1:Nx,1:Ny))
enstr=enstr+0.5d0*REAL(abs(temp2(1,1)),kind(0d0))/REAL(Nx*Ny,kind(0d0))

! .. Kinetic Energy ..
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
temp1(1:Nx,j)=( abs(u(1:Nx,j)-uold(1:Nx,j))/dt )**2
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planfxy,temp1(1:Nx,1:Ny),temp2(1:Nx,1:Ny))
enkin=0.5d0*REAL(abs(temp2(1,1)),kind(0d0))/REAL(Nx*Ny,kind(0d0))

! .. Potential Energy ..
!$OMP PARALLEL DO PRIVATE(j) SCHEDULE(static)
DO j=1,Ny
temp1(1:Nx,j)=0.5d0*(abs((u(1:Nx,j)+uold(1:Nx,j))*0.50d0))**2&
-0.125d0*Es*(abs(u(1:Nx,j))**4+abs(uold(1:Nx,j))**4)
END DO
!$OMP END PARALLEL DO
CALL dfftw_execute_dft_(planfxy,temp1(1:Nx,1:Ny),temp2(1:Nx,1:Ny))
enpot=REAL(abs(temp2(1,1)),kind(0d0))/REAL(Nx*Ny,kind(0d0))

en=enpot+enkin+enstr

END SUBROUTINE enercalc

( K)

一个 Fortran 子程序来保存求解二维克莱因-戈登方程后的最终结果代码下载

SUBROUTINE saveresults(Nt,plotgap,time,en,enstr,enkin,enpot)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine saves the energy and times stored during the
! computation for the nonlinear Klein-Gordon equation
!
! INPUT
!
! .. Parameters ..
! Nx = number of modes in x - power of 2 for FFT
! Ny = number of modes in y - power of 2 for FFT
! .. Vectors ..
! time = times at which save data
! en = total energy
! enstr = strain energy
! enpot = potential energy
! enkin = kinetic energy
!
! OUTPUT
!
!
! LOCAL VARIABLES
!
! .. Scalars ..
! n = loop counter
! .. Arrays ..
! name_config = array to hold the filename
!
! REFERENCES
!
! ACKNOWLEDGEMENTS
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
!--------------------------------------------------------------------
! External routines required
!
! External libraries required
IMPLICIT NONE	
! Declare variables
INTEGER(kind=4), INTENT(IN)	:: plotgap,Nt
REAL(KIND=8), DIMENSION(1+Nt/plotgap), INTENT(IN)	:: enpot, enkin
REAL(KIND=8), DIMENSION(1+Nt/plotgap), INTENT(IN)	:: en,enstr,time
INTEGER(kind=4)	:: j
CHARACTER*100	:: name_config

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

name_config = 'enkin.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt/plotgap
WRITE(11,*) enkin(j)
END DO
CLOSE(11)

name_config = 'enpot.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt/plotgap
WRITE(11,*) enpot(j)
END DO
CLOSE(11)

name_config = 'enstr.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,1+Nt/plotgap
WRITE(11,*) enstr(j)
END DO
CLOSE(11)

END SUBROUTINE saveresults

( L)

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

#define the complier
COMPILER = mpif90
# compilation settings, optimization, precision, parallelization
FLAGS = -O0 -mp

# libraries
LIBS = -L${FFTW_LINK} -lfftw3_threads -lfftw3 -lm
# source list for main program
SOURCES = KgSemiImp2d.f90 initialdata.f90 savedata.f90 getgrid.f90 \
storeold.f90 saveresults.f90 enercalc.f90

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

clean:
rm *.o

( M)

一个 Matlab 程序，用于绘制清单 E 生成的场代码下载

% A program to create a video of the computed results

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');
TIME=load('./tdata.dat');
% find number of grid points
Nx=length(X);
Ny=length(Y);

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

nplots=length(TIME);

for i =1:nplots
    %
    % Open file and dataset using the default properties.
    %
    FILE=['./data/u',num2str(9999999+i),'.datbin'];
    FILEPIC=['./data/pic',num2str(9999999+i),'.jpg'];
    fid=fopen(FILE,'r');
    [fname,mode,mformat]=fopen(fid);
    u=fread(fid,Nx*Ny,'real*8');
    u=reshape(u,Nx,Ny);
    % close files
    fclose(fid);
    %
    % Plot data on the screen.
    %
    figure(2);clf; mesh(xx,yy,real(u)); xlabel x; ylabel y;
    title(['Time ',num2str(TIME(i))]); colorbar; axis square;
    drawnow; frame=getframe(2); saveas(2,FILEPIC,'jpg');
end

使用 2DECOMP&FFT 的三维 MPI Fortran 程序

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

我们现在为三维非线性克莱因-戈登方程提供一个程序。该程序使用与二维代码相同的子程序结构。为了使程序易于重用，创建了清单 U 中列出的子程序来读取 INPUTFILE，该文件指定了程序要使用的参数，因此程序无需每次运行时都重新编译。为了使程序更好地扩展，保存傅里叶频率和网格点的数组也被分解，以便只创建并存储在每个处理器上使用的数组部分。该程序在清单 W 中列出，要使用此程序，只需使用 MPI 包装器和 gfortran（其他编译器，如 Intel fortran、PGI fortran 或 IBM XL fortran 也应该可以工作）编译它，然后在存储 INPUTFILE 和数据的目录中运行它。示例提交脚本在清单 Y 中，这需要一个名为“inputfile”的环境变量，该变量设置为指定问题参数（如网格大小和时间步数）的 INPUTFILE 的位置。

另一个补充是一个简短的后处理程序，用于创建头文件，以便可以利用 bov（值块）格式，该格式允许使用并行可视化软件 VisIt。VisIt 程序可以从 https://wci.llnl.gov/codes/visit/home.html 下载。此程序也可以在笔记本电脑、台式机以及并行计算机集群上运行。有关使用 VisIt 的文档，请参见 https://wci.llnl.gov/codes/visit/manuals.html 和 http://www.visitusers.org/index.php?title=Main_Page。有关如何远程使用 VisIt 的简短视频教程，请参见 http://cac.engin.umich.edu/resources/software/visit.html。另一种创建图像的程序是 Paraview，它可以从 http://www.paraview.org/ 下载。要使用 Paraview，可以使用 xdmf 格式头文件使 Paraview 读取二进制数据输出文件。有关使用 Paraview 进行渲染的简短视频教程，请参见 http://www.youtube.com/watch?v=rqCXFr73sfk&list=PL4FFYwZhtT0BlDBtjsSkQ4LcM1v9K24No。一个 Fortran 程序来创建这些头文件，在清单 X 中。

( N)

一个 Fortran 程序来求解三维克莱因-戈登方程代码下载

!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This program solves nonlinear Klein-Gordon equation in 3 dimensions
! u_{tt}-(u_{xx}+u_{yy}+u_{zz})+u=Es*|u|^2u
! using a second order implicit-explicit time stepping scheme.
!
! The boundary conditions are u(x=-Lx*pi,y,z)=u(x=Lx*\pi,y,z),
! u(x,y=-Ly*pi,z)=u(x,y=Ly*pi,z),u(x,y,z=-Ly*pi)=u(x,y,z=Ly*pi),
! The initial condition is u=0.5*exp(-x^2-y^2-z^2)*sin(10*x+12*y)
!
! .. 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
! ierr = error code
! plotnum = number of plot
! myid = Process id
! 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
! .. Arrays ..
! unew = approximate solution
! vnew = Fourier transform of approximate solution
! u = approximate solution
! v = Fourier transform of approximate solution
! uold = approximate solution
! vold = Fourier transform of approximate solution
! nonlin = nonlinear term, u^3
! nonlinhat = Fourier transform of nonlinear term, u^3
! .. 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
! en = total energy
! enstr = strain energy
! enpot = potential energy
! enkin = kinetic energy
! name_config = array to store filename for data to be saved
! fftfxyz = array to setup 2D Fourier transform
! fftbxyz = array to setup 2D Fourier transform
! .. 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
! getgrid.f90 -- Get initial grid of points
! initialdata.f90 -- Get initial data
! enercalc.f90 -- Subroutine to calculate the energy
! savedata.f90 -- Save initial data
! storeold.f90 -- Store old data
! External libraries required
! 2DECOMP&FFT -- Domain decomposition and Fast Fourier Library
! (http://www.2decomp.org/index.html)
! MPI library
PROGRAM Kg
IMPLICIT NONE
USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
INCLUDE 'mpif.h'
! Declare variables
INTEGER(kind=4)	:: Nx, Ny, Nz, Nt, plotgap	
REAL(kind=8), PARAMETER	:: &
pi=3.14159265358979323846264338327950288419716939937510d0
REAL(kind=8)	:: Lx,Ly,Lz,Es,dt,starttime,modescalereal	
COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE	:: kx,ky,kz
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: x,y,z	
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE:: u,nonlin
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE:: v,nonlinhat
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE:: uold
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE:: vold
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE:: unew
COMPLEX(kind=8), DIMENSION(:,:,:), ALLOCATABLE:: vnew
REAL(kind=8), DIMENSION(:,:,:), ALLOCATABLE :: savearray
REAL(kind=8), DIMENSION(:), ALLOCATABLE	:: time,enkin,enstr,enpot,en
INTEGER(kind=4)	:: ierr,i,j,k,n,allocatestatus,myid,numprocs
INTEGER(kind=4)	:: start, finish, count_rate, plotnum
TYPE(DECOMP_INFO)	:: decomp
INTEGER(kind=MPI_OFFSET_KIND) :: filesize, disp
INTEGER(kind=4)	:: p_row=0, p_col=0	
CHARACTER*100	:: name_config
! 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)

CALL readinputfile(Nx,Ny,Nz,Nt,plotgap,Lx,Ly,Lz, &
Es,DT,starttime,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
ALLOCATE(kx(decomp%zst(1):decomp%zen(1)),&
ky(decomp%zst(2):decomp%zen(2)),&
kz(decomp%zst(3):decomp%zen(3)),&
x(decomp%xst(1):decomp%xen(1)),&
y(decomp%xst(2):decomp%xen(2)),&
z(decomp%xst(3):decomp%xen(3)),&
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)),&
nonlin(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
nonlinhat(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
uold(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
vold(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
unew(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
vnew(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
savearray(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
time(1:1+Nt/plotgap),enkin(1:1+Nt/plotgap),&
enstr(1:1+Nt/plotgap),enpot(1:1+Nt/plotgap),&
en(1:1+Nt/plotgap),stat=allocatestatus)	
IF (allocatestatus .ne. 0) stop
IF (myid.eq.0) THEN
PRINT *,'allocated arrays'
END IF	
! setup fourier frequencies
CALL getgrid(myid,Nx,Ny,Nz,Lx,Ly,Lz,pi,name_config,x,y,z,kx,ky,kz,decomp)
IF (myid.eq.0) THEN
PRINT *,'Setup grid and fourier frequencies'
END IF
CALL initialdata(Nx,Ny,Nz,x,y,z,u,uold,decomp)
plotnum=1	
name_config = 'data/u'
savearray=REAL(u)
CALL savedata(Nx,Ny,Nz,plotnum,name_config,savearray,decomp)

CALL decomp_2d_fft_3d(u,v,DECOMP_2D_FFT_FORWARD)
CALL decomp_2d_fft_3d(uold,vold,DECOMP_2D_FFT_FORWARD)

modescalereal=1.0d0/REAL(Nx,KIND(0d0))
modescalereal=modescalereal/REAL(Ny,KIND(0d0))
modescalereal=modescalereal/REAL(Nz,KIND(0d0))

CALL enercalc(myid,Nx,Ny,Nz,dt,Es,modescalereal,&
enkin(plotnum),enstr(plotnum),&
enpot(plotnum),en(plotnum),&
kx,ky,kz,nonlin,nonlinhat,&
v,vold,u,uold,decomp)

IF (myid.eq.0) THEN
PRINT *,'Got initial data, starting timestepping'
END IF
time(plotnum)=0.0d0+starttime
   CALL system_clock(start,count_rate)	
DO n=1,Nt	
DO k=decomp%xst(3),decomp%xen(3)
DO j=decomp%xst(2),decomp%xen(2)
DO i=decomp%xst(1),decomp%xen(1)
nonlin(i,j,k)=(abs(u(i,j,k))*2)*u(i,j,k)
END DO
END DO
END DO
CALL decomp_2d_fft_3d(nonlin,nonlinhat,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)
vnew(i,j,k)=&
( 0.25*(kx(i)*kx(i) + ky(j)*ky(j)+ kz(k)*kz(k)-1.0d0)&
*(2.0d0*v(i,j,k)+vold(i,j,k))&
+(2.0d0*v(i,j,k)-vold(i,j,k))/(dt*dt)&
+Es*nonlinhat(i,j,k) )&
/(1/(dt*dt)-0.25*(kx(i)*kx(i)+ ky(j)*ky(j)+ kz(k)*kz(k)-1.0d0))
END DO
END DO
END DO
CALL decomp_2d_fft_3d(vnew,unew,DECOMP_2D_FFT_BACKWARD)
! normalize result
DO k=decomp%xst(3),decomp%xen(3)
DO j=decomp%xst(2),decomp%xen(2)
DO i=decomp%xst(1),decomp%xen(1)
unew(i,j,k)=unew(i,j,k)*modescalereal
END DO
END DO
END DO
IF (mod(n,plotgap)==0) THEN
plotnum=plotnum+1
time(plotnum)=n*dt+starttime
IF (myid.eq.0) THEN
PRINT *,'time',n*dt+starttime
END IF
CALL enercalc(myid,Nx,Ny,Nz,dt,Es,modescalereal,&
enkin(plotnum),enstr(plotnum),&
enpot(plotnum),en(plotnum),&
kx,ky,kz,nonlin,nonlinhat,&
vnew,v,unew,u,decomp)
savearray=REAL(unew,kind(0d0))
CALL savedata(Nx,Ny,Nz,plotnum,name_config,savearray,decomp)
END IF
! .. Update old values ..
CALL storeold(Nx,Ny,Nz,unew,u,uold,vnew,v,vold,decomp)
END DO
CALL system_clock(finish,count_rate)
IF (myid.eq.0) THEN
PRINT *,'Finished time stepping'
PRINT*,'Program took ',&
REAL(finish-start,kind(0d0))/REAL(count_rate,kind(0d0)),&
'for Time stepping'
CALL saveresults(Nt,plotgap,time,en,enstr,enkin,enpot)	
! Save times at which output was made in text format
PRINT *,'Saved data'
END IF
CALL decomp_2d_fft_finalize
   CALL decomp_2d_finalize

DEALLOCATE(kx,ky,kz,x,y,z,u,v,nonlin,nonlinhat,savearray,&
uold,vold,unew,vnew,time,enkin,enstr,enpot,en,&
stat=allocatestatus)	
IF (allocatestatus .ne. 0) STOP
IF (myid.eq.0) THEN
PRINT *,'Deallocated arrays'
PRINT *,'Program execution complete'
END IF
CALL MPI_FINALIZE(ierr)

END PROGRAM Kg

( O)

一个 Fortran 子程序来获取在三维克莱因-戈登方程上求解的网格代码下载

SUBROUTINE getgrid(myid,Nx,Ny,Nz,Lx,Ly,Lz,pi,name_config,x,y,z,kx,ky,kz,decomp)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine gets grid points and fourier frequencies for a
! pseudospectral simulation of the 2D nonlinear Klein-Gordon equation
!
! u_{tt}-(u_{xx}+u_{yy}+u_{zz})+u=Es*u^3
!
! The boundary conditions are u(x=-Lx*pi,y,z)=u(x=Lx*\pi,y,z),
! u(x,y=-Ly*pi,z)=u(x,y=Ly*pi,z),u(x,y,z=-Ly*pi)=u(x,y,z=Ly*pi),
!
! INPUT
!
! .. Scalars ..
! Nx = number of modes in x - power of 2 for FFT
! Ny = number of modes in y - power of 2 for FFT
! Ny = number of modes in z - power of 2 for FFT
! pi = 3.142....
! Lx = width of box in x direction
! Ly = width of box in y direction
! Lz = width of box in z direction
! myid = processor id
! .. Special Structures ..
! decomp = contains information on domain decomposition
! see http://www.2decomp.org/ for more information
!
! OUTPUT
!
! .. 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
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
! k = loop counter in z direction
!
! 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
IMPLICIT NONE
USE decomp_2d
INCLUDE 'mpif.h'
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: myid,Nx,Ny,Nz
REAL(kind=8), INTENT(IN)	:: Lx,Ly,Lz,pi
TYPE(DECOMP_INFO), INTENT(IN)	:: decomp
REAL(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1)), INTENT(OUT) :: x
REAL(KIND=8), DIMENSION(decomp%xst(2):decomp%xen(2)), INTENT(OUT) :: y
REAL(KIND=8), DIMENSION(decomp%xst(3):decomp%xen(3)), INTENT(OUT) :: z
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1)), INTENT(OUT):: kx
COMPLEX(KIND=8), DIMENSION(decomp%zst(2):decomp%zen(2)), INTENT(OUT):: ky
COMPLEX(KIND=8), DIMENSION(decomp%zst(3):decomp%zen(3)), INTENT(OUT):: kz
CHARACTER*100, INTENT(OUT)	:: name_config
INTEGER(kind=4)	:: i,j,k


DO i = 1,1+ Nx/2
IF ((i.GE.decomp%zst(1)).AND.(i.LE.decomp%zen(1))) THEN
kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/Lx
END IF
END DO
IF ((Nx/2 + 1 .GE.decomp%zst(1)).AND.(Nx/2 + 1 .LE.decomp%zen(1))) THEN
kx( Nx/2 + 1 ) = 0.0d0
ENDIF
DO i = Nx/2+2, Nx
IF ((i.GE.decomp%zst(1)).AND.(i.LE.decomp%zen(1))) THEN
Kx( i) = cmplx(0.0d0,-1.0d0)*REAL(1-i+Nx,KIND(0d0))/Lx
ENDIF
END DO
DO i=decomp%xst(1),decomp%xen(1)
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
IF ((j.GE.decomp%zst(2)).AND.(j.LE.decomp%zen(2))) THEN
ky(j)= cmplx(0.0d0,1.0d0)*REAL(j-1,kind(0d0))/Ly
END IF
END DO
IF ((Ny/2 + 1 .GE.decomp%zst(2)).AND.(Ny/2 + 1 .LE.decomp%zen(2))) THEN
ky( Ny/2 + 1 ) = 0.0d0
ENDIF
DO j = Ny/2+2, Ny
IF ((j.GE.decomp%zst(2)).AND.(j.LE.decomp%zen(2))) THEN
ky(j) = cmplx(0.0d0,-1.0d0)*REAL(1-j+Ny,KIND(0d0))/Ly
ENDIF
END DO
DO j=decomp%xst(2),decomp%xen(2)
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
IF ((k.GE.decomp%zst(3)).AND.(k.LE.decomp%zen(3))) THEN
kz(k)= cmplx(0.0d0,1.0d0)*REAL(k-1,kind(0d0))/Lz
END IF
END DO
IF ((Nz/2 + 1 .GE.decomp%zst(3)).AND.(Nz/2 + 1 .LE.decomp%zen(3))) THEN
kz( Nz/2 + 1 ) = 0.0d0
ENDIF
DO k = Nz/2+2, Nz
IF ((k.GE.decomp%zst(3)).AND.(k.LE.decomp%zen(3))) THEN
kz(k) = cmplx(0.0d0,-1.0d0)*REAL(1-k+Nz,KIND(0d0))/Lz
ENDIF
END DO
DO k=decomp%xst(3),decomp%xen(3)
z(k)=(-1.0d0 + 2.0d0*REAL(k-1,kind(0d0))/REAL(Nz,kind(0d0)))*pi*Lz
END DO
IF (myid.eq.0) THEN
! Save x grid points in text format
name_config = 'xcoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO i=1,Nx
WRITE(11,*) (-1.0d0 + 2.0d0*REAL(i-1,kind(0d0))/REAL(Nx,kind(0d0)))*pi*Lx
END DO
CLOSE(11)
! Save y grid points in text format
name_config = 'ycoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO j=1,Ny
WRITE(11,*) (-1.0d0 + 2.0d0*REAL(j-1,kind(0d0))/REAL(Ny,kind(0d0)))*pi*Ly
END DO
CLOSE(11)
! Save z grid points in text format
name_config = 'zcoord.dat'
OPEN(unit=11,FILE=name_config,status="UNKNOWN")
REWIND(11)
DO k=1,Nz
WRITE(11,*) (-1.0d0 + 2.0d0*REAL(k-1,kind(0d0))/REAL(Nz,kind(0d0)))*pi*Lz
END DO
CLOSE(11)
END IF
END SUBROUTINE getgrid

( P)

一个 Fortran 子程序来获取求解三维克莱因-戈登方程的初始数据代码下载

SUBROUTINE initialdata(Nx,Ny,Nz,x,y,z,u,uold,decomp)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine gets initial data for nonlinear Klein-Gordon equation
! in 3 dimensions
! u_{tt}-(u_{xx}+u_{yy}+u_{zz})+u=Es*u^3+
!
! The boundary conditions are u(x=-Lx*\pi,y,z)=u(x=Lx*\pi,y,z),
! u(x,y=-Ly*\pi,z)=u(x,y=Ly*\pi,z),u(x,y,z=-Ly*\pi)=u(x,y,z=Ly*\pi)
! The initial condition is u=0.5*exp(-x^2-y^2-z^2)*sin(10*x+12*y)
!
! INPUT
!
! .. 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
! .. Vectors ..
! x = x locations
! y = y locations
! z = z locations
! .. Special Structures ..
! decomp = contains information on domain decomposition
! see http://www.2decomp.org/ for more information
! OUTPUT
!
! .. Arrays ..
! u = initial solution
! uold = approximate solution based on time derivative of
! initial solution
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
! k = loop counter in z direction
!
! 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
IMPLICIT NONE
USE decomp_2d
INCLUDE 'mpif.h'
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny,Nz
TYPE(DECOMP_INFO), INTENT(IN)	:: decomp
REAL(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1)), INTENT(IN) :: x
REAL(KIND=8), DIMENSION(decomp%xst(2):decomp%xen(2)), INTENT(IN) :: y
REAL(KIND=8), DIMENSION(decomp%xst(3):decomp%xen(3)), INTENT(IN) :: z
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
    INTENT(OUT) :: u,uold
INTEGER(kind=4)	:: i,j,k

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.5d0*exp(-1.0d0*(x(i)**2 +y(j)**2+z(k)**2))!*&
!sin(10.0d0*x(i)+12.0d0*y(j))
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)
uold(i,j,k)=0.5d0*exp(-1.0d0*(x(i)**2 +y(j)**2+z(k)**2))!*&
!sin(10.0d0*x(i)+12.0d0*y(j))
END DO
END DO
END DO
END SUBROUTINE initialdata

( Q)

一个 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
IMPLICIT NONE
USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
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

( R)

一个 Fortran 子程序来更新求解三维克莱因-戈登方程时的数组代码下载

SUBROUTINE storeold(Nx,Ny,Nz,unew,u,uold,vnew,v,vold,decomp)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine copies arrays for a
! pseudospectral simulation of the 2D nonlinear Klein-Gordon equation
!
! u_{tt}-(u_{xx}+u_{yy}+u_{zz})+u=Es*u^3
!
! INPUT
!
! .. 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
! .. Arrays ..
! unew = approximate solution
! vnew = Fourier transform of approximate solution
! u = approximate solution
! v = Fourier transform of approximate solution
! uold = approximate solution
! vold = Fourier transform of approximate solution
! .. Special Structures ..
! decomp = contains information on domain decomposition
! see http://www.2decomp.org/ for more information
! OUTPUT
!
! u = approximate solution
! v = Fourier transform of approximate solution
! uold = approximate solution
! vold = Fourier transform of approximate solution
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
! k = loop counter in z direction
!
! 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
IMPLICIT NONE
USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
INCLUDE 'mpif.h'	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny,Nz
TYPE(DECOMP_INFO), INTENT(IN)	:: decomp
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)), INTENT(OUT):: uold
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)), INTENT(OUT):: vold
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)), INTENT(INOUT):: v
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)), INTENT(INOUT):: u
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)), INTENT(IN):: unew
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)), INTENT(IN):: vnew
INTEGER(kind=4)	:: i,j,k

DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
vold(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)
uold(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)
u(i,j,k)=unew(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)
v(i,j,k)=vnew(i,j,k)
END DO
END DO
END DO
END SUBROUTINE storeold

( S)

一个 Fortran 子程序来计算求解三维克莱因-戈登方程时的能量代码下载

SUBROUTINE enercalc(myid,Nx,Ny,Nz,dt,Es,modescalereal,enkin,enstr,&
enpot,en,kx,ky,kz,tempu,tempv,v,vold,u,uold,decomp)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine program calculates the energy for the nonlinear
! Klein-Gordon equation in 3 dimensions
! u_{tt}-(u_{xx}+u_{yy}+u_{zz})+u=Es*|u|^2u
!
! The energy density is given by
! 0.5u_t^2+0.5u_x^2+0.5u_y^2+0.5u_z^2+0.5u^2+Es*0.25u^4
!
! 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
! dt = timestep
! Es = +1 for focusing, -1 for defocusing
! modescalereal = parameter to scale after doing backward FFT
! myid = Process id
! .. Arrays ..
! u = approximate solution
! v = Fourier transform of approximate solution
! uold = approximate solution
! vold = Fourier transform of approximate solution
! tempu = array to hold temporary values - real space
! tempv = array to hold temporary values - fourier space
! .. Vectors ..
! kx = fourier frequencies in x direction
! ky = fourier frequencies in y direction
! kz = fourier frequencies in z direction
! .. Special Structures ..
! decomp = contains information on domain decomposition
! see http://www.2decomp.org/ for more information
! OUTPUT
!
! .. Scalars ..
! enkin = Kinetic energy
! enstr = Strain energy
! enpot = Potential energy
! en = Total energy
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
! k = loop counter in z direction
!
! 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
IMPLICIT NONE
USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
INCLUDE 'mpif.h'	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny,Nz,myid
REAL(KIND=8), INTENT(IN)	:: dt,Es,modescalereal
TYPE(DECOMP_INFO), INTENT(IN)	:: decomp
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1)),INTENT(IN)	:: kx
COMPLEX(KIND=8), DIMENSION(decomp%zst(2):decomp%zen(2)),INTENT(IN)	:: ky
COMPLEX(KIND=8), DIMENSION(decomp%zst(3):decomp%zen(3)),INTENT(IN)	:: kz
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
    INTENT(IN)	:: u,uold
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
             INTENT(IN)	:: v,vold
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
    INTENT(INOUT) :: tempu
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
             INTENT(INOUT):: tempv
REAL(KIND=8), INTENT(OUT)	:: enkin,enstr
REAL(KIND=8), INTENT(OUT)	:: enpot,en	
INTEGER(KIND=4)	:: i,j,k

!.. Strain energy ..
DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
tempv(i,j,k)=0.5d0*kx(i)*(vold(i,j,k)+v(i,j,k))
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempv,tempu,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)
tempu(i,j,k)=abs(tempu(i,j,k)*modescalereal)**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enstr=0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
tempv(i,j,k)=0.5d0*ky(j)*(vold(i,j,k)+v(i,j,k))
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempv,tempu,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)
tempu(i,j,k)=abs(tempu(i,j,k)*modescalereal)**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enstr=enstr+0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
tempv(i,j,k)=0.5d0*kz(k)*(vold(i,j,k)+v(i,j,k))
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempv,tempu,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)
tempu(i,j,k)=abs(tempu(i,j,k)*modescalereal)**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enstr=enstr+0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
! .. Kinetic Energy ..
DO k=decomp%xst(3),decomp%xen(3)
DO j=decomp%xst(2),decomp%xen(2)
DO i=decomp%xst(1),decomp%xen(1)
tempu(i,j,k)=( abs(u(i,j,k)-uold(i,j,k))/dt )**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enkin=0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
! .. Potential Energy ..
DO k=decomp%xst(3),decomp%xen(3)
DO j=decomp%xst(2),decomp%xen(2)
DO i=decomp%xst(1),decomp%xen(1)
tempu(i,j,k)=0.5d0*(abs((u(i,j,k)+uold(i,j,k))*0.50d0))**2&
-0.125d0*Es*(abs(u(i,j,k))**4+abs(uold(i,j,k))**4)
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enpot=REAL(abs(tempv(1,1,1)),kind(0d0))	
en=enpot+enkin+enstr
END IF
END SUBROUTINE enercalc

( T)

一个 Fortran 子程序来保存求解三维克莱因-戈登方程后的最终结果代码下载

SUBROUTINE enercalc(myid,Nx,Ny,Nz,dt,Es,modescalereal,enkin,enstr,&
enpot,en,kx,ky,kz,tempu,tempv,v,vold,u,uold,decomp)
!--------------------------------------------------------------------
!
!
! PURPOSE
!
! This subroutine program calculates the energy for the nonlinear
! Klein-Gordon equation in 3 dimensions
! u_{tt}-(u_{xx}+u_{yy}+u_{zz})+u=Es*|u|^2u
!
! The energy density is given by
! 0.5u_t^2+0.5u_x^2+0.5u_y^2+0.5u_z^2+0.5u^2+Es*0.25u^4
!
! 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
! dt = timestep
! Es = +1 for focusing, -1 for defocusing
! modescalereal = parameter to scale after doing backward FFT
! myid = Process id
! .. Arrays ..
! u = approximate solution
! v = Fourier transform of approximate solution
! uold = approximate solution
! vold = Fourier transform of approximate solution
! tempu = array to hold temporary values - real space
! tempv = array to hold temporary values - fourier space
! .. Vectors ..
! kx = fourier frequencies in x direction
! ky = fourier frequencies in y direction
! kz = fourier frequencies in z direction
! .. Special Structures ..
! decomp = contains information on domain decomposition
! see http://www.2decomp.org/ for more information
! OUTPUT
!
! .. Scalars ..
! enkin = Kinetic energy
! enstr = Strain energy
! enpot = Potential energy
! en = Total energy
!
! LOCAL VARIABLES
!
! .. Scalars ..
! i = loop counter in x direction
! j = loop counter in y direction
! k = loop counter in z direction
!
! 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
IMPLICIT NONE
USE decomp_2d
USE decomp_2d_fft
USE decomp_2d_io
INCLUDE 'mpif.h'	
! Declare variables
INTEGER(KIND=4), INTENT(IN)	:: Nx,Ny,Nz,myid
REAL(KIND=8), INTENT(IN)	:: dt,Es,modescalereal
TYPE(DECOMP_INFO), INTENT(IN)	:: decomp
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1)),INTENT(IN)	:: kx
COMPLEX(KIND=8), DIMENSION(decomp%zst(2):decomp%zen(2)),INTENT(IN)	:: ky
COMPLEX(KIND=8), DIMENSION(decomp%zst(3):decomp%zen(3)),INTENT(IN)	:: kz
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
    INTENT(IN)	:: u,uold
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
             INTENT(IN)	:: v,vold
COMPLEX(KIND=8), DIMENSION(decomp%xst(1):decomp%xen(1),&
    decomp%xst(2):decomp%xen(2),&
    decomp%xst(3):decomp%xen(3)),&
    INTENT(INOUT) :: tempu
COMPLEX(KIND=8), DIMENSION(decomp%zst(1):decomp%zen(1),&
             decomp%zst(2):decomp%zen(2),&
             decomp%zst(3):decomp%zen(3)),&
             INTENT(INOUT):: tempv
REAL(KIND=8), INTENT(OUT)	:: enkin,enstr
REAL(KIND=8), INTENT(OUT)	:: enpot,en	
INTEGER(KIND=4)	:: i,j,k

!.. Strain energy ..
DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
tempv(i,j,k)=0.5d0*kx(i)*(vold(i,j,k)+v(i,j,k))
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempv,tempu,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)
tempu(i,j,k)=abs(tempu(i,j,k)*modescalereal)**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enstr=0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
tempv(i,j,k)=0.5d0*ky(j)*(vold(i,j,k)+v(i,j,k))
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempv,tempu,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)
tempu(i,j,k)=abs(tempu(i,j,k)*modescalereal)**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enstr=enstr+0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
DO k=decomp%zst(3),decomp%zen(3)
DO j=decomp%zst(2),decomp%zen(2)
DO i=decomp%zst(1),decomp%zen(1)
tempv(i,j,k)=0.5d0*kz(k)*(vold(i,j,k)+v(i,j,k))
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempv,tempu,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)
tempu(i,j,k)=abs(tempu(i,j,k)*modescalereal)**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enstr=enstr+0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
! .. Kinetic Energy ..
DO k=decomp%xst(3),decomp%xen(3)
DO j=decomp%xst(2),decomp%xen(2)
DO i=decomp%xst(1),decomp%xen(1)
tempu(i,j,k)=( abs(u(i,j,k)-uold(i,j,k))/dt )**2
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enkin=0.5d0*REAL(abs(tempv(1,1,1)),kind(0d0))
END IF
! .. Potential Energy ..
DO k=decomp%xst(3),decomp%xen(3)
DO j=decomp%xst(2),decomp%xen(2)
DO i=decomp%xst(1),decomp%xen(1)
tempu(i,j,k)=0.5d0*(abs((u(i,j,k)+uold(i,j,k))*0.50d0))**2&
-0.125d0*Es*(abs(u(i,j,k))**4+abs(uold(i,j,k))**4)
END DO
END DO
END DO
CALL decomp_2d_fft_3d(tempu,tempv,DECOMP_2D_FFT_FORWARD)
IF(myid.eq.0) THEN
enpot=REAL(abs(tempv(1,1,1)),kind(0d0))	
en=enpot+enkin+enstr
END IF
END SUBROUTINE enercalc

( U)

用于读取求解三维 Klein-Gordon 方程参数的 Fortran 子程序代码下载

SUBROUTINE readinputfile(Nx,Ny,Nz,Nt,plotgap,Lx,Ly,Lz, &
Es,DT,starttime,myid,ierr)
!--------------------------------------------------------------------------------
!
!
! PURPOSE
!
! Read inputfile intialize parameters, which are stocked in the Input File
!
! .. INPUT ..
! Nx = number of modes in the x direction
! Ny = number of modes in the y direction
! Nz = number of modes in the z direction
! Nt = the number of timesteps
! plotgap = the number of timesteps to take before plotting
! myid = number of MPI process
! ierr = MPI error output variable
! Lx = size of the periodic domain of computation in x direction
! Ly = size of the periodic domain of computation in y direction
! Lz = size of the periodic domain of computation in z direction
! DT = the time step
! starttime = initial time of computation
! InputFileName = name of the Input File
! REFERENCES
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
!---------------------------------------------------------------------------------
! EXTERNAL ROUTINES REQUIRED
IMPLICIT NONE
INCLUDE 'mpif.h'
! .. Scalar Arguments ..
INTEGER(KIND=4), INTENT(IN)	:: myid	
INTEGER(KIND=4), INTENT(OUT)	:: Nx,Ny,Nz,Nt
INTEGER(KIND=4), INTENT(OUT)	:: plotgap, ierr
REAL(KIND=8), INTENT(OUT)	:: Lx, Ly, Lz, DT, starttime, Es
! .. Local scalars ..
INTEGER(KIND=4)	:: stat
! .. Local Arrays ..
CHARACTER*40	:: InputFileName
INTEGER(KIND=4), DIMENSION(1:5)	:: intcomm
REAL(KIND=8), DIMENSION(1:6)	:: dpcomm

IF(myid.eq.0) THEN
CALL GET_ENVIRONMENT_VARIABLE(NAME="inputfile",VALUE=InputFileName, STATUS=stat)
IF(stat.NE.0) THEN
PRINT*,"Set environment variable inputfile to the name of the"
PRINT*,"file where the simulation parameters are set"
STOP
END IF
OPEN(unit=11,FILE=trim(InputFileName),status="OLD")
REWIND(11)
READ(11,*) intcomm(1), intcomm(2), intcomm(3), intcomm(4), intcomm(5), &
dpcomm(1), dpcomm(2), dpcomm(3), dpcomm(4), dpcomm(5), dpcomm(6)
CLOSE(11)
PRINT *,"NX ",intcomm(1)
PRINT *,"NY ",intcomm(2)
PRINT *,"NZ ",intcomm(3)
PRINT *,"NT ",intcomm(4)
PRINT *,"plotgap ",intcomm(5)
PRINT *,"Lx ",dpcomm(1)
PRINT *,"Ly ",dpcomm(2)
PRINT *,"Lz ",dpcomm(3)
PRINT *,"Es ",dpcomm(4)	
PRINT *,"Dt ",dpcomm(5)
PRINT *,"strart time ",dpcomm(6)
PRINT *,"Read inputfile"
END IF
CALL MPI_BCAST(dpcomm,6,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
CALL MPI_BCAST(intcomm,5,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)

Nx=intcomm(1)
Ny=intcomm(2)
Nz=intcomm(3)
Nt=intcomm(4)
plotgap=intcomm(5)
Lx=dpcomm(1)
Ly=dpcomm(2)
Lz=dpcomm(3)
Es=dpcomm(4)
DT=dpcomm(5)
starttime=dpcomm(6)

END SUBROUTINE readinputfile

( V)

用于编译列表 N 中的 MPI 程序的示例 Makefile 代码下载

	# All settings here for use on FLUX, a cluster at the University of Michigan
# Center for Advanced Computing (CAC), using INTEL nehalem hardware,
# Need to load fftw module

COMPILER = mpif90
decompdir=../2decomp_fft
# compilation settings, optimization, precision, parallelization
FLAGS = -O0 -fltconsistency
LIBS = -L${FFTW_LINK} -lfftw3

DECOMPLIB = -I${decompdir}/include -L${decompdir}/lib -l2decomp_fft


# libraries
# source list for main program
SOURCES = KgSemiImp3d.f90 initialdata.f90 savedata.f90 getgrid.f90 \
storeold.f90 saveresults.f90 enercalc.f90 readinputfile.f90

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


clean:
rm -f *.o
clobber:
rm -f Kg

( W)

用于在求解三维 Klein-Gordon 方程后创建 BOV (Brick of Values) 头文件的 Fortran 子程序代码下载

PROGRAM BovCreate
!--------------------------------------------------------------------------------
! .. Purpose ..
! BovCreate is a postprocessing program which creates header files for VisIt
! It uses the INPUTFILE and assumes that the filenames in the program are
! consistent with those in the current file.
!
! .. PARAMETERS .. INITIALIZED IN INPUTFILE
! time = start time of the simulation
! Nx = power of two, number of modes in the x direction
! Ny = power of two, number of modes in the y direction
! Nz = power of two, number of modes in the z direction
! Nt = the number of timesteps
! plotgap = the number of timesteps to take before plotting
! Lx = definition of the periodic domain of computation in x direction
! Ly = definition of the periodic domain of computation in y direction
! Lz = definition of the periodic domain of computation in z direction
! Es = focusing or defocusing
! Dt = the time step
!
! REFERENCES
!
! ACCURACY
!
! ERROR INDICATORS AND WARNINGS
!
! FURTHER COMMENTS
!------------------------------------------------------------------------------------
! EXTERNAL ROUTINES REQUIRED
IMPLICIT NONE
! .. Scalar Arguments ..
INTEGER(KIND=4)	:: Nx, Ny, Nz, Nt, plotgap
REAL(KIND=8)	:: Lx, Ly, Lz, DT, time, Es
! .. Local scalars ..
INTEGER(KIND=4)	:: stat,plotnum,ind,n,numplots
! .. Local Arrays ..
CHARACTER*50	:: InputFileName, OutputFileName, OutputFileName2
CHARACTER*10	:: number_file
InputFileName='INPUTFILE'
OPEN(unit=11,FILE=trim(InputFileName),status="OLD")
REWIND(11)
READ(11,*) Nx, Ny, Nz, Nt, plotgap, Lx, Ly, Lz, Es, DT, time
CLOSE(11)	

plotnum=1
numplots=1+Nt/plotgap
DO n=1,numplots
OutputFileName = 'data/u'
ind = index(OutputFileName,' ') - 1
WRITE(number_file,'(i0)') 10000000+plotnum
OutputFileName = OutputFileName(1:ind)//number_file
ind = index(OutputFileName,' ') - 1
OutputFileName = OutputFileName(1:ind)//'.bov'
OutputFileName2='u'
ind = index(OutputFileName2,' ') - 1
OutputFileName2 = OutputFileName2(1:ind)//number_file
ind = index(OutputFileName2,' ') - 1
OutputFileName2 = OutputFileName2(1:ind)//'.datbin'
OPEN(unit=11,FILE=trim(OutputFileName),status="UNKNOWN")
REWIND(11)
WRITE(11,*) 'TIME: ',time
WRITE(11,*) 'DATA_FILE: ',trim(OutputFileName2)
WRITE(11,*) 'DATA_SIZE: ', Nx, Ny, Nz
WRITE(11,*) 'DATA_FORMAT: DOUBLE'
WRITE(11,*) 'VARIABLE: u'
WRITE(11,*) 'DATA_ENDIAN: LITTLE'
WRITE(11,*) 'CENTERING: ZONAL'
WRITE(11,*) 'BRICK_ORIGIN:', -Nx/2, -Ny/2, -Nz/2
WRITE(11,*) 'BRICK_SIZE:', Nx, Ny, Nz
WRITE(11,*) 'DIVIDE_BRICK: true'
WRITE(11,*) 'DATA_BRICKLETS:', Nx/2, Ny/2, Nz/2
CLOSE(11)

time=time+plotgap*DT
plotnum=plotnum+1
END DO
END PROGRAM BovCreate

( X)

用于在求解三维 Klein-Gordon 方程后创建 XDMF 头文件以帮助 Paraview 和 VisIt 读取数据的 Fortran 子程序代码下载

PROGRAM XdmfCreate
        !--------------------------------------------------------------------------------
        ! .. Purpose ..
        ! XdmfCreate is a postprocessing program which creates header files for
        ! Paraview and Visit
        ! It uses the INPUTFILE and assumes that the filenames in the program
        ! are
        ! consistent with those in the current file.
        !
        ! .. PARAMETERS .. INITIALIZED IN INPUTFILE
        ! time = start time of the simulation
        ! Nx = power of two, number of modes in the x direction
        ! Ny = power of two, number of modes in the y direction
        ! Nz = power of two, number of modes in the z direction
        ! Nt = the number of timesteps
        ! plotgap = the number of timesteps to take before plotting
        ! Lx = definition of the periodic domain of computation in x direction
        ! Ly = definition of the periodic domain of computation in y direction
        ! Lz = definition of the periodic domain of computation in z direction
        ! Dt = the time step
        ! Es = focusing or defocusing
        !
        ! REFERENCES
        !
        ! www.xdmf.org
        ! Paraview users mailing list
        ! VisIt users mailing list
        !
        ! ACCURACY
        !
        ! ERROR INDICATORS AND WARNINGS
        !
        ! FURTHER COMMENTS
        !------------------------------------------------------------------------------------
        ! EXTERNAL ROUTINES REQUIRED
        IMPLICIT NONE
        ! .. Scalar Arguments ..
        INTEGER(KIND=4) :: Nx, Ny, Nz, Nt, plotgap
        REAL(KIND=8)    :: Lx, Ly, Lz, DT, Es, time
        ! .. Local scalars ..
        INTEGER(KIND=4) :: stat,plotnum,ind,n,numplots
        REAL(KIND=8), PARAMETER :: pi=3.14159265358979323846264338327950288419716939937510d0
        ! .. Local Arrays ..
        CHARACTER*50    :: InputFileName, OutputFileName, OutputFileName2
        CHARACTER*10    :: number_file
        InputFileName='inputfile'
        OPEN(unit=11,FILE=trim(InputFileName),status="OLD")
        REWIND(11)
        READ(11,*) Nx, Ny, Nz, Nt, plotgap, Lx, Ly, Lz, Es, DT, time
        CLOSE(11)       
        time=0.0 
        plotnum=1
        numplots=1+Nt/plotgap
        OutputFileName='udata.xmf'
        OPEN(unit=11,FILE=trim(OutputFileName),status="UNKNOWN")
        REWIND(11)
        WRITE(11,'(a)') "<?xml version=""1.0"" ?>"
        WRITE(11,'(a)') "<!DOCTYPE Xdmf SYSTEM ""Xdmf.dtd"" []>"
        WRITE(11,'(a)') "<Xdmf xmlns:xi=""http://www.w3.org/2001/XInclude"" Version=""2.0"">"
        WRITE(11,'(a)') "<Domain>"
        WRITE(11,'(a)') " <Topology name=""topo"" TopologyType=""3DCoRectMesh"""
        WRITE(11,*) "  Dimensions=""",Nx,Ny,Nz,""">"
        WRITE(11,'(a)') " </Topology>"
        WRITE(11,'(a)') " <Geometry name=""geo"" Type=""ORIGIN_DXDYDZ"">"
        WRITE(11,'(a)') "  <!-- Origin -->"
        WRITE(11,'(a)') "  <DataItem Format=""XML"" Dimensions=""3"">"
        WRITE(11,*) -pi*Lx, -pi*Ly, -pi*Lz 
        WRITE(11,'(a)') "  </DataItem>"
        WRITE(11,'(a)') "  <!-- DxDyDz -->"
        WRITE(11,'(a)') "  <DataItem Format=""XML"" Dimensions=""3"">"
        WRITE(11,*) REAL(2.0*pi*Lx/Nx,kind(0d0)), REAL(2.0*pi*Ly/Ny,kind(0d0)), REAL(2.0*pi*Lz/Nz,kind(0d0))  
        WRITE(11,'(a)') "   </DataItem>"
        WRITE(11,'(a)') "  </Geometry>"
        WRITE(11,'(a)')
        WRITE(11,'(a)') "       <Grid Name=""TimeSeries"" GridType=""Collection"" CollectionType=""Temporal"">"
        WRITE(11,'(a)') "               <Time TimeType=""HyperSlab"">"
        WRITE(11,'(a)') "                       <DataItem Format=""XML"" NumberType=""Float"" Dimensions=""3"">"
        WRITE(11,*) "                   0.0", dt," 2" 
        WRITE(11,'(a)') "                       </DataItem>"
        WRITE(11,'(a)') "               </Time>"

        DO n=1,numplots

                OutputFileName = 'data/u'
                ind = index(OutputFileName,' ') - 1
                WRITE(number_file,'(i0)') 10000000+plotnum
                OutputFileName = OutputFileName(1:ind)//number_file
                ind = index(OutputFileName,' ') - 1
                OutputFileName = OutputFileName(1:ind)//'.datbin'
                
                OutputFileName2 = "<Grid Name=""T"
                WRITE(number_file,'(i0)') n 
                OutputFileName2 = OutputFileName2(1:13)//number_file
                ind =  index(number_file,' ') - 1
                OutputFileName2 = OutputFileName2(1:13+ind)//'" GridType="Uniform">'
                plotnum=plotnum+1
                
                WRITE(11,'(a)')
                WRITE(11,'(3X,a)') "   ",OutputFileName2
                WRITE(11,'(a)') "       <Topology Reference=""/Xdmf/Domain/Topology[1]""/>"
                WRITE(11,'(a)') "       <Geometry Reference=""/Xdmf/Domain/Geometry[1]""/>"
                WRITE(11,'(a)') "        <Attribute Name=""Displacement"" Center=""Node"">"
                WRITE(11,'(a)') "               <DataItem Format=""Binary"" "
                WRITE(11,'(a)') "                DataType=""Float"" Precision=""8"" Endian=""Little"" "
                WRITE(11,*) "            Dimensions=""",Nx, Ny, Nz,""">"
                WRITE(11,'(16X,a)') "           ",OutputFileName
                WRITE(11,'(a)') "               </DataItem>"
                WRITE(11,'(a)') "       </Attribute>"
                WRITE(11,'(3X,a)') "</Grid>"
        END DO
        WRITE(11,'(a)')"        </Grid>"
        WRITE(11,'(a)') "</Domain>"
        WRITE(11,'(a)') "</Xdmf>"
        CLOSE(11)
        END PROGRAM XdmfCreate

( Y)

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

#!/bin/bash
#PBS -N 3Dkg
#PBS -l nodes=1:ppn=1,walltime=04:15:00
#PBS -q flux
#PBS -l qos=muite_flux
#PBS -A muite_flux
#PBS -M [email protected]
#PBS -m abe
#PBS -V
#

# Go to location of executable
cd /nobackup/muite/KG3D

# The export command below assumes that the inputfile is in the same location as the executable.
# For best performance, the executable should be in the fast large file system. 
export inputfile=INPUTFILE

# make a directory where binary field data is output, the program may crash if this does not exist 
# or is not created 
mkdir data

# run the code

mpirun ./kg

练习

比较公式 2 和 3 中隐式和半隐式时间步进方案的精度。哪个方案在最短的实际时间内产生最精确的结果？
编写串行 Fortran 程序，使用公式 3 中的完全隐式时间步进方案来求解二维和三维 Klein-Gordon 方程。
编写 OpenMP 并行 Fortran 程序，使用公式 3 中的完全隐式时间步进方案来求解二维和三维 Klein-Gordon 方程。
MPI 命令 MPI_BCAST 用于列表 N 中的子程序 readinputfile。查找此命令（可能使用编程部分介绍中列出的参考资料之一），并解释它的功能。
编写 MPI 并行 Fortran 程序，使用公式 3 中的完全隐式时间步进方案来求解二维和三维 Klein-Gordon 方程。
比较具有周期性边界条件（ $\mathbb {T} ^{3}$ ）的完全三维模拟结果与 Donninger 和 Schlag^[9] 中总结的整个真实空间 ( $\mathbb {R} ^{3}$ ) 上爆发的分析预测。
Grenier^[10] 解释说线性 Klein-Gordon 方程可以写成两个耦合的薛定谔方程。可以将这种公式扩展到非线性 Klein-Gordon 方程。如果我们令

u=\phi +\xi \quad {\text{and}}\quad {\frac {\partial u}{\partial t}}=\phi -\xi

( 5)

那么这两个耦合方程

i{\frac {\partial }{\partial t}}{\begin{bmatrix}\phi \\\xi \end{bmatrix}}={\begin{bmatrix}-\Delta -1&-\Delta \\\Delta &\Delta +1\end{bmatrix}}{\begin{bmatrix}\phi \\\xi \end{bmatrix}}\pm {\begin{bmatrix}1\\-1\end{bmatrix}}{\frac {|\phi +\xi |^{2}(\phi +\xi )}{2}}

( 6)

等效于非线性 Klein-Gordon 方程

{\frac {\partial ^{2}u}{\partial t^{2}}}-\Delta u+u=\pm u^{3}

.

( 7)

1. 请补充细节以解释为什么公式 5 和 6 等价于公式 7 。具体来说，请展示通过将公式 5 和 7 中的两个方程相加减，我们可以得到 $i{\frac {\partial }{\partial t}}\left(\phi +\xi \right)=-\left(\phi -\xi \right)$ 和 $i{\frac {\partial }{\partial t}}\left(\phi -\xi \right)=-\Delta \left(\phi +\xi \right)-\left(\phi +\xi \right)\pm \left|\phi +\xi \right|^{2}\left(\phi +\xi \right)$ 。对第一个方程求导并代入第二个方程，然后回忆我们定义了 $u=\phi +\xi$ 在公式 5 中，我们得到了公式 7 中的克莱因-戈登方程。
2. 使用隐式中点法或克兰克-尼科尔森法求解这两个方程。

备注

↑ 关于此方程的一个不完整但易于理解的数学介绍可以在 http://wiki.math.toronto.edu/DispersiveWiki/index.php/Semilinear_NLW 中找到。
↑ Donninger 和 Schlag (2011)
↑ Yang (2006)
↑ Landau (1996)
↑ Grenier (1994)
↑ Donninger 和 Schlag (2011)
↑ Nakanishi 和 Schlag (2011)
↑ 目前，Matlab 的视频命令无法可靠地从非常长的模拟中生成单个视频，因此最好使用 Matlab 创建静止图像。
↑ Donninger 和 Schlag (2011)
↑ Grenier (1994)

参考文献

Donninger, R.; Schlag, W. (2011). "聚焦三次非线性克莱因-戈登方程的爆破/全局存在二分法的数值研究". Nonlinearity. 24: 2547–2562. doi:10.1088/0951-7715/24/9/009.

Grenier, W. (1994). 相对论量子力学. 施普林格出版社.

Landau, R.H. (1996). 量子力学 II. 约翰威立父子公司.

Schlag, W. (2011). 不变流形和色散哈密顿演化方程. 欧洲数学学会.

Yang, L. (2006), 克莱因-戈登-薛定谔方程的数值研究, 新加坡国立大学 {{citation}}: Cite has empty unknown parameter: |month= (help)

[1] 关于此方程的一个不完整但易于理解的数学介绍可以在 http://wiki.math.toronto.edu/DispersiveWiki/index.php/Semilinear_NLW 中找到。

[2] Donninger 和 Schlag (2011)

[3] Yang (2006)

[4] Landau (1996)

[5] Grenier (1994)

[6] Donninger 和 Schlag (2011)

[7] Nakanishi 和 Schlag (2011)

[8] 目前，Matlab 的视频命令无法可靠地从非常长的模拟中生成单个视频，因此最好使用 Matlab 创建静止图像。

[9] Donninger 和 Schlag (2011)

[10] Grenier (1994)

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]