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

From Wikibooks, open books for an open world
Jump to: navigation, search

Fortran Programs and Getting Started on Windows[edit]

Example Programs[edit]

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 Cohen[1]. We recognize that most people will be unfamiliar with Fortran and probably more familiar with Matlab[2], 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 Wagenbreth[3]. 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.

  1. 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
    
  2. 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
    


  3. 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}
    
  4. A Matlab plotting script[4] 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');
    
                            ( 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[edit]

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[edit]

  1. Please read the resources on the web page http://cac.engin.umich.edu/started/index.html to learn how to use the Flux cluster.
  2. 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.
  3. 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[edit]

  1. Metcalf, Reid and Cohen (2009)
  2. Although Matlab is written in C, it was originally written in Fortran and so has a similar style to Fortran.
  3. Levesque and Wagenbreth (2009)
  4. For 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.

References[edit]

Levesque, J.; Wagenbreth, G. (2011). High Performance Computing: Programming and Applications. CRC Press. 

Metcalf, M.; Reid, J.; Cohen, M. (2011). Modern Fortran Explained. Oxford University Press.