PROGRAM jacobi IMPLICIT NONE !HPF$ NOSEQUENCE ! All these are parallel variables with "field" and "mask" distributed ! blockwise. INTEGER, PARAMETER :: n = 200, iterations = 2000 REAL, DIMENSION(n), PARAMETER :: north_boundary = 30.0, & east_boundary = 80.0, west_boundary = 80.0, south_boundary = 250.0 REAL, DIMENSION(n, n) :: field = 30.0 LOGICAL, DIMENSION(n, n) :: mask = .TRUE. !HPF$ ALIGN mask(:,:) with field !HPF$ DISTRIBUTE (*, BLOCK) :: field ! The following are HDF related variables. INTERFACE EXTRINSIC(hpf_serial) FUNCTION d8aimg(file_name, raster, nx, ny, mode) INTEGER :: d8aimg INTEGER :: nx, ny, mode INTEGER(kind=1) :: raster(nx,ny) CHARACTER(len=*) :: file_name END FUNCTION d8aimg EXTRINSIC(hpf_serial) FUNCTION d8pimg(file_name, raster, nx, ny, mode) INTEGER :: d8pimg INTEGER :: nx, ny, mode INTEGER(kind=1) :: raster(nx,ny) CHARACTER(len=*) :: file_name END FUNCTION d8pimg END INTERFACE INTEGER(kind=1), DIMENSION(n, n) :: raster INTEGER :: i, j, status INTEGER, PARAMETER :: comp_rle = 11 field(UBOUND(field, dim=1), :) = east_boundary mask(UBOUND(mask, dim=1), :) = .FALSE. field(LBOUND(field, dim=1), :) = west_boundary mask(LBOUND(mask, dim=1), :) = .FALSE. field(:, UBOUND(field, dim=2)) = north_boundary mask(:, UBOUND(mask, dim=2)) = .FALSE. field(:, LBOUND(field, dim=2)) = south_boundary mask(:, LBOUND(mask, dim=2)) = .FALSE. raster=INT(field, 1) status=d8pimg('movie.hdf', raster, n, n, COMP_RLE) DO i = 1, iterations WHERE (mask) field = (EOSHIFT(field, 1, dim=1) + EOSHIFT(field, -1, dim=1) & + EOSHIFT(field, 1, dim=2) + EOSHIFT(field, -1, dim=2)) * 0.25 END WHERE IF (MOD(i, 10) .EQ. 0) THEN raster = INT(field, 1) status = d8aimg('movie.hdf', raster, n, n, comp_rle) WRITE(*,*) i, ' -> ', (INT(field(j, 10)), j = 1, 5), & (INT(field(j, 10)), j=96, 105), & (INT(field(j, 10)), j=196, 200) END IF END DO END PROGRAM jacobiHere is how to compile, link, and run this program:
gustav@sp20:../jacobi 11:14:45 !523 $ make jacobi-raster xlhpf90 -o jacobi-raster jacobi-raster.f \ -L/afs/ovpit.indiana.edu/@sys/HDF/lib -lmfhdf -ldf -ljpeg -lz -lm ** jacobi === End of Compilation 1 === 1501-510 Compilation successful for file jacobi-raster.f. gustav@sp20:../jacobi 11:14:51 !524 $ make run-raster poe jacobi-raster -euidevice css0 -procs 8 -euilib ip -ilevel 3 \ > jacobi-raster.log 2> jacobi-raster.err gustav@sp20:../jacobi 11:25:33 !525 $
/usr/lpp/pessl/example/hpfand which we've been working with. You must fix this bug, because otherwise you will not be able to write large images.
main.f you will find the following code:
do it = 1, dif_ntemps do i = 1, dif_npts iy = (i-1)/dif_nx+1 ix = mod(i-1,dif_nx)+1 do k = 1, dif_npts temp(k,it) = temp(k,it) + xsines(k,ix) * ysines(k,iy) & * x(i,it) enddo enddo enddotowards the end of the file. This is the part of the program where final temperatures are calculated for various times. The index
iin the second
DOloop runs from
dif_npts, which is incorrect. It should run from
n, because that index numbers harmonics, and we have harmonics altogether. The variable that is used in the code,
dif_npts, stores the numbers of points of the physical grid used, e.g., for writing images.
If your image is, say,
dif_npts will be 40,000. On the
other hand the default number of harmonics is 7 in each direction.
So n = 49, which is a lot less than 40,000. Trying to run that
second loop from 1 through 40,000 will invariably result in incorrect
In summary that part of the code should look as follows:
do it = 1, dif_ntemps do i = 1, n iy = (i-1)/dif_nx+1 ix = mod(i-1,dif_nx)+1 do k = 1, dif_npts temp(k,it) = temp(k,it) + xsines(k,ix) * ysines(k,iy) & * x(i,it) enddo enddo enddo
I'd like to thank Xin Hong, for having drawn my attention to this problem.
Code fragments presented in this document and in
are correct in this respect.
You should avoid pixel colours 0 and 1, because these often correspond to black and white, and you want a neatly graded display from red through dark blue, without any white and black spots or lines. Other low value colours may be bad too, e.g., they may or may not correspond to pure yellow, green, red, etc. You can find the lowest "safe" number simply by experimenting (or reading HDF manuals). In our jacobi.f example, I've used 30, which is safe enough.
The movie should look as follows. Initially there should be a reddish region somewhere inside the square dissolving gradually towards blue around the edges. Then as time goes on, the red should become gradually bluer until towards the end you'll have blue everywhere, i.e., all the heat should flow away from the beam.
Observe that in program jacobi.f I didn't have to do any scaling, because I have forced my temperatures to stay between 30 and 250 by the virtue of imposed boundary conditions. I could therefore map from temperature onto the raster image without any scaling at all. Of course I did it on purpose, to avoid cluttering the program with that kind of stuff.
But you will have to scale your data in the Spectral Method example.
/usr/lpp/pessl/doc/postscript/pessl.00.pson the SP. You can view it on-line with
/afs/ovpit.indiana.edu/@sys/X11R6.4/bin/ghostviewYou can also copy it to a more convenient location, e.g., to your AFS directory.
/usr/lpp/xlhpf/postscriptThere are three files there, which are all compressed. Use the command
uncompressto uncompress them. You will have to copy them to your AFS directory to do that. There may not be enough space in your SP directory. The files are:
/afs/ovpit.indiana.edu/common/HDF/doc/Users_Guide/ug41r2.pdf /afs/ovpit.indiana.edu/common/HDF/doc/Ref_Manual/RM_4.1r2b.pdfThese have to be read with Acrobat Reader,
Sometimes AFS doesn't come up properly on some of the nodes that you
have declared in your
host.list file. If POE hits any of those
nodes during execution, and the run has been started up from AFS,
it'll bump out and complain. As it happened, AFS didn't run on
nodes 40 and 42 the last couple of days and that has caused some
Anyhow, here is how it works for me:
In my AFS home directory
/afs/ovpit.indiana.edu/home/gustavI have created a directory called "jacobi" and made it AFS-open to all:
$ cd /afs/ovpit.indiana.edu/home/gustav $ fs setacl jacobi -acl system:anyuser all $ fs listacl jacobi Access list for jacobi is Normal rights: system:anyuser rlidwka gustav rlidwka $ cd jacobi $ make clean $ make jacobi-raster $ make run-raster poe jacobi-raster -euidevice css0 -procs 8 -euilib ip -ilevel 3 \ > jacobi-raster.log 2> jacobi-raster.err $ cd .. $ fs setacl jacobi -acl system:anyuser none $ fs listacl jacobi Access list for jacobi is Normal rights: gustav rlidwka
Perhaps a better way would be to create a directory called
your AFS home and make it
-acl system:anyuser none, and then create
work and make it
-acl system:anyuser all. I've tested that too and it worked just
fine. That way you
don't have to be so paranoid about running
fs setacl on your
jacobi directory after you have finished your POE run.
Your AFS quota is currently 100MB per student. There is plenty of space on the AFS and that quota can be easily increased if you need more.
Please feel free to see me if you would like me to help you set up your AFS and SP environments.
You need to study the discussion of the FFT subroutine in the PESSL manual
View that manual with
ghostview and go to pages 623-632. Observe that
the Fourier transform of the diffusion coefficient function is going
to be complex in general (in the lecture I have demonstrated that the
F is real and symmetric, but that is not the same
transform of D(x,y) alone).
On the other hand, a Fourier transform of the initial condition T0(x,y)is going to be real (can you explain why?). The sizing condition in this case is somewhat different than in a real-to-complex case. Compare with notes on page 626 of the PESSL manual.