down

Section 7: Work Integrals in the Sun and Stirling engines
Project 1 Part 2

In the second part of Project 1 the main task is to write a Fortran subroutine for reading tables of data from text files, a Fortran module that calculates a work integral, given a set of values of temperatures and densities (to be read with the first subroutine), and then to use these to compute work integrals in the context of the Sun and in the context of Stirling engines.

You will need a correctly functioning splines.f90, with spline derivative and integral, from last weeks Section 6 exercise, so please make sure that you have completed that exercise.

In addition to the material from Section 6, a number of other files that are needed for Project 1 will appear in your 7_Stirling directory when you run the "cvs update -AdP" in your ComputerPhysics directory. You may assume that the pieces you get have been tested and are working correctly, but it is your responsibility to make sure, by testing, that the pieces you add or modify are working as intended.

Your Project1 directory should not readable by others and, in order to make the test realistic, you should attempt to complete the project on your own (no groups).

Here are your bonus points ...   Credits: 9/9

Note: More informtion about the Stirling engine kit demonstrated at the May 12 lecture is availble on an Absalon page, under "Course material".

Subsections
  approx time  
Getting the exercise files
1 min
Table reading Fortran subroutine
45 min
Makefile
15 min
Plots
30 min
The work integral Fortran module
45 min
Tests, results
45 min
Stirling engine optimization
1 hour
Home Work
2 hours


down up Getting the exercise files

[about 1 minute]

To extract additional files that you will need for Project 1, do

    cd ~/ComputerPhysics
    cvs update -AdP

In case of problems, see the CVS update help page.

