Parallel Spectral Numerical Methods/Fortran Programs and Getting Started on Windows

Fortran Programs and Getting Started on Windows Example Programs To do parallel programming using OpenMP or MPI (Message passing interface), we typically need to use a lower level language than Matlab such as Fortran. Another possible choice of language is C, however Fortran has superior array handling capabilities compared to C, and has a similar syntax to Matlab, so is typically easier to use for scientific computations which make heavy use of regular arrays. It is therefore useful to introduce a few simple programs in Fortran before we begin studying how to create parallel programs. A good recent reference on Fortran is Metcalf, Reid and CohenMetcalf, Reid and Cohen (2009). We recognize that most people will be unfamiliar with Fortran and probably more familiar with MatlabAlthough Matlab is written in C, it was originally written in Fortran and so has a similar style to Fortran., C or C++, but we expect that the example codes will make it easy for anyone with some introductory programming background. A recent guide which describes how to write efficient parallel Fortran code is Levesque and WagenbrethLevesque and Wagenbreth (2009). Our programs are written to be run on the Flux cluster at the University of Michigan. More information on this cluster can be found at http://cac.engin.umich.edu/resources/systems/flux/ and at http://cac.engin.umich.edu/started/index.html. Below are four files you will need to run this. <li>A makefile to compile the Fortran code on Flux in listing A . This should be saved as makefile. Before using the makefile to compile the code, you will need to type module load fftw/3.2.1-intel at the command line prompt once logged into Flux. Then place the makefile and heat.f90 in the same directory, the example files below assume this directory is {${HOME}/ParallelMethods/Heat} and type make to compile the file. Once the file is compiled type qsub fluxsubscript to get the cluster to run your program and then output the results. The programs that follow use the library FFTW to do the fast Fourier Transforms. More information on this library can be found at http://www.fftw.org/. A An example makefile for compiling a Fourier spectral Fortran heat equation program Code download #define the complier COMPILER = mpif90 # compilation settings, optimization, precision, parallelization FLAGS = -O0 # libraries LIBS = -L${FFTW_LINK} -lfftw3 -lm # source list for main program SOURCES = heat.f90 test: $(SOURCES) ${COMPILER} -o heat $(FLAGS) $(SOURCES) $(LIBS) clean: rm *.o <li>The Fortran program in listing B – this should be saved as heat.f90 B A Fortran Fourier spectral program to solve the heat equation using backward Euler timestepping Code download !-------------------------------------------------------------------- ! ! ! PURPOSE ! ! This program solves heat equation in 1 dimension ! u_t=\alpha*u_{xx} ! using a the backward Euler method for x\in[0,2\pi] ! ! The boundary conditions are u(0)=u(2\pi) ! The initial condition is u=sin(x) ! ! .. Parameters .. ! Nx = number of modes in x - power of 2 for FFT ! Nt = number of timesteps to take ! Tmax = maximum simulation time ! plotgap = number of timesteps between plots ! FFTW_IN_PLACE = value for FFTW input ! FFTW_MEASURE = value for FFTW input ! FFTW_EXHAUSTIVE = value for FFTW input ! FFTW_PATIENT = value for FFTW input ! FFTW_ESTIMATE = value for FFTW input ! FFTW_FORWARD = value for FFTW input ! FFTW_BACKWARD = value for FFTW input ! pi = 3.14159265358979323846264338327950288419716939937510d0 ! L = width of box ! alpha = heat conductivity ! .. Scalars .. ! i = loop counter in x direction ! n = loop counter for timesteps direction ! allocatestatus = error indicator during allocation ! start = variable to record start time of program ! finish = variable to record end time of program ! count_rate = variable for clock count rate ! planfx = Forward 1d fft plan in x ! planbx = Backward 1d fft plan in x ! dt = timestep ! .. Arrays .. ! u = approximate REAL solution ! v = Fourier transform of approximate solution ! vna = temporary field ! .. Vectors .. ! kx = fourier frequencies in x direction ! x = x locations ! time = times at which save data ! name_config = array to store filename for data to be saved ! ! REFERENCES ! ! ACKNOWLEDGEMENTS ! ! ACCURACY ! ! ERROR INDICATORS AND WARNINGS ! ! FURTHER COMMENTS ! Check that the initial iterate is consistent with the ! boundary conditions for the domain specified !-------------------------------------------------------------------- ! External routines required ! ! External libraries required ! FFTW3 -- Fast Fourier Transform in the West Library ! (http://www.fftw.org/) PROGRAM main ! Declare variables IMPLICIT NONE INTEGER(kind=4), PARAMETER :: Nx=64 INTEGER(kind=4), PARAMETER :: Nt=20 REAL(kind=8), PARAMETER & :: pi=3.14159265358979323846264338327950288419716939937510d0 REAL(kind=8), PARAMETER :: L=5.0d0 REAL(kind=8), PARAMETER :: alpha=0.50d0 REAL(kind=8) :: dt=0.2d0/REAL(Nt,kind(0d0)) COMPLEX(KIND=8), DIMENSION(:),ALLOCATABLE :: kx REAL(kind=8), DIMENSION(:),ALLOCATABLE :: x COMPLEX(KIND=8), DIMENSION(:,:),ALLOCATABLE :: u,v REAL(kind=8), DIMENSION(:),ALLOCATABLE :: time COMPLEX(KIND=8), DIMENSION(:),ALLOCATABLE :: vna INTEGER(kind=4) :: i,j,k,n INTEGER(kind=4) :: start, finish, count_rate, AllocateStatus INTEGER(kind=4), PARAMETER :: FFTW_IN_PLACE = 8, FFTW_MEASURE = 0, & FFTW_EXHAUSTIVE = 8, FFTW_PATIENT = 32, FFTW_ESTIMATE = 64 INTEGER(kind=4), PARAMETER :: FFTW_FORWARD = -1, FFTW_BACKWARD=1 COMPLEX(KIND=8), DIMENSION(:),ALLOCATABLE :: fftfx,fftbx INTEGER(kind=8) :: planfx,planbx CHARACTER*100 :: name_config CALL system_clock(start,count_rate) ALLOCATE(kx(1:Nx),x(1:Nx),u(1:Nx,1:1+Nt),v(1:Nx,1:1+Nt),& time(1:1+Nt),vna(1:Nx),fftfx(1:Nx),fftbx(1:Nx),& stat=AllocateStatus) IF (AllocateStatus .ne. 0) STOP ! set up ffts CALL dfftw_plan_dft_1d(planfx,Nx,fftfx(1:Nx),fftbx(1:Nx),& FFTW_FORWARD,FFTW_ESTIMATE) CALL dfftw_plan_dft_1d(planbx,Nx,fftbx(1:Nx),fftfx(1:Nx),& FFTW_BACKWARD,FFTW_ESTIMATE) PRINT *,'Setup FFTs' ! setup fourier frequencies DO i=1,1+Nx/2 kx(i)= cmplx(0.0d0,1.0d0)*REAL(i-1,kind(0d0))/L END DO kx(1+Nx/2)=0.00d0 DO i = 1,Nx/2 -1 kx(i+1+Nx/2)=-kx(1-i+Nx/2) END DO DO i=1,Nx x(i)=(-1.00d0 + 2.00d0*REAL(i-1,kind(0d0))/REAL(Nx,KIND(0d0)))*pi*L END DO PRINT *,'Setup grid and fourier frequencies and splitting coefficients' u(1:Nx,1)=sin(x(1:Nx)) ! transform initial data CALL dfftw_execute_dft_(planfx,u(1:Nx,1),v(1:Nx,1)) PRINT *,'Got initial data, starting timestepping' time(1)=0.0d0 vna(1:Nx)=v(1:Nx,1) PRINT *,'Starting timestepping' DO n=1,Nt DO i=1,Nx vna(i)=vna(i)/(1-dt*kx(i)*kx(i)) END DO PRINT *,'storing plot data ',n time(n+1)=time(n)+dt v(1:Nx,n+1)=vna(1:Nx) CALL dfftw_execute_dft_(planbx,v(1:Nx,n+1),u(1:Nx,n+1)) u(1:Nx,n+1)=u(1:Nx,n+1)/REAL(Nx,KIND(0d0)) ! normalize END DO PRINT *,'Finished time stepping' CALL system_clock(finish,count_rate) PRINT*,'Program took ',REAL(finish-start)/REAL(count_rate),'for execution' ! Write data out to disk name_config = 'u.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO j=1,1+Nt DO i=1,Nx WRITE(11,*) REAL(u(i,j)) END DO END DO CLOSE(11) name_config = 'tdata.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO j=1,1+Nt WRITE(11,*) time(j) END DO CLOSE(11) name_config = 'xcoord.dat' OPEN(unit=11,FILE=name_config,status="UNKNOWN") REWIND(11) DO i=1,Nx WRITE(11,*) x(i) END DO CLOSE(11) PRINT *,'Saved data' DEALLOCATE(kx,x,u,v,& time,vna,fftfx,fftbx,& stat=AllocateStatus) IF (AllocateStatus .ne. 0) STOP CALL dfftw_destroy_plan(planbx) CALL dfftw_destroy_plan(planfx) CALL dfftw_cleanup() PRINT *,'Program execution complete' END PROGRAM main <li> An example submission script to use on the cluster in Listing C – this should be saved as fluxsubscript. More examples can be found at http://cac.engin.umich.edu/resources/software/pbs.html. To use it, please change the email address from your_uniqname@umich.edu to an email address at which you can receive notifications of when jobs start and are finished. C An example submission script for use on Flux Code download #!/bin/bash #PBS -N heatequation #PBS -l nodes=1,walltime=00:10:00 #PBS -l qos=math471f11_flux #PBS -A math471f11_flux #PBS -q flux #PBS -M your_uniqname@umich.edu #PBS -m abe #PBS -V # Create a local directory to run and copy your files to local. # Let PBS handle your output mkdir /tmp/${PBS_JOBID} cp ${HOME}/ParallelMethods/Heat/heatequation /tmp/${PBS_JOBID}/heatequation cd /tmp/${PBS_JOBID} ./heatequation #Clean up your files cd cd ParallelMethods/Heat # Retrieve your output cp /tmp/${PBS_JOBID}/u.dat ${HOME}/ParallelMethods/Heat/u.dat cp /tmp/${PBS_JOBID}/xcoord.dat ${HOME}/ParallelMethods/Heat/xcoord.dat cp /tmp/${PBS_JOBID}/tdata.dat ${HOME}/ParallelMethods/Heat/tdata.dat /bin/rm -rf /tmp/${PBS_JOBID} A Matlab plotting scriptFor many computational problems, one can visualize the results with 10-100 times less computational power than was needed to generate the results, so for problems which are not too large, it is much easier to use a high level language like Matlab to post-process the data. to generate Fig. i is in listing D . D A Matlab program to plot the computed results Code download % A Matlab program to plot the computed results clear all; format compact, format short, set(0,'defaultaxesfontsize',18,'defaultaxeslinewidth',.9,... 'defaultlinelinewidth',3.5,'defaultpatchlinewidth',5.5); % Load data load('./u.dat'); load('./tdata.dat'); load('./xcoord.dat'); Tsteps = length(tdata); Nx = length(xcoord); Nt = length(tdata); u = reshape(u,Nx,Nt); % Plot data figure(3); clf; mesh(tdata,xcoord,u); xlabel t; ylabel x; zlabel('u'); center|thumb| i The solution to the heat equation computed by the Fortran program in B and post-processed by the Matlab program in D . Getting Started on Windows To run these Programs on Windows you need Cygwin installed. Cygwin emulates a bash shell on windows. You also need the following packages for Cygwin: make, mpi, fftw3 (list incomplete). Below is the adapted make file, the others remain the same. A makefile to compile the Fortran code in listing E . This should be saved as makefile instead of ( A ). Then place the makefile and heat.f90 ( B ) in the same directory. type make to compile the file. Once the file is compiled type ./heat.exe to get your computer to run the program. use the Matlab plotting script from above to get the picture. E An example makefile for compiling a Fourier spectral Fortran heat equation program Code download #define the complier COMPILER = mpif90 # compilation settings, optimization, precision, parallelization FLAGS = -O0 # libraries LIBS = -lfftw3 -lm # source list for main program SOURCES = heat.f90 test: $(SOURCES) ${COMPILER} -o heat $(FLAGS) $(SOURCES) $(LIBS) clean: rm *.o rm heat Exercises Please read the resources on the web page http://cac.engin.umich.edu/started/index.html to learn how to use the Flux cluster. Modify the Fortran program for the 1-D heat equation to solve the Allen-Cahn equation, with your choice of time stepping scheme. Create a plot of the output of your run. Include the source code and plot in your solutions. Modify the Fortran program for the 1-D heat equation to solve the 2-D heat equation with your choice of time stepping scheme. Your program should save the field at each time step rather than putting all the fields in a single large array. Create a plot of the initial and final states of your run. Include the source code and plots in your solutions. Notes References LevWag11 Levesque J. Wagenbreth G. High Performance Computing: Programming and Applications CRC Press 2011 MetReiCoh11 Metcalf M. Reid J. Cohen M. Modern Fortran Explained Oxford University Press 2011