As a first example we'll look briefly at Reid Guenther's 3D boson star code. In order to convert a code we must know what grid functions it uses and what the update routine looks like. In the case of Reid's code we have a complex scalar field, a potential, and some other functions. The declarations are:
# This program solves 3D time dependent schrodinger equation # coupled with a newtonian potential: # phi_t = i/2*(phi_xx + phi_yy + phi_zz) - i*V*phi # V_xx + V_yy + V_zz = phi*conjg(phi) constant parameter float xmin := 0 constant parameter float xmax := 100 constant parameter float ymin := 0 constant parameter float ymax := 100 constant parameter float zmin := 0 constant parameter float zmax := 100 constant parameter int solve_ncycle0 := 10 parameter int solve_ncycle := 3 # these parameters are for determining nu - boundary layer function constant parameter int order := 2 constant parameter float edgdst := 0 constant parameter float delta := 0 constant parameter float height := 0 constant parameter float epsilon := 1.0e-6 rect coordinates t,x,y,z uniform rect grid g1 float phre on g1 at 0,1 float phim on g1 at 0,1 float v on g1 at -1,0,1 alias float rho on g1 float nu on g1 float vnph on g1 looper standard
Parameters that are declared as constant will be passed to update routines in a global common block. Other parameters can be passed in the header. Since the update routine already exists, there is no reason to define operators or residuals. All that is left is the update declaration.
The update routine has the following header:
subroutine evolve3d(phren, phimn, phrenp1, phimnp1,
* vnm1, vn, vnph, vnp1,
* rho, nu, wksp,
* nwksp,
* ds,
* x, y, z,
* dt,
* solve_ncycle)
integer ds(3), solve_ncycle, nwksp
real*8 phren(ds(1),ds(2),ds(3))
real*8 phimn(ds(1),ds(2),ds(3))
real*8 phrenp1(ds(1),ds(2),ds(3))
real*8 phimnp1(ds(1),ds(2),ds(3))
real*8 vnm1(ds(1),ds(2),ds(3))
real*8 vn(ds(1),ds(2),ds(3))
real*8 vnph(ds(1),ds(2),ds(3))
real*8 vnp1(ds(1),ds(2),ds(3))
real*8 rho(ds(1),ds(2),ds(3))
real*8 nu(ds(1),ds(2),ds(3))
real*8 x(ds(1)), y(ds(2)), z(ds(3))
real*8 wksp(nwksp)
real*8 dt
We must get RNPL to generate a header as close to this as possible. The
first thing we do is take the code from the inside of evolve3d
(everything after the header and before the return) and save it to a file
which we'll call evolve3d.inc. The variable wksp is a temporary
work array. RNPL names such arrays work0, work1, ... This
means we must change the name of the work array in the code. This can be
done with a global search and replace (for instance with sed
s/wksp/work0/g evolve3d.inc. This will change nwksp to nwork0
as it should. This is the only modification necessary to the existing
code. Other codes may need other changes. For instance, the shape of the
grid is called ds above but g1_shp by RNPL . If the shape is
needed in the actual code, this name will have to be changed as well. The
RNPL code to generate the header is:
evolve3d.inc evolve3d UPDATES phre, phim, v
HEADER phre[phrenp1,phren], phim[phimnp1,phimn], v[vnp1,vn,vnm1], rho,
vnph,
nu, AUTO work#0(5*nx*ny*nz), x, y, z, solve_ncycle, dt
The time levels that RNPL generates for phre will be phre_n and phre_np1 by default. However, in Reid's code these grid functions are called phren and phrenp1. Likewise for the other grid functions. This naming problem is resolved with the bracket notation above.
This completes the RNPL program. Save this to a file called bos0.rnpl.
When the RNPL compiler is run ( rnpl -f77 bos0.rnpl), it pulls the contents of evolve3d.inc into the stub it produces and calls the routine evolve3d. The result is placed in a file called updates.f. The header in this file is:
!----------------------------------------------------------------------
! This routine updates the following grid functions
! phre phim v
!----------------------------------------------------------------------
subroutine evolve3d(phrenp1,phren,phimnp1,phimn,vnp1,vn,vnm1,rho,
& nu,vnph,g1_shp,g1_bds,x,y,z,dt,solve_ncycle,work0,nwork0)
implicit none
include 'globals.inc'
integer g1_shp(3)
integer g1_bds(6)
real*8 phre_np1(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 phre_n(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 phim_np1(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 phim_n(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 v_np1(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 v_n(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 v_nm1(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 rho(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 nu(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 vnph(g1_bds(1):g1_bds(2),g1_bds(3):g1_bds(4),
& g1_bds(5):g1_bds(6))
real*8 x(*)
real*8 y(*)
real*8 z(*)
real*8 dt
integer solve_ncycle
integer nwork0
real*8 work0(nwork0)
When RNPL is run, it produces the following warnings:
warning: calc_resid: Update exists for phre but no residual. warning: calc_resid: Update exists for phim but no residual. warning: calc_resid: Update exists for v but no residual. warning: calc_resid: No residual is being evaluated. warning: calc_resid: Update exists for phre but no residual. warning: calc_resid: Update exists for phim but no residual. warning: calc_resid: Update exists for v but no residual. warning: calc_resid: No residual is being evaluated.
These would be a problem if we were expecting RNPL to produce its own updates.