next up previous
Next: Initial Data Up: Converting Existing Programs Previous: Converting Existing Programs

Reid's 3D Boson Stars

 

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.



next up previous
Next: Initial Data Up: Converting Existing Programs Previous: Converting Existing Programs



Robert Marsa
Fri Jul 14 13:58:46 CDT 1995