Parallel Spectral Numerical Methods/The Klein-Gordon Equation

The Klein-Gordon Equation

Background

The focusing/defocusing nonlinear Klein-Gordon equation describes the evolution of a possible complex scalar field $u$ according to,

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

( 1)

where $+$ is the focusing case and $-$ the defocusing case in a similar manner to the nonlinear Schrödinger equation. Blow up of three dimensional radially symmetric real solutions to this equation have recently been numerically studied by Donninger and Schlag. Two dimensional simulations of the Klein-Gordon equation can be found in Yang. The linear Klein-Gordon equation occurs as a modification of the linear Schrödinger equation that is consistent with special relativity, see for example Landau or Grenier . At the present time, there have been no numerical studies of blow up of solutions to this equation without the assumption of radial symmetry. This equation has generated a large mathematical literature and is still poorly understood. Most of this mathematical literature has concentrated on analyzing the equation on an infinite three dimensional space with initial data that either decays exponentially as one tends to infinity or is nonzero on a finite set of the domain. Here, we will simulate this equation in a periodic setting. Since this equation is a wave equation, it has a finite speed of propagation of information, much as a sound wave in air takes time to move from one point to another. Consequently for short time simulations, a simulation of a solution that is only nonzero on a finite part of the domain is similar to a simulation on an infinite domain. However, over long times, the solution can spread out and interact with itself on a periodic domain, whereas on an infinite domain, the interaction over long times is significantly reduced and the solution primarily spreads out. Understanding the interactions in a periodic setting is an interesting mathematical problem. The Klein-Gordon equation has a conserved energy given by

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

The equation is also time reversible. For long time simulations, one wants to construct numerical methods that approximately conserve this energy and are also time reversible. When using Fourier spectral methods, we primarily need to ensure that the time discretization preserves these properties, since the spectral spatial discretization will typically automatically satisfy these properties. Following Donninger and Schlag, we use two schemes. First, an implicit-explicit time stepping scheme which is time reversible but only conserves the energy approximately and is given by

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

and second, a fully implicit time stepping scheme with fixed point iteration

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

which conserves a discrete energy exactly

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

As before, the superscript $n$ denotes the time step and $k$ denotes the iterate in the fixed point iteration scheme. Iterations are stopped once the difference between two successive iterates falls below a certain tolerance.

Matlab and Python Programs

Listings A , B , C and D demonstrate Matlab implementations of these time stepping schemes. In one dimension, the Klein-Gordon equation has easily computable exact solutions, (see for example Nakanishi and Schlag) which can be used to test the accuracy of the numerical schemes. These equations seem to display three possibilities for the behavior of solutions which are dependent on the initial conditions:

• the solutions could disperse or thermalize, that is a given localized initial condition spreads out over the entire space
• the solutions blow up or become infinite
• a portion of the solution travels around as a localized particle while the rest of the solution disperses.

Since the equations are reversible, there is also the possibility that a solution which is initially distributed over the spatial domain localizes itself.

( A)

A Matlab program to solve the 1-dimensional Klein Gordon equation 1 using the time discretization in eq. 2 Code download

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

A Python program to solve the 1-dimensional Klein Gordon equation 1 using the time discretization in eq. 2 . Code download

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

A Python program to solve the 1-dimensional Klein Gordon equation 1 using the time discretization in eq. 3 . Code download

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

A Matlab program to solve the 2-dimensional Klein Gordon equation 1 using the time discretization in eq. 3 Code download

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

A Python program to solve the 2-dimensional Klein Gordon equation 1 using the time discretization in eq. 3 . Code download

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

A Python program to solve the 2-dimensional Klein Gordon equation 1 using the time discretization in eq. 2 . Code download

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

A Matlab program to solve the 3-dimensional Klein Gordon equation 1 using the time discretization in eq. 2 Code download

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

A Python program to solve the 3-dimensional Klein Gordon equation, eq. 1 , using the time discretization in eq. 2 Code download

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

A Two-Dimensional OpenMP Fortran Program

The programs that we have developed in Fortran have become rather long. Here we add subroutines to make the programs shorter and easier to maintain. Listing E is the main Fortran program which uses OpenMP to solve the 2D Klein-Gordon equation. Notice that by using subroutines, we have made the main program significantly shorter and easier to read. It is still not as simple to read as the Matlab program, but is significantly better than some of the previous Fortran programs. It is also much easier to maintain, and once the subroutines have been written and debugged, they may be reused in other programs. The only drawback in using too many subroutines is that one may encounter a slight decrease in performance due to the overhead of calling a subroutine and passing data to it. The subroutines are in listings F , G , H , I , J , K and an example makefile is in listing L . Finally listing M contains a Matlab program which produces pictures from the binary files that have been computed. One can then use another program to take the images and create a video.

( E)

A Fortran program to solve the 2D Klein-Gordon equation Code download

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

A Fortran subroutine to get the grid to solve the 2D Klein-Gordon equation on Code download

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)

A Fortran program to save a field from the solution of the 2D Klein-Gordon equation Code download

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)

A Fortran subroutine to update arrays when solving the 2D Klein-Gordon equation Code download

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)

A Fortran subroutine to calculate the energy when solving the 2D Klein-Gordon equation Code download

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)

A Fortran subroutine to save final results after solving the 2D Klein-Gordon equation Code download

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)

An example makefile for compiling the OpenMP program in listing E Code download

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

A Matlab program to plot the fields produced by listing E Code download

% 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

A Three-Dimensional MPI Fortran Program using 2DECOMP&FFT

For this section, we will use the library 2DECOMP&FFT available from http://www.2decomp.org/index.html. The website includes some examples which indicate how this library should be used, in particular the sample code at http://www.2decomp.org/case_study1.html is a very helpful indication of how one converts a code that uses FFTW to one that uses MPI and the aforementioned library.