The files will be extracted to a 7_Stirling directory. Since you need these files in the Project1 directory, and you normally will not need to change these files, we make soft links so the files can be seen from both directories:

    cd Project1
    ln -s ../7_Stirling/* ./         ! with Windows/CygWin remove the '-s' (to avoid editor problems)
    ls -l
    ls -lL
Note the output from the two last commands; the ls -l shows that some of the entries are soft links (marked with -> ). The ls -lL "resolves" the soft links, and shows the files pointed to.

Finally, let's make sure the Project1 directory is not readable by others:

    chmod -R o-r ~/ComputerPhysics/Project1


down up A Fortran subroutine for reading tables

[about 45 minutes]

The first task is to write a general Fortran subroutine read_table, in a file read_table.f90, which may be used for reading a table of values from a text file. The intended use is shown in the main.f90 and output.f90 files, which contain code such as:

  integer, parameter:: mcol=2, mrow=100
  integer:: nrow
  real, dimension(mcol,mrow):: table
  ...
  call read_table ('sun.dat', table, mcol, mrow, nrow)
The subroutine should open the file as a text file, read mcol values per line into each row of the array table, up to a maximum of mrow rows. The actual number of rows read should be returned in the nrow argument.

NOTE: The subroutine should not be inside a module, and hence the size of the table array needs to be specified as (mcol,mrow).

Pieces of code that may be used for reading the rows of data, while checking for the error that occurs at the end of the file, are shown in the PDF lecture notes.

Here is a hint about how to declare the file name in the subroutine (this information may, of course, also be found in your Fortran book)   Credits: -3/-3

When you are done with the read_table.f90 file, please confirm here:

Here is my working version of the read_table.f90 routine Locate and upload your read_table.f90 file: Credits: 10/-5
OK

The next step is to work on the Makefile, so the read_table subroutine can be tested and used:


down top Makefile

[about 15 minutes]

As the first step in modifying the Makefile add the lines in Makefile that are needed to make the output.x program. Those lines are analogous to the corresponding lines for the integral_test.x program, so use analogy reasoning to work it out.

Use a macro OUTPUT (or WORK as per an earlier version of these instructions -- the actual name doesn't matter) to list the object files needed for the output.x program :

OUTPUT = ....
... the rest of the Makefile is here ...

# by convention, put these dependency lines at the bottom of the Makefile:
output.o    : eqn_of_state_sun.o .......
....        : ...

When you are done, try it out. The Makefile should allow you to compile the program output.x with

    make clean					
    make output.x
and then, via these lines (where sun.dat is there because it is being read by output.x):
plot.dat: sun.dat output.x
        ./output.x
it should run the program when you give the command
    make plot.dat

Now put plot.dat (and only plot.dat for now!) as the default: target, so everything is done by the single command make. Try after a make clean:

    make clean
    make
NOTE: The test of your uploaded Makefile expects that the only other files that need to be present are the ones involved in making plot.dat -- that's why no other targets should be on the default line at this point.

Here is my working version of the Makefile Locate and upload your Makefile file: Credits: 10/-5
OK


down top Plotting

[about 30 minutes]

The output.x program writes out data in a plot.dat file, intended for making plots. The next task is to write an IDL procedure that can make plots from the plot.dat file.

figure1.pro

The first figure should correspond to the one that was shown during the lecture; a plot illustrating the path integral in the specific-volume / pressure plane.

The procedure should be called figure1, and should be in a file figure1.pro, which should contain

PRO figure1
  openr,1,'plot.dat'
  ... [ commands for reading the data ]
  close,1
  V = 1./rho
  plot, V, P, xticklen=1, yticklen=1, ... [ other options ]
END
To figure out what the commands for reading the data should be, look in the lecture notes.

Try it out by just giving the command

IDL> figure1
at an IDLDE or IDL command prompt. If you edit the file (for example for adding axis text options) you need to recompile before the change becomes effective:
IDL> .com figure1            [ or push the compile button in IDLDE ]
IDL> figure1                 [ or push the run button in IDLDE ]

Now add a call to the polyfill procedure for filling an area, after the plot command:

  plot, V, P, ...
  polyfill, V, P

By adding an option color=100 to the polyfill command you can make the area gray, and by loading a different color table (with xloadct or for example loadct,39) you can fill the area with color (use the online help if needed).

polyfill overwrites the grid. To restore it, repeat the plot command a second time, now with a , /noerase option:

plot, ....
polyfill, ...
plot, ...., /noerase
Control the number of squares shown by adding xticks and yticks options to the plot commands (look up these keywords in the on-line help, or just try them out):
plot, ...., xticks=8, yticks=8
polyfill, ...
plot, ...., xticks=8, yticks=8, /noerase

Finally, replace the , xticks=.., yticks=.. string in both plot commands with , _extra=extra, and add the same thing as an argument to the first line:

PRO figure1, _extra=extra
...
plot, ...., _extra=extra
polyfill, ...
plot, ...., _extra=extra, /noerase
END
The _extra parameter makes it possible to pass any number of extra parametes on to procedures that are called inside other procedures.

Use it to experiment with different plot options directly on the command line:

IDL> figure1,xticks=9,yticks=8
IDL> figure1,xticks=7,yticks=8
IDL> figure1,xticks=7,yticks=8,title='figure1',xtitle='V',ytitle='P'
IDL> figure1,xticks=7,yticks=8,title='figure1',charsize=2    ; makes it easier to read the numbers!
IDL> figure1,xticks=7,yticks=8,title='figure1',line=1
IDL> figure1,xticks=7,yticks=8,title='figure1',thick=3

Here is my working version of the figure1.pro routine Locate and upload your figure1.pro file: Credits: 10/-5
OK

You may use the method from the lecture (estimating the surface of the colored area on the plot) to give an estimate of the surface (within 50%) in the field below.

Give an estimate of the surface (within 50%): Credits: 5/-2
OK



down top The work integral Fortran module

[about 45 minutes]

The next task is writing a Fortran module for the work integral (see the lecture notes for an explanation of work integrals).

Creating the file

The file should be called work.f90 and the module should be called work.

If you are uncertain about the structure of the file you may choose to buy this hint   Credits: -3/-3

or just look in any other Fortran module file, and reason by analogy.

To figure out the argument list of the work_integral function, look at how it is called in the main.f90 program, and look at how arguments are handled in your own spline_integral function.

Local declarations

In addition to the argument arrays (T and P) you may need to declare a local scratch array (or you may choose to write the code such that it does not need a scratch array). In any case you need an array declaration for the return value. In case of doubt, look in your spline module how to do this.

The executable code

The executable code should produce the return value; i.e., an array with the result of applying spline_integral to the integrand.

What about the integrand? It's definition is in the lecture notes, and you have access to the spline_derivative and pressure functions.

Makefile

The Makefile for the full Work.x program should contain

FC = gfortran
...
WORK =  ... [ object files with correct dependencies at the bottom of the Makefile ] ... main.o
work.x: $(WORK)
        $(FC) $(FFLAGS) $(WORK) -o work.x
work:
        ./work.x
[ remember that there must be a TAB, not 8 blanks, in the rule lines ]
Add such a block to the Makefile and add work to the default: target in the Makefile as a dependency (to force the work.x program to be executed). After this you can make the work.x program and attempt to run it, by just typing
   make

Concerning the dependencies of the object files, consider the following:

Think about what this means for the order of the .o files. Also, consider which .o files should be in the list.

Here is my updated version of Makefile Locate and upload your Makefile file: Credits: 10/-5
OK



down top General path integrals, test, result

[about 45 minutes]

We need to make sure that the work integral routine gives the correct result. As constructed, this is a bit difficult, because the module is too specialized. By changing it to a module that can do general path integrals one can apply it, as a special case, to a case with a known result. Let's do that:

path.f90

Make a copy of work.f90 as path.f90 and change the names inside from work_integral to path_integral and so on. Then make the function general, by moving the equation-of-state call out of the routine, and making path.f90 compute only the path integral of A dB, again rewritten as A(i) dB(i)/di di.

sun.f90

Since you moved the equation-of-state call out of path.f90 you need a slightly different main program, which should be called sun.f90. So, copy main.f90 to sun.f90, add a call to equation_of_state, and change the call to work_integral into a call to path_integral. So,
  cp main.f90 sun.f90

and update the program so it computes

  work(1:nrow) = path_integral (pressure(T(1:nrow),rho(1:nrow)), 1./rho(1:nrow))
  print *,'Result =',work(nrow)

circle.f90

Now it is easy to construct a test case with a known solution. Let A and B be defined as
  pi = 4.*atan(1.)
  A(:) = cos((/(i,i=0,mrow-1)/)*2.*pi/(mrow-1))
  B(:) = sin((/(i,i=0,mrow-1)/)*2.*pi/(mrow-1))
Then it is not difficult to figure out what the result should be. Try drawing the path integral with pen and paper (or write the numbers out in a file and make a new figure_test.pro to plot the enclosed surface with).

To make the test, copy sun.f90 to circle.f90, replace the reading of data with the loop above, and add the appropriate lines in your Makefile so that you can say

    make circle
to run the test.

When you have done this, compare the last value returned from the path_integral(A,B) with the expected result, and have the program print OK if the test succeeds.

Here is my working version of the path.f90 routine Locate and upload your path.f90 file: Credits: 10/-5
OK

The work integral result

Once you are certain that the path module is correct, use the sun.f90 program, compiled into a sun.x executable file (add another section to the Makefile, in the same way as before, to make sun.x), and use it to compute (what should be -- check it!) the same result as from the work.x program.

When done, confirm below. Please note that you must get it right with at most one retry, so use copy & paste, carefully, from the output of the work.x or path.x programs -- and then remove any exponent notation and move the decimal point instead (since the web form input unfortunately cannot use exponent notation).

There should be no chance of making a mistake here, since you already have a ball-park estimate (from the plot area estimate), and you have tested the path integral module with the test.x program.

My result for the work integral is (within 0.1%): Credits: 50/-5
OK

Additional points in case you can determine the maximum efficiency (interval 0-1, within 0.001) of the ideal Carnot process that just encloses the process discussed in this exercise.... Credits: 20/-5
OK



down top Stirling engine optimization

[about 1 hour]

Now that you have a setup that can compute the work integral for one set of (solar) data, it is easy to make a similar setup that reads data for a Stirling engine instead. These are data for normal air (cf. the "coffee-cup Stirling engine" ;-), so you need to change to eqn_of_state_air.o in the Makefile setup for a new target:

STIRLING = eqn_of_state_air.o splines.o path.o stirling.o
stirling.x: $(STIRLING)
        $(FC) $(FFLAGS) $(STIRLING) -o stirling.x $(LFLAGS)
stirling: stirling.x stirling.dat
        ./stirling.x
Create the file stirling.f90 as a copy of sun.f90, initially with the only difference that it reads stirling.dat instead of sun.dat. Run it with "make stirling" and enter the result here (1% accuracy):

The work integral for stirling.dat is: Credits: 10/-5
OK

Optimizing the Stirling engine

You can get extra points by improving the efficiency of the Stirling engine: Add a loop with an index that goes from -5 to +5 and shift a copy of the temperature data that many indices (use the built-in Fortran function CSHIFT(ARRAY,I)). Print the work integral for each case. For what shift does the work integral reach a maximum?

NOTE: ARRAY in this case is not the whole (1:mrow) array, of course, but also not even the whole (1:nrow) array. You have essentially two alternatives to get it right; 1) be very smart and think very carefylly, or 2) "cheat"; -- print the values out before and after the shift and look at what they are (and think about what they should be).

The work integral is largest for I = Credits: 10/-5
OK

What is the value of the work integral then (0.1% accuracy)?

The maximized work integral for stirling.dat is: Credits: 10/-5
OK



down top Home Work

[about 2 hours]

If you did not finish the project during this exercise you may want to look into how to use CVS to access you project (repository) from home. It is worth the effort, since you can use this technique for a lot of things in the future.

Otherwise, spend the hours at home studying the examples in your Fortran 90/95 book and try some of the examples from the book (most people have part of their memory "in the fingers"; i.e., typing something in from a book is, rather than just a boring exercise, a surprisingly efficient way of learning syntax).


$Id: index.php,v 1.13 2010/05/12 06:05:21 aake Exp $ aake@astro.ku.dk