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

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
```

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

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

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

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

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.

[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.

[1]

[2]

[3]

[4]