We now give a program for the three-dimensional nonlinear Klein-Gordon equation. The program uses the same subroutine structure as the two-dimensional code. To make the program easy to reuse, the subroutine listed in listing U has been created to read an INPUTFILE which specifies the parameters to use for the program and so the program does not need to be recompiled every time it is run. To enable the program to scale better, the arrays which hold the Fourier frequencies and grid points have also been decomposed so that only the portions of the arrays used on each processor are created and stored on the processor. The program is listed in listing W , to use this program simply compile it using MPI wrappers and gfortran (other compilers such as Intel fortran, PGI fortran or IBM XL fortran should also work), and then run it in the directory from which the INPUTFILE and data are stored. An example submission script is in listing Y, this requires an environment variable named "inputfile" to be set with the location of the INPUTFILE that specifies problem parameters, such as grid size and number of timesteps.

A further addition is a short postprocessing program to create header files to allow one to use the bov (brick of values) format that allows one to use the parallel visualization software VisIt. The program VisIt can be downloaded from https://wci.llnl.gov/codes/visit/home.html. This program also run on laptops, desktops as well as parallel computer clusters. Documentation on using VisIt is available here https://wci.llnl.gov/codes/visit/manuals.html and here http://www.visitusers.org/index.php?title=Main_Page. A short video tutorial on how to use VisIt remotely is available at http://cac.engin.umich.edu/resources/software/visit.html. An alternative program to create images is Paraview which can be downloaded from http://www.paraview.org/. To use Paraview, xdmf format header files can be used to get Paraview to read in the binary data output files. A short video tutorial on rendering with Paraview can be found at http://www.youtube.com/watch?v=rqCXFr73sfk&list=PL4FFYwZhtT0BlDBtjsSkQ4LcM1v9K24No. A Fortran program to create these header files is in listing X.

( N)

A Fortran program to solve the 3D Klein-Gordon equation Code download

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

A Fortran subroutine to get the grid to solve the 3D Klein-Gordon equation on Code download

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)

A Fortran subroutine to get the initial data to solve the 3D Klein-Gordon equation for Code download

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)

A Fortran program to save a field from the solution of the 3D Klein-Gordon equation Code download

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)

A Fortran subroutine to update arrays when solving the 3D Klein-Gordon equation Code download

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)

A Fortran subroutine to calculate the energy when solving the 3D Klein-Gordon equation Code download

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)

A Fortran subroutine to save final results after solving the 3D Klein-Gordon equation Code download

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)

A Fortran subroutine to read in the parameters to use when solving the 3D Klein-Gordon equation Code download

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)

An example makefile for compiling the MPI program in listing N Code download

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

A Fortran subroutine to create BOV (Brick of Values) header files after solving the 3D Klein-Gordon equation Code download

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)

A Fortran subroutine to create XDMF header files after solving the 3D Klein-Gordon equation for helping Paraview and VisIt read in the data Code download

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)

An example submission script for use on Flux. Change {your_username} appropriately. Code download

#!/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 muite@umich.edu
#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

Exercises

1. Compare the accuracy of the implicit and semi-implicit time stepping schemes in eqs. 2 and 3 . Which scheme produces the most accurate results in the least amount of real time?
2. Write serial Fortran programs to solve the two- and three-dimensional Klein-Gordon equations using the fully implicit time stepping scheme in eq. 3 .
3. Write OpenMP parallel Fortran programs to solve the two- and three-dimensional Klein-Gordon equations using the fully implicit time stepping scheme in eq. 3 .
4. The MPI command MPI_BCAST is used in the subroutine readinputfile, listed in list N . Look up this command (possibly using one of the references listed in the introduction to programming section) and explain what it does.
5. Write an MPI parallel Fortran program to solve the two- and three-dimensional Klein-Gordon equations using the fully implicit time stepping scheme in eq. 3 .
6. Compare the results of fully three-dimensional simulations with periodic boundary conditions ($\mathbb {T} ^{3}$) with analytical predictions for blow up on the entire real space ($\mathbb {R} ^{3}$) summarized in Donninger and Schlag.
7. Grenier explains that the linear Klein-Gordon equation can be written as two coupled Schrödinger equations. One can extend this formulation to the nonlinear Klein-Gordon equation. If we let

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

( 5)

then the two coupled equations

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

are equivalent to the nonlinear Klein-Gordon equation

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

( 7)

1. 1. Fill in the details to explain why eqs. 5 and 6 are equivalent to eq. 7 . In particular show that by adding and subtracting the two equations in eqs. 5 and 7 , we get $i{\frac {\partial }{\partial t}}\left(\phi +\xi \right)=-\left(\phi -\xi \right)$ and $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)$. Differentiating the first of these equations and substituting it into the second, then recalling that we defined $u=\phi +\xi$ in eq. 5 gives us the Klein-Gordon equation in eq. 7 .
   2. Solve these two equations using either the implicit midpoint rule or the Crank Nicolson method.

References

Donninger, R.; Schlag, W. (2011). "Numerical study of the blowup/global existence dichotomy for the focusing cubic nonlinear Klein-Gordon equation". Nonlinearity. 24: 2547–2562. doi:10.1088/0951-7715/24/9/009.

Grenier, W. (1994). Relativistic Quantum Mechanics. Springer.

Landau, R.H. (1996). Quantum Mechanics II. Wiley.

Schlag, W. (2011). Invariant Manifolds and Dispersive Hamiltonian Evolution Equations. European Mathematical Society.

Yang, L. (2006), Numerical studies of the Klein-Gordon-Schrodinger equations, National University of Singapore {{citation}}: Cite has empty unknown parameter: |month= (help)