Parallel Spectral Numerical Methods/The Klein-Gordon Equation
The Klein-Gordon Equation
Background
[1]The focusing/defocusing nonlinear Klein-Gordon equation describes the evolution of a possible complex scalar field 
according to,
-
,(
,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[2]. Two dimensional simulations of the Klein-Gordon equation can be found in Yang[3]. 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[4] or Grenier [5]. 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
.
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[6], 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
-
(
and second, a fully implicit time stepping scheme with fixed point iteration
-
(
which conserves a discrete energy exactly
-
.(
.As before, the superscript 
denotes the time step and 
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[7]) 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 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
|
( |
""" 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()
|
( |
""" 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()
|
( |
% 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
|
( |
#!/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()
|
( |
""" 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()
|
( |
% 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
|
( |
#!/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[8].
|
( |
!-------------------------------------------------------------------- ! ! ! 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
|
( |
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
|
( |
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
|
( |
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
|
( |
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
|
( |
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
|
( |
#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
|
( |
% 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
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. 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 is listed in listing W , to use this program simply compile it using gfortran, no special flags are required, and then run it in the directory from which the INPUTFILE and data are stored. 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 . To use Paraview, there is a better format of header file to use. A Fortran program to create these header files is in listing X.
|
( |
!-------------------------------------------------------------------- ! ! ! 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
|
( |
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
|
( |
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
|
( |
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
|
( |
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
|
( |
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
|
( |
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
|
( |
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
|
( |
# 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
|
( |
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
|
( |
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
Exercises
- 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?
- Write serial Fortran programs to solve the two- and three-dimensional Klein-Gordon equations using the fully implicit time stepping scheme in eq. 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 .
- 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.
- 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 .
- Compare the results of fully three-dimensional simulations with periodic boundary conditions (
) with analytical predictions for blow up on the entire real space (
) summarized in Donninger and Schlag[9]. - Grenier[10] 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
) with analytical predictions for blow up on the entire real space (
) summarized in Donninger and Schlag-
(
then the two coupled equations
-
(
are equivalent to the nonlinear Klein-Gordon equation
-
.(
.-
- 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
and 
. Differentiating the first of these equations and substituting it into the second, then recalling that we defined 
in eq. 5 gives us the Klein-Gordon equation in eq. 7 . - Solve these two equations using either the implicit midpoint rule or the Crank Nicolson method.
- 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
and 
. Differentiating the first of these equations and substituting it into the second, then recalling that we defined 
in eq. Notes
- ↑ An incomplete but easily accessible mathematical introduction to this equation can be found at http://wiki.math.toronto.edu/DispersiveWiki/index.php/Semilinear_NLW.
- ↑ Donninger and Schlag (2011)
- ↑ Yang (2006)
- ↑ Landau (1996)
- ↑ Grenier (1994)
- ↑ Donninger and Schlag (2011)
- ↑ Nakanishi and Schlag (2011)
- ↑ At the present time, Matlab’s video commands cannot reliably produce a single video from a very long simulation, so it is better to use Matlab to create still images.
- ↑ Donninger and Schlag (2011)
- ↑ Grenier (1994)
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.
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, http://scholarbank.nus.edu.sg/handle/10635/15515