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

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

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

    chmod -R o-r ~/ComputerPhysics/Project1

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)

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:

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.

• This makes sense to you?
  • Y -> good, continue
  • N -> back to Exercise 6, Makefiles...

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 temperature and density data are for the Sun, so you need the eqn_of_state_sun.f90 equation-of-state routine.
• To make the full compilation the files need to be compiled in the correct order. For instance the eqn_of_state_sun.o must come before files which USE this module, so that the module info (.mod file) has been made when the files that USE it are compiled.
• The compilor does not know this order so you have to provide this information in the Makefile.
• The most general way it to make a target for each object file (*.o) in the project and list the .o names of the files it depends on (only "use" dependencies are needed -- check which modules are being references with "use" statements, and where the module is defined):

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

plot.dat: sun.dat output.x
        ./output.x

    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

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

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

Try it out by just giving the command

IDL> figure1

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

• See IDL startup if colors don't work in IDL.

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

plot, ....
polyfill, ...
plot, ...., /noerase

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

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

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 ]

   make

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

• The module information file (which often has a .mod file name extension) is created at the same time as the .o object file.
• The module information is needed when compiling any subroutine that makes use of the module.

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

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

sun.f90

  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

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

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

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

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

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

Optimizing the Stirling engine

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 = -5 -4 -3 -2 -1 1 2 3 4 5 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

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