Docstoc

Tutorial introduction to Fortran 90

Document Sample
Tutorial introduction to Fortran 90 Powered By Docstoc
					Tutorial introduction to Fortran 90
A simple program
The study of physics usually begins with mechanics, so it makes sense to start our study of computational physics in the same place. The intent of this tutorial is to allow you to get some familiarity with the programming process — creating a code file, compiling and debugging a program — and not to teach you to program. The program examples can be typed in as is, and along the way you will get some idea of the language. Of course, there are other good reasons. Mechanics has evolved methods of extreme sophistication in the 350 years since Newton. One of the main impetuses has been the study of orbital motion in the solar system. But these sophisticated methods are there to provide a good starting point for perturbation treatments; often obscured in the beauty of the formalism is the fact that there are few exact solutions. So numerical treatments — starting already with Euler — have provided a powerful alternative. Let’s start off with something very simple: dropping a ball from, say, the leaning tower of Pisa. The equation of motion is m d2 y = −mg dt2

with y the height, t the time, m the mass of the ball and g the acceleration due to gravity. Of course, this has the solution at time t, assuming an initial height y0 and velocity zero at time t = 0, 1 y = y 0 − g t2 2 We can begin our programming in the simplest way, using the computer to evaluate this mathematical expression. We will evaluate the height for a set of times 0.1, 0.2, . . . , 1.0 s. We will do this pretty much in the same way we would do the task with a scientific calculator, so think of the program as a list of instructions for using such a tool. Here is a simple program in Fortran to calculate the height as a function of time:
program drop0 ! drop0.f90 - rwf ! Simple Fortran program ! straightforward evaluation of expressions implicit none

2

4

6

1

8

10

12

14

real : : real : : real : : real : : real : : real : : integer real : :

g =9.8 ! m =1.0 ! y0 = 1 0 . 0 ! v0 = 0 . 0 ! tmax = 1 . 0 ! dt = 0 . 1 ! : : Nsteps , step y

accel due to gravity (m/s^2) mass of ball ( kg) initial height 10 m initial vertical velocity component (m/s) drop time (s) time step

16

18

20

22

24

! determine number of time steps: Nsteps=tmax / dt+1 ! calculation loop do step =1, Nsteps t = step ∗ dt y = y0 −0.5∗ g∗t ∗∗2 print ∗ , t , y end do end program drop0

26

(Note that the line numbers printed here are just for reference below, and are not part of the program!) • Comments, which are there to provide explanations of what is going on, start from an exclamation mark and continue to the end of the line. • Note that the program is “bracketed” by program ... and end program ... statements. The name given in the statement is arbitrary, and does not have to be the same as the file. • implicit none is a statement that turns off some features of older versions of Fortran. • Note that all variables used in the program are declared in the first section of the program; some are defined to have initial values. • Lines 8-15 define various constants and variables. Numerical values could simply be used instead in the appropriate places, but it is usually better to give names to these and not have “magic numbers” floating about in your program. • This program uses a do loop to iterate over a range of integers which are used to calculate the time and hence the height of the ball. • In line 18, the number of steps to be taken is calculated. • Line 20 introduces the do loop: The variable step is given values in succession from 1 to Nsteps and the indented portion below is executed. Note that the statement ends with a statement end do. Do loops can be contained in other do loops, but cannot overlap. 2

