Typically one would also specify , but we'll ignore this fact for
now. To set this problem up with ` RNPL` we must identify our requirements.
We'll need one grid function . We'll need four difference operators
for , , , and .
We'll also need parameters for the initial data, namely , ,
, **A**, , , and as well as parameters
for the domain boundaries, **xmin**, **xmax**, **ymin**, **ymax**, **zmin**, and
**zmax**.

We begin our ` RNPL` program by specifying the parameters. While ` RNPL` \
statements may appear in any order, it is good for users to keep things
organized. Our parameter declarations look like this:

parameter float xmin := 0 parameter float xmax := 100 parameter float ymin := 0 parameter float ymax := 100 parameter float zmin := 0 parameter float zmax := 100 parameter float A := 1.0 parameter float c_x := 50.0 parameter float c_y := 50 parameter float c_z := 50 parameter float delta_x parameter float delta_y parameter float delta_z

Each parameter can be given a default value. There are other parameters
the ` RNPL` program will need, but they are automatically defined, so we'll
discuss them later.

Next, we define the coordinate system. The declaration looks like this:

rect coordinates t,x,y,z

The name ` rect` is the name of the coordinate system. We can now
define a grid which uses these coordinates.

uniform rect grid g1 [1:Nx][1:Ny][1:Nz] {xmin:xmax}{ymin:ymax}{zmin:zmax}

We will define our grid function to have three time levels so we can use the standard leap-frog operators to solve the equation. The definition is:

float phi on g1 at -1,0,1

Now comes the operator definitions. We need four second derivatives, one for each coordinate.

operator D_LF(f,t,t) := (<1>f[0][0][0] - 2*<0>f[0][0][0] + <-1>f[0][0][0])/(dt*dt) operator D_LF(f,x,x) := (<0>f[1][0][0] - 2*<0>f[0][0][0] + <0>f[-1][0][0])/(dx*dx) operator D_LF(f,y,y) := (<0>f[0][1][0] - 2*<0>f[0][0][0] + <0>f[0][-1][0])/(dy*dy) operator D_LF(f,z,z) := (<0>f[0][0][1] - 2*<0>f[0][0][0] + <0>f[0][0][-1])/(dz*dz)

Since we wish ` RNPL` to produce the complete program, we must specify the
partial-differential equations. This is done by defining the residual.

evaluate residual phi { [1:Nx][1:Ny][1:1] := <1>phi[0][0][0] = 0; [1:Nx][1:Ny][Nz:Nz] := <1>phi[0][0][0] = 0; [1:Nx][1:1][1:Nz] := <1>phi[0][0][0] = 0; [1:Nx][Ny:Ny][1:Nz] := <1>phi[0][0][0] = 0; [1:1][1:Ny][1:Nz] := <1>phi[0][0][0] = 0; [Nx:Nx][1:Ny][1:Nz] := <1>phi[0][0][0] = 0; [2:Nx-1][2:Ny-1][2:Nz-1] := D_LF(phi,t,t) = D_LF(phi,x,x) + D_LF(phi,y,y) + D_LF(phi,z,z) }

The boundary conditions could also have been stated with a time derivative of phi, but this would have required another operator definition. The above method is the simplest.

To get ` RNPL` to generate the initial data, we must provide an
initialization for phi.

initialize phi { [1:Nx][1:Ny][1:Nz] := A*exp(-(x-c_x)^2/delta_x^2)* exp(-(y-c_y)^2/delta_y^2)* exp(-(z-c_z)^2/delta_z^2) }

We now instruct ` RNPL` to solve the equation iteratively and to
automatically generate the update routine.

looper iterative auto update phi

Save the program to a file named ` wave3d_rnpl`. The ` RNPL` compiler
can now produce C or FORTRAN code with one of the following commands.

rnpl -c wave3d_rnpl rnpl -f77 wave3d_rnpl

We can produce a makefile for this program using a utility included with
the ` RNPL` distribution.

Before we can run the program, we need to create a parameter file. Here is an example:

parameters for wave3d tag := "3d_" lambda := .3 Nx0 := 32 Ny0 := 32 Nz0 := 32 xmin := 0 xmax := 10 ymin := 0 ymax := 10 zmin := 0 zmax := 10 ser := 0 fout := 1 iter := 10 epsiter := 1.0e-5 rmod := 1 A := 1.0 c_x := 5.0 c_y := 5 c_z := 5 delta_x := .8 delta_y := .8 delta_z := .8 in_file := "w3d_in.hdf" out_file := "w3d_out.hdf" level:=0

Along with the parameters we defined are some automatically defined
parameters. See the Reference Manual for more information. Save this to a
file called ` w3d_0`.

The program can now be executed with the following command:

wave3d w3d_0

As stated earlier, the initial data should be specified with values for
both and . Since we've only specified , we have
little control over the initial data. ` RNPL` generates a second time level
of data using an iterative procedure over which the user has no control.
So how could we get proper initial data? One was is to write a separate
program for generating initial data from user defined parameters. The only
(current) way for ` RNPL` to generate ``good'' initial data is to reduce the
problem to 1st order form and specify data for the two derivatives of
. This approach is taken in the next section.

Fri Jul 14 13:58:46 CDT 1995