• Lines 21-23 constitute the body of the loop. Note that the indentation is not required, but is useful for seeing the structure of the code. • Note the exponentiation operator is **. • The print statement will print the values of the variables in its argument list on the terminal. The asterisk * is a placeholder for a format statement, and implies a free format. You can specify the format of the data by using a proper format statement instead of the asterisk; this allows you to specify number of significant digits, etc. Note the comma after the *. Type this into a file. The comments are unimportant from the compiler’s point of view and can be omitted — they are there for you to keep track of things. Save the file with a .f90 extension, e.g. something like drop0.f90. Now, before we can run the program, it must be compiled. We will use the program ifc (Intel Fortran compiler) for this. The command for this is ifc -o drop0 drop0.f90 This takes the file drop0.f90 and creates an executable program drop0 (that’s the -o drop0 If the command prompt returns after this with no messages all went well. On the other hand it is probable that you got some stream of error messages and warnings. You will have to see what they say and try and fix your code, maybe with some help. OK, it compiles cleanly. Now run it: ./drop0 This line runs the program drop0 that exists in the current directory (that’s what the ./ indicates). (Technical note: the Unix PATH, which lists where executable programs exist, usually excludes the current directory for security reasons, so you need to explicitly set the path). You program runs! You get a list of numbers: do they seem sensible? Now make the problem a little less trivial by including the drag force, so that the equation of motion is d2 y 1 m 2 = −mg + CρAv 2 dt 2 Here, C is a constant; ρ the density of air and A the cross-sectional area of the body. 3

This actually has the solution, y = y0 + where γ = 1 ln γ2g 1 + (γv0 )2 cosh γgt

CρAv 2 /2mg and y0 is the initial height.

This expression has a number of mathematical functions in it, including the square root. There are suitable versions of these supplied by the system. Here is a simple program in Fortran to evaluate this expression. You can enter this, but it’s quicker to copy your first program to a new file (cp drop0.f90 drop1.f90) and modify the new file. Use what you have.
program drop1 ! drop1.f90 - rwf ! Simple Fortran program ! straightforward evaluation of expressions implicit none real : : real : : real : : real : : real : : real : : real : : real : : real : : real : : integer real : : g =9.8 ! m =1.0 ! y0 = 1 0 . 0 ! v0 = 0 . 0 ! tmax = 2 . 0 ! dt = 0 . 1 ! C =0.5 ! rho = 1 . 0 ! R =0.1 ! pi =3.1415926 : : Nsteps , step y , t , gamma accel due to gravity (m/s^2) mass of ball ( kg) initial height 10 m initial vertical velocity component (m/s) drop time (s) time step drag coefficient density of air ( kg/m^3) radius of ball (m)

2

4

6

8

10

12

14

16

18

20

22

24

26

28

! determine number of time steps: Nsteps=tmax / dt+1 ! calculation loop do step =1, Nsteps t = step ∗ dt gamma=sqrt ( C∗ rho ∗ pi∗R ∗ ∗ 2 / ( 2 ∗ m∗g ) ) y = y0 + ( 1 . 0 / ( g∗ gamma ∗ ∗ 2 ) ) ∗ log ( sqrt ( 1 . 0 + ( gamma ∗v0 ) ∗ ∗ 2 ) / ( cosh ( gamma ∗g∗t ) ) ) print ∗ , t , y end do end program drop1

30

• In many languages you will need to import the math functions from a library or module. Fortran has them built into the language. (You can of course import other things from a library or module). 4

• Note the way of calling a function — this is very similar to obvious mathematical notation. A function typically returns a value. • Note that trigonometric functions expect their argument in radians. • The system does not supply a value for π. The distance fallen in a given time is presumably smaller than in the original code, where there was no air resistance. How much smaller? Calculate it! Let’s add a line to calculate the vacuum value of free fall, say yvac = y0-0.5*g**2. Add this line to the do loop, and alter your program to print it out. Note that you have a new variable name. How about a graph? The simplest way to do this is to use a program like gnuplot. This requires the data to be plotted to be in a file. The simplest way to do this is to cut and paste from the screen to a text editor, and save the file, for example as drop1.out. The start gnuplot and plot it: gnuplot> plot ’drop1.out’ Hmm. This only plots one set of data. You have to tell the program which columns to plot: gnuplot> plot ’drop1.out’ using 1:2, ’drop1.out’ using 1:3 This can be abbreviated; you might like to use lines instead of points. gnuplot> plot ’drop1.out’ u 1:2 with lines, ’drop1.out’ u 1:3 w l Now let’s change the program so that it runs until one of the objects hits the ground at y = 0. We can do this by replacing the do loop with a do-while loop:
do while ( yvac > 0 . 0 ) step = step + 1

2

While we’re at it change the distance fallen to 50 m. Now, when we run the program, the data scrolls out of view. It’s time to write it to a file from within the program. We can do so by first opening a file
open(1, file=’drop2.out’)

5

Note that the file name is a character string enclosed in quotes. the “1” is a logical unit used to refer to the file. To write to the file we can replace the print statements with:
write(1,*) t,y,yvac

Here the “1” stands for the file; the “*” is still a placeholder for a format. Note that there is no comma before the t. Finally, we can close the file (after we have written everything to it.
close(1)

Plot the results using gnuplot.

6


				
DOCUMENT INFO
Shared By:
Stats:
views:87
posted:12/19/2009
language:English
pages:6
Description: Tutorial introduction to Fortran 90